ATLAS Offline Software
LArFCS_StepInfoSD.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration
3 */
4 
5 // class header
6 #include "LArFCS_StepInfoSD.h"
8 #include "G4Step.hh"
9 #include "G4ThreeVector.hh"
10 #include "G4TouchableHistory.hh"
11 
13 #include <string>
14 #include <utility>
15 
16 #include "CaloDetDescr/CaloDetDescrElement.h"
18 #include "GaudiKernel/MsgStream.h"
20 
23  : FCS_StepInfoSD(std::move(a_name), config),
24  m_calculator(m_config.m_LArCalculator) {}
25 
27 
28 G4bool LArFCS_StepInfoSD::ProcessHits(G4Step* a_step, G4TouchableHistory*) {
29  G4bool result(false);
30  // If there's no energy, there's no hit. (Aside: Isn't this energy
31  // the same as the energy from the calculator? Not necessarily.
32  // The calculator may include detector effects such as
33  // charge-collection which are not modeled by Geant4.)
34  if (a_step->GetTotalEnergyDeposit() <= 0.) {
35  return result;
36  }
37 
38  if (m_calculator) {
39  if (!m_calo_dd_man.get()) {
41  }
42 
43  const double StepLength = a_step->GetStepLength() / CLHEP::mm;
44  const G4ThreeVector preStepPosition =
45  a_step->GetPreStepPoint()
46  ->GetPosition(); // pre step is the position we're interested in
47  const G4ThreeVector postStepPosition =
48  a_step->GetPostStepPoint()->GetPosition();
49  std::vector<const G4Step*> steps;
50  bool madeSubSteps(false);
52  // create smaller substeps instead
53  G4int nsub_step = (int)(StepLength / m_config.substpsize) + 1;
54  G4double delta = 1. / ((double)nsub_step);
55  // G4cout <<"Orig prestep: "<<preStepPosition<<std::endl;
56  for (G4int i = 0; i < nsub_step; i++) {
57  G4double fraction1 = ((G4double)i) * delta;
58  G4double fraction2 = (((G4double)i) + 1.) * delta;
59  G4ThreeVector subpoint1 =
60  preStepPosition * (1 - fraction1) + postStepPosition * fraction1;
61  G4ThreeVector subpoint2 =
62  preStepPosition * (1 - fraction2) + postStepPosition * fraction2;
63  G4StepPoint* startpoint = new G4StepPoint(*(a_step->GetPreStepPoint()));
64  G4StepPoint* endpoint = new G4StepPoint(*(a_step->GetPostStepPoint()));
65  startpoint->SetPosition(subpoint1);
66  endpoint->SetPosition(subpoint2);
67 
68  G4Step* newstep = new G4Step(*a_step);
69  newstep->SetPreStepPoint(startpoint);
70  newstep->SetPostStepPoint(endpoint);
71  newstep->SetStepLength((subpoint1 - subpoint2).mag());
72  newstep->SetTotalEnergyDeposit(a_step->GetTotalEnergyDeposit() /
73  nsub_step);
74  steps.push_back(newstep);
75  madeSubSteps = true;
76  }
77  } else {
78  steps.push_back(a_step);
79  }
80 
81  if (m_config.verboseLevel > 4) {
82  G4cout << "LArFCS_StepInfoSD::ProcessHits(): original step in Volume: "
83  << a_step->GetPreStepPoint()->GetPhysicalVolume()->GetName()
84  << " position: " << a_step->GetPreStepPoint()->GetPosition()
85  << " Length=" << a_step->GetStepLength() / CLHEP::mm
86  << " E=" << a_step->GetTotalEnergyDeposit() << G4endl;
87  std::vector<LArHitData> processedHits;
88  if (m_calculator->Process(a_step, processedHits)) {
89  if (m_config.verboseLevel > -1) {
90  G4cout << " #LArHitData=" << processedHits.size();
91  for (const auto& lhd : processedHits) {
92  G4cout << " ; id=" << (std::string)lhd.id << " E=" << lhd.energy
93  << G4endl;
94  }
95  G4cout << G4endl;
96  }
97  }
98  G4cout << "LArFCS_StepInfoSD::ProcessHits(): #substep=" << steps.size()
99  << G4endl;
100  }
101 
102  double et_all = 0; // Total collected charge in all substeps
103  for (const G4Step* substep : steps) {
104  double et(0.); // Total collected charge in this substep
105  G4ThreeVector stepPosition =
106  0.5 * (substep->GetPreStepPoint()->GetPosition() +
107  substep->GetPostStepPoint()->GetPosition());
108  std::vector<LArHitData> processedHits;
109  if (m_config.verboseLevel > 4) {
110  G4cout << " LArFCS_StepInfoSD::ProcessHits(): substep in Volume: "
111  << substep->GetPreStepPoint()->GetPhysicalVolume()->GetName()
112  << " position: " << substep->GetPreStepPoint()->GetPosition()
113  << " Length=" << substep->GetStepLength() / CLHEP::mm
114  << " E=" << substep->GetTotalEnergyDeposit() << G4endl;
115  }
116  if (m_calculator->Process(substep, processedHits)) {
117  for (const auto& larhit : processedHits) {
118  et += larhit.energy;
119  et_all += larhit.energy;
120  }
121  if (m_config.verboseLevel > 4) {
122  G4cout << " substep #LArHitData=" << processedHits.size();
123  for (const auto& lhd : processedHits) {
124  G4cout << " ; id=" << (std::string)lhd.id << " E=" << lhd.energy
125  << G4endl;
126  }
127  G4cout << G4endl;
128  }
129  } else {
130  if (m_config.verboseLevel > 4) {
131  // Maybe 0 hits or something like that...
132  G4cout << this->GetName()
133  << " WARNING ProcessHits: Call to ILArCalculatorSvc::Process "
134  "failed! Details:"
135  << G4endl << " " << "Volume: "
136  << a_step->GetPreStepPoint()->GetPhysicalVolume()->GetName()
137  << " " << m_calculator << " position: " << stepPosition
138  << " SL: " << StepLength << G4endl << " "
139  << "Orig position: "
140  << substep->GetPreStepPoint()->GetPosition() << " / "
141  << substep->GetPostStepPoint()->GetPosition() << G4endl
142  << " " << "StepLength: " << StepLength
143  << " step: " << preStepPosition << " / " << postStepPosition
144  << G4endl;
145  }
146  continue;
147  }
148 
149  // drop hits with zero deposited energy (could still happen with negative
150  // corrections from calculator)
151  // Or if total energy is <0
152  if (et <= 0.) {
153  if (m_config.verboseLevel > 4) {
154  G4cout << this->GetName()
155  << " WARNING ProcessHits: Total negative energy: " << et
156  << " not processing..." << G4endl;
157  }
158  continue;
159  }
160 
161  const size_t numberOfProcessedHits = processedHits.size();
162  const G4ThreeVector originalStepPosition = stepPosition;
163  double maxSubHitEnergy = -999.;
164  int maxSubHitEnergyindex = -1;
165  if (numberOfProcessedHits > 0) {
166  maxSubHitEnergy = processedHits[0].energy;
167  maxSubHitEnergyindex = 0;
168  }
169  // Figure out the subhit with most energy
170  for (size_t i = 1; i < numberOfProcessedHits; ++i) {
171  if (maxSubHitEnergy < processedHits[i].energy) {
172  maxSubHitEnergy = processedHits[i].energy;
173  maxSubHitEnergyindex = i;
174  }
175  }
176  if (maxSubHitEnergyindex ==
177  -1) { // because there were no hits; numberOfProcessedHits ==0
178  G4cout << this->GetName()
179  << " WARNING ProcessHits: numberOfProcessedHits is zero"
180  << G4endl;
181  continue;
182  }
183  // Identifier for the subhit with max energy
184  Identifier maxEnergyIdentifier =
185  this->ConvertID(processedHits[maxSubHitEnergyindex].id);
186  const CaloDetDescrElement* maxEnergyCell =
187  m_calo_dd_man.get()->get_element(maxEnergyIdentifier);
188 
189  Identifier invalidIdentifier;
190  // for (size_t i=0; i<numberOfProcessedHits; ++i) {
191  for (const auto& larhit : processedHits) {
192  Identifier id = this->ConvertID(larhit.id);
193  if (larhit.id[0] == 10) {
194  if (m_config.verboseLevel > 9) {
195  G4cout << this->GetName()
196  << " VERBOSE ProcessHits: Dead Material LArG4Identifier: "
197  << (std::string)larhit.id << G4endl << " " << id
198  << G4endl << " " << id.getString() << G4endl
199  << " "
200  << a_step->GetPreStepPoint()
201  ->GetPhysicalVolume()
202  ->GetLogicalVolume()
203  ->GetName()
204  << G4endl << " numberOfProcessedHits: "
205  << numberOfProcessedHits << G4endl;
206  G4cout << " " << invalidIdentifier << G4endl;
207  }
208  } else if (id == invalidIdentifier) {
209  G4cout
210  << this->GetName()
211  << " WARNING ProcessHits: Something wrong with LArG4Identifier: "
212  << (std::string)larhit.id << G4endl << " " << id
213  << G4endl << " " << id.getString() << G4endl
214  << " "
215  << a_step->GetPreStepPoint()
216  ->GetPhysicalVolume()
217  ->GetLogicalVolume()
218  ->GetName()
219  << G4endl
220  << " numberOfProcessedHits: " << numberOfProcessedHits
221  << G4endl;
222  G4cout << " " << invalidIdentifier << G4endl;
223  }
224  // need to get the cell information
225  if (numberOfProcessedHits > 1) {
226  if (!m_larEmID->is_em_barrel(id)) {
227  // It didn't seem to happen outside em_barrel, so flag up if it
228  // does:
229  G4cout << this->GetName()
230  << " WARNING ProcessHits: Outside LAr barrel, but "
231  "numberOfProcessedHits="
232  << numberOfProcessedHits
233  << ", LArG4Identifier: " << (std::string)larhit.id << G4endl;
234  } else {
236  // find subhit with largest energy
237  if (maxSubHitEnergyindex == -1) {
238  G4cout
239  << this->GetName()
240  << " WARNING ProcessHits: no subhit index with e>-999??? "
241  << G4endl;
242  continue;
243  }
244  if (m_config.verboseLevel > 9) {
245  G4cout << this->GetName()
246  << " VERBOSE ProcessHits: shifting subhits: largest "
247  "energy subhit index is "
248  << maxSubHitEnergyindex << " E: " << maxSubHitEnergy
249  << " identifier: " << maxEnergyIdentifier.getString()
250  << G4endl;
251  }
252  // from identifier
253  const CaloDetDescrElement* thiscell =
254  m_calo_dd_man.get()->get_element(id);
255  if (!maxEnergyCell) {
256  // How often does this happen? Do not shift.
257  G4cout << this->GetName()
258  << " WARNING ProcessHits: maxEnergyCell failed: "
259  << maxEnergyIdentifier.getString() << G4endl
260  << " "
261  << m_calo_dd_man.get()->get_element(id)->getSampling()
262  << G4endl << " " << originalStepPosition.eta()
263  << G4endl << " " << originalStepPosition.phi()
264  << G4endl;
265  stepPosition = originalStepPosition;
266  } else if (maxEnergyCell == thiscell) {
267  // The cells match, so do not shift this hit.
268  if (m_config.verboseLevel > 9) {
269  G4cout << this->GetName()
270  << " VERBOSE ProcessHits: Original step position: "
271  << originalStepPosition.x() << " "
272  << originalStepPosition.y() << " "
273  << originalStepPosition.z() << G4endl << " "
274  << "This cell: " << thiscell->x() << " "
275  << thiscell->y() << " " << thiscell->z() << G4endl
276  << " " << "No shift" << G4endl;
277  }
278  stepPosition = originalStepPosition;
279  } else {
280  // the two cells do not match => shift
281  G4ThreeVector diff(thiscell->x() - maxEnergyCell->x(),
282  thiscell->y() - maxEnergyCell->y(),
283  thiscell->z() - maxEnergyCell->z());
284  stepPosition = originalStepPosition + diff;
285  if (m_config.verboseLevel > 9) {
286  const CaloDetDescrElement* bestcell =
287  m_calo_dd_man.get()->get_element(
288  m_calo_dd_man.get()->get_element(id)->getSampling(),
289  originalStepPosition.eta(),
290  originalStepPosition.phi());
291  G4cout << this->GetName()
292  << " VERBOSE ProcessHits: Original step position: "
293  << originalStepPosition.x() << " "
294  << originalStepPosition.y() << " "
295  << originalStepPosition.z() << G4endl << " "
296  << "This cell: " << thiscell->x() << " "
297  << thiscell->y() << " " << thiscell->z() << G4endl
298  << " "
299  << "Highest E cell: " << maxEnergyCell->x() << " "
300  << maxEnergyCell->y() << " " << maxEnergyCell->z()
301  << G4endl << " "
302  << "(Best cell: " << bestcell->x() << " "
303  << bestcell->y() << " " << bestcell->z() << ")"
304  << G4endl << " "
305  << "Shifted step position: " << stepPosition.x() << " "
306  << stepPosition.y() << " " << stepPosition.z()
307  << G4endl;
308  }
309  }
310  }
311  }
312  }
313  // Finalize time for LAr hits?: NO
314  // double time = larhit.energy)==0 ? 0. : (double)
315  // larhit.time/larhit.energy/CLHEP::ns;
316  double time = larhit.time;
317  double energy = larhit.energy / CLHEP::MeV;
318  // Drop any hits that don't have a good identifier attached
319  if (!m_calo_dd_man.get()->get_element(id)) {
320  if (m_config.verboseLevel > 4) {
321  G4cout << this->GetName()
322  << " DEBUG update_map: bad identifier: " << id.getString()
323  << " skipping this hit." << G4endl;
324  }
325  continue;
326  }
327  this->update_map(
328  stepPosition, id, energy, time, true,
329  numberOfProcessedHits); // store numberOfProcessedHits as info
330  } // numberOfProcessedHits
331  } // istep
332  if (madeSubSteps) {
333  // only delete steps when doing substeps. Do not delete the original
334  // G4Step!
335  while (!steps.empty()) {
336  delete steps.back();
337  steps.pop_back();
338  }
339  }
340  if (m_config.verboseLevel > 4) {
341  G4cout << "LArFCS_StepInfoSD::ProcessHits(): Etotal substeps=" << et_all
342  << G4endl << G4endl << G4endl;
343  }
344  }
345 
346  return result;
347 }
348 
350  Identifier id;
351  if (a_ident[0] == 4) {
352  // is LAr
353  if (a_ident[1] == 1) {
354  // is LAr EM
355  try {
356  id = m_larEmID->channel_id(a_ident[2], // barrel_ec
357  a_ident[3], // sampling
358  a_ident[4], // region
359  a_ident[5], // eta
360  a_ident[6]); // phi
361  } catch (LArID_Exception& e) {
362  G4ExceptionDescription description;
363  description << "ConvertID: LArEM_ID error code " << e.code() << " "
364  << (std::string)e;
365  G4Exception("LArFCS_StepInfoSD", "ConvertIDEM", FatalException,
366  description);
367  abort();
368  }
369  } else if (a_ident[1] == 2) {
370  // is EM HEC
371  try {
372  id = m_larHecID->channel_id(a_ident[2], // zSide
373  a_ident[3], // sampling
374  a_ident[4], // region
375  a_ident[5], // eta
376  a_ident[6]); // phi
377  } catch (LArID_Exception& e) {
378  G4ExceptionDescription description;
379  description << "ConvertID: LArHEC_ID error code " << e.code() << " "
380  << (std::string)e;
381  G4Exception("LArFCS_StepInfoSD", "ConvertIDHEC", FatalException,
382  description);
383  abort();
384  }
385  } else if (a_ident[1] == 3) {
386  // FCAL
387  if (a_ident[3] > 0) {
388  // is EM FCAL
389  try {
390  id = m_larFcalID->channel_id(a_ident[2], // zSide
391  a_ident[3], // sampling
392  a_ident[4], // eta
393  a_ident[5]); // phi
394  } catch (LArID_Exception& e) {
395  G4ExceptionDescription description;
396  description << "ConvertID: LArFCAL_ID error code " << e.code() << " "
397  << (std::string)e;
398  G4Exception("LArFCS_StepInfoSD", "ConvertIDFCAL", FatalException,
399  description);
400  abort();
401  }
402  } else {
403  G4ExceptionDescription description;
404  description << "ConvertID: Unsupported ID. ";
405  G4Exception("LArFCS_StepInfoSD", "ConvertIDMiniFCAL", FatalException,
406  description);
407  abort();
408  }
409  }
410  }
411  return id;
412 }
AllowedVariables::e
e
Definition: AsgElectronSelectorTool.cxx:37
LArG4Identifier
Definition: LArG4Identifier.h:121
et
Extra patterns decribing particle interation process.
get_generator_info.result
result
Definition: get_generator_info.py:21
python.SystemOfUnits.mm
float mm
Definition: SystemOfUnits.py:98
LArFCS_StepInfoSD::m_calculator
ILArCalculatorSvc * m_calculator
Member variable - the calculator we'll use.
Definition: LArFCS_StepInfoSD.h:44
CaloDetDescrElement::y
float y() const
cell y
Definition: Calorimeter/CaloDetDescr/CaloDetDescr/CaloDetDescrElement.h:365
xAOD::et
et
Definition: TrigEMCluster_v1.cxx:25
LArFCS_StepInfoSD::ConvertID
Identifier ConvertID(const LArG4Identifier &a_ident) const
Helper function for making "real" identifiers from LArG4Identifiers.
Definition: LArFCS_StepInfoSD.cxx:349
FCS_StepInfoSD::getCaloDDManager
void getCaloDDManager()
Keep a map instead of trying to keep the full vector.
Definition: FCS_StepInfoSD.cxx:49
ILArCalculatorSvc::Process
virtual G4bool Process(const G4Step *, std::vector< LArHitData > &) const =0
CaloDetDescrElement
This class groups all DetDescr information related to a CaloCell. Provides a generic interface for al...
Definition: Calorimeter/CaloDetDescr/CaloDetDescr/CaloDetDescrElement.h:66
FCS_StepInfoSD::update_map
void update_map(const CLHEP::Hep3Vector &l_vec, const Identifier &l_identifier, double l_energy, double l_time, bool l_valid, int l_detector)
Definition: FCS_StepInfoSD.cxx:62
mc.diff
diff
Definition: mc.SFGenPy8_MuMu_DD.py:14
FCS_Param::Config::shift_lar_subhit
bool shift_lar_subhit
Definition: FCS_StepInfoSD.h:40
FCS_StepInfoSD::m_calo_dd_man
CxxUtils::CachedPointer< const CaloDetDescrManager > m_calo_dd_man
Definition: FCS_StepInfoSD.h:106
LArEM_Base_ID::channel_id
Identifier channel_id(const ExpandedIdentifier &exp_id) const
Build a cell identifier from an expanded identifier.
FCS_Param::Config::substpsize
double substpsize
Definition: FCS_StepInfoSD.h:42
CaloDetDescrManager.h
Definition of CaloDetDescrManager.
config
Definition: PhysicsAnalysis/AnalysisCommon/AssociationUtils/python/config.py:1
python.SystemOfUnits.MeV
float MeV
Definition: SystemOfUnits.py:172
LArFCS_StepInfoSD.h
beamspotman.steps
int steps
Definition: beamspotman.py:503
ParticleGun_FastCalo_ChargeFlip_Config.energy
energy
Definition: ParticleGun_FastCalo_ChargeFlip_Config.py:78
lumiFormat.i
int i
Definition: lumiFormat.py:85
LArFCAL_Base_ID::channel_id
Identifier channel_id(const ExpandedIdentifier &exp_id) const
cell identifier for a channel from ExpandedIdentifier
FCS_StepInfoSD::m_larHecID
const LArHEC_ID * m_larHecID
Definition: FCS_StepInfoSD.h:104
LArFCS_StepInfoSD::LArFCS_StepInfoSD
LArFCS_StepInfoSD(G4String a_name, const FCS_Param::Config &config)
Constructor.
Definition: LArFCS_StepInfoSD.cxx:21
xAOD::double
double
Definition: CompositeParticle_v1.cxx:159
FCS_StepInfoSD::m_larEmID
const LArEM_ID * m_larEmID
Pointers to the identifier helpers.
Definition: FCS_StepInfoSD.h:102
id
SG::auxid_t id
Definition: Control/AthContainers/Root/debug.cxx:239
LArHEC_Base_ID::channel_id
Identifier channel_id(const ExpandedIdentifier &exp_id) const
channel identifier for a channel from ExpandedIdentifier
CaloDetDescrElement::x
float x() const
cell x
Definition: Calorimeter/CaloDetDescr/CaloDetDescr/CaloDetDescrElement.h:363
Identifier::getString
std::string getString() const
Provide a string form of the identifier - hexadecimal.
Definition: Identifier.cxx:25
ILArCalculatorSvc.h
python.CaloAddPedShiftConfig.int
int
Definition: CaloAddPedShiftConfig.py:45
CaloSwCorrections.time
def time(flags, cells_name, *args, **kw)
Definition: CaloSwCorrections.py:242
LArEM_Base_ID::is_em_barrel
bool is_em_barrel(const Identifier id) const
test if the id belongs to the EM barrel
CaloDetDescrElement::z
float z() const
cell z
Definition: Calorimeter/CaloDetDescr/CaloDetDescr/CaloDetDescrElement.h:367
CxxUtils::CachedPointer::get
pointer_t get() const
Return the current value of the element.
FCS_StepInfoSD
Common sensitive detector class for LAr systems.
Definition: FCS_StepInfoSD.h:71
FCS_StepInfoSD::m_config
FCS_Param::Config m_config
Definition: FCS_StepInfoSD.h:100
LArFCS_StepInfoSD::ProcessHits
G4bool ProcessHits(G4Step *a_step, G4TouchableHistory *) override
Main processing method.
Definition: LArFCS_StepInfoSD.cxx:28
FCS_Param::Config
Definition: FCS_StepInfoSD.h:37
Trk::StepLength
@ StepLength
Definition: MaterialAssociationType.h:17
FCS_Param::Config::shorten_lar_step
bool shorten_lar_step
Definition: FCS_StepInfoSD.h:41
FCS_Param::Config::verboseLevel
int verboseLevel
Helper to keep the same verbosity everywhere.
Definition: FCS_StepInfoSD.h:39
G4StepHelper::preStepPosition
G4ThreeVector preStepPosition(const G4Step *theStep)
TODO.
Definition: StepHelper.cxx:14
mag
Scalar mag() const
mag method
Definition: AmgMatrixBasePlugin.h:26
LArID_Exception
Exception class for LAr Identifiers.
Definition: LArID_Exception.h:20
G4StepHelper::postStepPosition
G4ThreeVector postStepPosition(const G4Step *theStep)
TODO.
Definition: StepHelper.cxx:19
description
std::string description
glabal timer - how long have I taken so far?
Definition: hcg.cxx:88
LArFCS_StepInfoSD::~LArFCS_StepInfoSD
virtual ~LArFCS_StepInfoSD()
Destructor.
Definition: LArFCS_StepInfoSD.cxx:26
FCS_StepInfoSD::m_larFcalID
const LArFCAL_ID * m_larFcalID
Definition: FCS_StepInfoSD.h:103
Identifier
Definition: IdentifierFieldParser.cxx:14