PFChargedFlowElementCreatorAlgorithm.cxx 8.52 KB
Newer Older
1
//eflowRec includes
2
#include "eflowRec/PFChargedFlowElementCreatorAlgorithm.h"
3
#include "eflowRec/eflowRecCluster.h"
4
#include "eflowRec/eflowRecTrack.h"
5
#include "eflowRec/eflowTrackClusterLink.h"
6

7
8
//EDM includes
#include "xAODBase/IParticleContainer.h"
9
#include "xAODCaloEvent/CaloClusterContainer.h"
10
#include "xAODPFlow/FlowElementAuxContainer.h"
11
12
#include "xAODPFlow/PFODefs.h"
#include "xAODCore/AuxStoreAccessorMacros.h"
13

14
15
16
17
18
19
20
21
22
23
24
25
26
StatusCode PFChargedFlowElementCreatorAlgorithm::initialize(){

  ATH_CHECK(m_eflowCaloObjectContainerReadHandleKey.initialize());

  ATH_CHECK(m_chargedFlowElementContainerWriteHandleKey.initialize());
  
  return StatusCode::SUCCESS;
}

StatusCode PFChargedFlowElementCreatorAlgorithm::execute(const EventContext& ctx) const {

  ATH_MSG_DEBUG("Starting PFOChargedCreatorAlgorithm::execute");

27
28
29
30
31
  SG::WriteHandle<xAOD::FlowElementContainer> chargedFlowElementContainerWriteHandle(m_chargedFlowElementContainerWriteHandleKey,ctx);
  ATH_CHECK(chargedFlowElementContainerWriteHandle.record(std::make_unique<xAOD::FlowElementContainer>(),std::make_unique<xAOD::FlowElementAuxContainer>()));
  
  /* Create Charged FlowElements from all eflowCaloObjects */
  SG::ReadHandle<eflowCaloObjectContainer> eflowCaloObjectContainerReadHandle(m_eflowCaloObjectContainerReadHandleKey,ctx);
scott snyder's avatar
scott snyder committed
32
  for (const auto thisEflowCaloObject : *eflowCaloObjectContainerReadHandle) createChargedFlowElements(*thisEflowCaloObject,true,chargedFlowElementContainerWriteHandle);
33
34
35

  std::sort(chargedFlowElementContainerWriteHandle->begin(), chargedFlowElementContainerWriteHandle->end(), [] (const xAOD::FlowElement* flowElement1, const xAOD::FlowElement* flowElement2) {return flowElement1->pt()>flowElement2->pt();});

36
37
38
39
40
41
  return StatusCode::SUCCESS;  

}

