ATLAS Offline Software
Egamma1BaselineAlgTool.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2025 CERN for the benefit of the ATLAS collaboration
3 */
4 
6 #include "../dump.h"
7 #include "../dump.icc"
10 
11 namespace GlobalSim {
12 
14  const std::string& name,
15  const IInterface* parent) :
16  base_class(type, name, parent){
17  }
18 
20 
21  CHECK(m_nbhdContainerReadKey.initialize());
24 
25  return StatusCode::SUCCESS;
26  }
27 
29  Egamma1BaselineAlgTool::run(const EventContext& ctx) const {
30  ATH_MSG_DEBUG("run()");
31 
32 
33  // read in LArStrip neighborhoods from the event store
34  // there is one neighborhood per EFex RoI
35  auto in =
37  ctx);
38  CHECK(in.isValid());
39 
40  ATH_MSG_DEBUG("read in " << (*in).size() << " neighborhoods");
41 
42  ap_int<16> peak = 0;
43  ap_int<16> secondMax = 0;
44 
46  CHECK(h_eRatio.record(std::make_unique<std::vector<float> >()));
48  CHECK(h_eRatioSimple.record(std::make_unique<std::vector<float> >()));
49 
50  for (const auto& nbhd : *in) {
51  auto c_phi = combine_phi(nbhd);
52  if (msgLevel() <= MSG::DEBUG) {
53  std::stringstream ss;
54  ss << "Baseline input: ";
55  for (const auto& i : c_phi) {ss << i << ' ';}
56  ATH_MSG_DEBUG(ss.str());
57  }
58  //if (c_phi.empty()) {continue;} // corner case: not all phi have len 17
59  auto input = digitizer::digitize16(c_phi);
60  //Do want to ap_int them?
61  //ap_int<10>* c_input = &input[0]; // vector->array
62 
63  if (msgLevel() <= MSG::DEBUG) {
64  std::stringstream ss;
65  ss << "Baseline input: ";
66  for (const auto& i : input) {ss << i << ' ';}
67  ATH_MSG_DEBUG(ss.str());
68  }
69 
70  //Central peak is easy, it's column 9 in the middle row, i.e. entry 25.
71  peak = input.at(25);
72  ATH_MSG_DEBUG("Peak " << peak);
73 
74  //Find the 6 secondary peaks, lets do it like the VHDL
75  //holder of the current second peak (we will have 6)
76  std::vector<ap_int<16>> secondPeak;
77  secondPeak.resize(6);
78 
79  //Noise Margin: This should be 2sigma... set to 1 for now...
80  int noiseMargin = 1;
81  //Route one: <<< middle row, 24 down to 17
82  secondPeak[0] = secondPeakSearch(input, peak, 24, 17, noiseMargin);
83  //Route two: middle row >>>, 26 up to 33
84  secondPeak[1] = secondPeakSearch(input, peak, 26, 33, noiseMargin);
85  //Route three: <<< top row, 8 down to 0
86  secondPeak[2] = secondPeakSearch(input, peak, 8, 0, noiseMargin);
87  //Route four: top row >>>, 8 up to 16
88  secondPeak[3] = secondPeakSearch(input, peak, 8, 16, noiseMargin);
89  //Route five: <<< bottom row, 42 down to 34
90  secondPeak[4] = secondPeakSearch(input, peak, 42, 34, noiseMargin);
91  //Route five: bottom row >>>, 42 up to 50
92  secondPeak[5] = secondPeakSearch(input, peak, 42, 50, noiseMargin);
93 
94  auto result = std::max_element(secondPeak.begin(), secondPeak.end());
95  ATH_MSG_DEBUG("Max element found at index "
96  << std::distance(secondPeak.begin(), result)
97  << " has value " << *result);
98  //tidy up.
99  secondMax = *result;
100  ATH_MSG_DEBUG("Peak " << peak << " second " << secondMax);
101 
102  //I am not confident on waht div_gen_0 does. So I am casting to floats for now.
103  if(peak > 0 || secondMax > 0){
104  auto eRatio = static_cast< float >(peak - secondMax)/static_cast< float >(peak + secondMax);
105  h_eRatio->push_back(eRatio);
106  ATH_MSG_DEBUG("eRatio (p-sp/p+sp) is " << eRatio);
107 
108  auto eRatioSimple = static_cast< float >(secondMax)/static_cast< float >(peak);
109  h_eRatioSimple->push_back(eRatioSimple);
110  ATH_MSG_DEBUG("eRatio (sp/p) is " << eRatio);
111  }
112 
113  }
114 
115  return StatusCode::SUCCESS;
116  }
117 
119  const ap_int<16> peak,
120  const int startCell,
121  const int endCell,
122  const ap_int<16> noiseMargin) const {
123  //First set a series of counters to "descending" = 0
124  int ascending = 0;
125 
126  //Setup holder of the last energy we checked... which starts at the peak...
127  ap_int<16> lastEnergy = peak;
128  ap_int<16> secondPeak = 0;
129 
130  //There must be a better way to do this?
131  int direction = 0;
132  if(startCell > endCell) {
133  direction =-1;
134  } else {
135  direction =1;
136  }
137 
138  //Route 1: go down in the middlle row
139  for (auto itr = input.begin() + startCell; itr != input.begin() + endCell + direction; itr+=direction){
140  //Going down hill... until we are not...
141  ATH_MSG_DEBUG("Input is " << *itr << " last energy is " << lastEnergy);
142  if(ascending==0 && *itr>lastEnergy && *itr-lastEnergy > noiseMargin){
143  ATH_MSG_DEBUG("We are going up now " << *itr << " is more then " << lastEnergy);
144  ascending=1;
145  lastEnergy=*itr;
146  //Now check if we are at the peak as we are going uphill
147  } else if(ascending==1 && lastEnergy>*itr && lastEnergy-*itr > noiseMargin){
148  ATH_MSG_DEBUG("We are past the top " << *itr << " is less than " << lastEnergy);
149  if(secondPeak<*itr) secondPeak = lastEnergy;
150  ATH_MSG_DEBUG("The peak was " << secondPeak);
151  //I think we can break here... don't need to find a third peak
152  break;
153  } else {
154  lastEnergy=*itr;
155  }
156  }
157 
158  return secondPeak;
159  }
160 
161  std::vector<double>
163  auto result = std::vector<double>();
164 
165  const auto& phi_low = nbhd->phi_low();
166  //if (phi_low.size() != s_required_phi_len) {return result;}
167 
168  const auto& phi_center = nbhd->phi_center();
169  //if (phi_center.size() != s_required_phi_len) {return result;}
170 
171  const auto& phi_high = nbhd->phi_high();
172  //if (phi_high.size() != s_required_phi_len) {return result;}
173 
174  result.reserve(s_combination_len);
175 
176  //Make a big vector as the VHDL does. Not strictly necessary, could just use the nbhd.
177  std::transform(std::begin(phi_high), std::end(phi_high), std::back_inserter(result), [](const auto& high) {
178  return high.m_e;
179  });
180  std::transform(std::begin(phi_center), std::end(phi_center), std::back_inserter(result), [](const auto& center) {
181  return center.m_e;
182  });
183  std::transform(std::begin(phi_low), std::end(phi_low), std::back_inserter(result), [](const auto& low) {
184  return low.m_e;
185  });
186 
187  return result;
188  }
189 
190  std::string Egamma1BaselineAlgTool::toString() const {
191 
192  std::stringstream ss;
193  ss << "Egamma1BaselineAlgTool. name: " << name() << '\n'
194  << m_nbhdContainerReadKey << '\n'
195  << '\n';
196  return ss.str();
197  }
198 }
199 
get_generator_info.result
result
Definition: get_generator_info.py:21
GlobalSim::Egamma1BaselineAlgTool::toString
virtual std::string toString() const override
Definition: Egamma1BaselineAlgTool.cxx:190
PowhegControl_ttHplus_NLO.ss
ss
Definition: PowhegControl_ttHplus_NLO.py:83
GlobalSim::Egamma1BaselineAlgTool::run
virtual StatusCode run(const EventContext &ctx) const override
Definition: Egamma1BaselineAlgTool.cxx:29
make_unique
std::unique_ptr< T > make_unique(Args &&... args)
Definition: SkimmingToolEXOT5.cxx:23
GlobalSim::LArStripNeighborhood::phi_center
const StripDataVector & phi_center() const
Definition: LArStripNeighborhood.h:40
SG::ReadHandle
Definition: StoreGate/StoreGate/ReadHandle.h:67
GlobalSim::Egamma1BaselineAlgTool::m_eRatioKey
SG::WriteHandleKey< std::vector< float > > m_eRatioKey
Definition: Egamma1BaselineAlgTool.h:58
GlobalSim::Egamma1BaselineAlgTool::initialize
StatusCode initialize() override
Definition: Egamma1BaselineAlgTool.cxx:19
PlotCalibFromCool.begin
begin
Definition: PlotCalibFromCool.py:94
GlobalSim::digitizer::digitize16
static std::vector< ap_int< 16 > > digitize16(const std::vector< double > &v)
Definition: Digitizer.h:42
python.CaloAddPedShiftConfig.type
type
Definition: CaloAddPedShiftConfig.py:42
mergePhysValFiles.end
end
Definition: DataQuality/DataQualityUtils/scripts/mergePhysValFiles.py:92
Egamma1BaselineAlgTool.h
GlobalSim
AlgTool to obtain a selection of eFex RoIs read in from the event store.
Definition: dump.h:8
MonitoredCollection.h
GlobalSim::ap_int
Definition: ap_int.h:23
lumiFormat.i
int i
Definition: lumiFormat.py:85
EL::StatusCode
::StatusCode StatusCode
StatusCode definition for legacy code.
Definition: PhysicsAnalysis/D3PDTools/EventLoop/EventLoop/StatusCode.h:22
ATH_MSG_DEBUG
#define ATH_MSG_DEBUG(x)
Definition: AthMsgStreamMacros.h:29
Amg::transform
Amg::Vector3D transform(Amg::Vector3D &v, Amg::Transform3D &tr)
Transform a point from a Trasformation3D.
Definition: GeoPrimitivesHelpers.h:156
test_pyathena.parent
parent
Definition: test_pyathena.py:15
GlobalSim::Egamma1BaselineAlgTool::m_eRatioSimpleKey
SG::WriteHandleKey< std::vector< float > > m_eRatioSimpleKey
Definition: Egamma1BaselineAlgTool.h:64
CHECK
#define CHECK(...)
Evaluate an expression and check for errors.
Definition: Control/AthenaKernel/AthenaKernel/errorcheck.h:422
GlobalSim::Egamma1BaselineAlgTool::m_nbhdContainerReadKey
SG::ReadHandleKey< LArStripNeighborhoodContainer > m_nbhdContainerReadKey
Definition: Egamma1BaselineAlgTool.h:51
python.TriggerAPI.TriggerAPISession.ascending
ascending
Definition: TriggerAPISession.py:435
SG::VarHandleKey::initialize
StatusCode initialize(bool used=true)
If this object is used as a property, then this should be called during the initialize phase.
Definition: AthToolSupport/AsgDataHandles/Root/VarHandleKey.cxx:103
Monitored.h
Header file to be included by clients of the Monitored infrastructure.
GlobalSim::LArStripNeighborhood::phi_low
const StripDataVector & phi_low() const
Definition: LArStripNeighborhood.h:39
GlobalSim::LArStripNeighborhood::phi_high
const StripDataVector & phi_high() const
Definition: LArStripNeighborhood.h:41
name
std::string name
Definition: Control/AthContainers/Root/debug.cxx:240
SG::WriteHandle
Definition: StoreGate/StoreGate/WriteHandle.h:73
SG::WriteHandle::record
StatusCode record(std::unique_ptr< T > data)
Record a const object to the store.
GlobalSim::Egamma1BaselineAlgTool::secondPeakSearch
ap_int< 16 > secondPeakSearch(const std::vector< ap_int< 16 >> &input, const ap_int< 16 > peak, const int startCell, const int endCell, const ap_int< 16 > noiseMargin) const
Definition: Egamma1BaselineAlgTool.cxx:118
DEBUG
#define DEBUG
Definition: page_access.h:11
GlobalSim::Egamma1BaselineAlgTool::s_combination_len
static constexpr int s_combination_len
Definition: Egamma1BaselineAlgTool.h:81
GlobalSim::LArStripNeighborhood
Definition: LArStripNeighborhood.h:28
GlobalSim::Egamma1BaselineAlgTool::combine_phi
std::vector< double > combine_phi(const LArStripNeighborhood *) const
Definition: Egamma1BaselineAlgTool.cxx:162
Amg::distance
float distance(const Amg::Vector3D &p1, const Amg::Vector3D &p2)
calculates the distance between two point in 3D space
Definition: GeoPrimitivesHelpers.h:54
GlobalSim::Egamma1BaselineAlgTool::Egamma1BaselineAlgTool
Egamma1BaselineAlgTool(const std::string &type, const std::string &name, const IInterface *parent)
Definition: Egamma1BaselineAlgTool.cxx:13