void PFChargedFlowElementCreatorAlgorithm::createChargedFlowElements(const eflowCaloObject& energyFlowCaloObject, bool addClusters, SG::WriteHandle<xAOD::FlowElementContainer>& chargedFlowElementContainerWriteHandle) const {

42
  /* Loop over all tracks in the eflowCaloObject */
43
  for (unsigned int iTrack = 0; iTrack < energyFlowCaloObject.nTracks(); ++iTrack) {
44
45
46
47

    eflowRecTrack* efRecTrack = energyFlowCaloObject.efRecTrack(iTrack);

    /* Skip tracks that haven't been subtracted */
48
    if (!m_eOverPMode && !efRecTrack->isSubtracted()){continue;}
49

50
    /* Create new xAOD::FlowElement and set the type to charged PFlow */
51
52
    xAOD::FlowElement* thisFE = new xAOD::FlowElement();
    chargedFlowElementContainerWriteHandle->push_back(thisFE);
53
    thisFE->setSignalType(xAOD::FlowElement::SignalType::ChargedPFlow);
54
55
56
57
58
59
60
61
62
63
64
65

    /* Get the track elementLink and add it to the xAOD:FE. 
    Note we first have to convert it to an IParticle ElementLink. */
    ElementLink<xAOD::TrackParticleContainer> theTrackLink = efRecTrack->getTrackElemLink();
    ElementLink< xAOD::IParticleContainer > theIParticleTrackLink; 
    theIParticleTrackLink.resetWithKeyAndIndex(theTrackLink.persKey(),theTrackLink.persIndex() ); 
    std::vector<ElementLink<xAOD::IParticleContainer> > vecIParticleTrackLinkContainer;
    vecIParticleTrackLinkContainer.push_back(theIParticleTrackLink);
    thisFE->setChargedObjectLinks(vecIParticleTrackLinkContainer);

    //Now set the charge
    thisFE->setCharge(efRecTrack->getTrack()->charge());
66
    thisFE->setSignalType(xAOD::FlowElement::ChargedPFlow);
67

68
69
70
71
72
73
    std::pair<double,double> etaPhi(0.0,0.0);

    if (m_eOverPMode){
      /* In EOverPMode want charged eflowObjects to have extrapolated eta,phi as coordinates
       * (needed for analysis of EOverP Data) */
      etaPhi = efRecTrack->getTrackCaloPoints().getEM2etaPhi();
74
75

      /*add information to xAOD*/
Mark Hodgkinson's avatar
Mark Hodgkinson committed
76
      const static SG::AuxElement::Accessor<int> accLHED("LayerHED");
77
78
      accLHED(*thisFE) = efRecTrack->getLayerHED();

Mark Hodgkinson's avatar
Mark Hodgkinson committed
79
      const static SG::AuxElement::Accessor<std::vector<int> > accCellOrderVector("LayerVectorCellOrdering");
80
81
      accCellOrderVector(*thisFE) = efRecTrack->getLayerCellOrderVector();

Mark Hodgkinson's avatar
Mark Hodgkinson committed
82
      const static SG::AuxElement::Accessor<std::vector<float> > accRadiusCellOrderVector("RadiusVectorCellOrdering");
83
84
      accRadiusCellOrderVector(*thisFE) = efRecTrack->getRadiusCellOrderVector();

Mark Hodgkinson's avatar
Mark Hodgkinson committed
85
      const static SG::AuxElement::Accessor<std::vector<float> > accAvgEDensityCellOrderVector("AvgEdensityVectorCellOrdering");
86
      accAvgEDensityCellOrderVector(*thisFE) = efRecTrack->getAvgEDensityCellOrderVector();
87
88
89
90
91
92
93
94
95
    } else {
      /* In normal mode we want the track eta,phi at the perigee */
      etaPhi.first = efRecTrack->getTrack()->eta();
      etaPhi.second = efRecTrack->getTrack()->phi();
    }

    /* Set the 4-vector of the xAOD::PFO */
    thisFE->setP4(efRecTrack->getTrack()->pt(), etaPhi.first, etaPhi.second, efRecTrack->getTrack()->m());

96
97
98
    ATH_MSG_DEBUG("Created charged PFO with E, pt, eta and phi of " << thisFE->e() << ", " << thisFE->pt() << ", " << thisFE->eta() << " and " << thisFE->phi());

    /* Add the amount of energy the track was expected to deposit in the calorimeter - this is needed to calculate the charged weight in the jet finding */
Mark Hodgkinson's avatar
Mark Hodgkinson committed
99
    const static SG::AuxElement::Accessor<float> accTracksExpectedEnergyDeposit("TracksExpectedEnergyDeposit");
100
101
102
103
104
    accTracksExpectedEnergyDeposit(*thisFE) = efRecTrack->getEExpect();
    ATH_MSG_DEBUG("Have set that PFO's expected energy deposit to be " << efRecTrack->getEExpect());

    /* Flag if this track was in a dense environment for later checking */
    //There is an issue using bools - when written to disk they convert to chars. So lets store the bool as an int.
Mark Hodgkinson's avatar
Mark Hodgkinson committed
105
    const static SG::AuxElement::Accessor<int> accIsInDenseEnvironment("IsInDenseEnvironment");
106
107
108
    accIsInDenseEnvironment(*thisFE) = efRecTrack->isInDenseEnvironment();

     /* Optionally we add the links to clusters to the xAOD::PFO */
109
    if (addClusters){
110
111
112
113
114
115
116
117
118
119
120
121
122
123

      std::vector<std::pair<eflowTrackClusterLink*,float> > trackClusterLinkPairs = energyFlowCaloObject.efRecLink();

      std::vector<eflowTrackClusterLink*> thisTracks_trackClusterLinks = efRecTrack->getClusterMatches();

      /** Each eflowCaloObject has a list of clusters for all the tracks it represents.
       *  We only want the subset of the clusters matched to this track, and collect these in thisTracks_trackClusterLinksSubtracted.
       */

      std::vector<eflowTrackClusterLink*> thisTracks_trackClusterLinksSubtracted;

      //Create vector of pairs which map each CaloCluster to the ratio of its new energy to unstracted energy
      std::vector<std::pair<ElementLink<xAOD::CaloClusterContainer>, double> > vectorClusterToSubtractedEnergies;

124
125
      for (auto& trackClusterLink : thisTracks_trackClusterLinks){
        for (auto& trackClusterLinkPair : trackClusterLinkPairs){
126
127
128
129
130
131
132
133
          if (!m_eOverPMode && trackClusterLinkPair.first == trackClusterLink && !std::isnan(trackClusterLinkPair.second)) {
            thisTracks_trackClusterLinksSubtracted.push_back(trackClusterLink);
            eflowRecCluster* efRecCluster = trackClusterLinkPair.first->getCluster();
            ElementLink<xAOD::CaloClusterContainer> theOriginalClusterLink = efRecCluster->getOriginalClusElementLink();
              ElementLink<xAOD::CaloClusterContainer> theSisterClusterLink = (*theOriginalClusterLink)->getSisterClusterLink();
            if (theSisterClusterLink.isValid()) vectorClusterToSubtractedEnergies.push_back(std::pair(theSisterClusterLink,trackClusterLinkPair.second));
            else vectorClusterToSubtractedEnergies.push_back(std::pair(theOriginalClusterLink,trackClusterLinkPair.second));
          }
134
          else if (m_eOverPMode && trackClusterLinkPair.first == trackClusterLink) thisTracks_trackClusterLinksSubtracted.push_back(trackClusterLink);
135
136
137
138
139
140
        }
      }

      //sort the vectorClusterToSubtractedEnergies in order of subtracted energy ratio from low (most subtracted) to high (least subtracted)
      std::sort(vectorClusterToSubtractedEnergies.begin(),vectorClusterToSubtractedEnergies.end(), [](auto const& a, auto const&b){return a.second < b.second;});
      //now split this into two vectors, ready to be used by the FlowElement
141
142
143
144
145
146
147
148
      std::vector<ElementLink<xAOD::IParticleContainer> > theClusters;
      std::vector<float> theClusterWeights;
      for (auto thePair : vectorClusterToSubtractedEnergies){
        ElementLink< xAOD::IParticleContainer > theIParticleTrackLink; 
        theIParticleTrackLink.resetWithKeyAndIndex(thePair.first.persKey(),thePair.first.persIndex()); 
        theClusters.push_back(theIParticleTrackLink);
        theClusterWeights.push_back(thePair.second);
      }
149

150
151
152
153
      //Set the vector of clusters and vector of corresponding energies.
      thisFE->setOtherObjectLinks(theClusters,theClusterWeights);

    }//if we add clusters
154

155
156
  }//loop over eflowRecTracks

157
}