Main processing method.
28 {
30
31
32
33
34 if (a_step->GetTotalEnergyDeposit() <= 0.) {
36 }
37
41 }
42
43 const double StepLength = a_step->GetStepLength() / CLHEP::mm;
45 a_step->GetPreStepPoint()
46 ->GetPosition();
48 a_step->GetPostStepPoint()->GetPosition();
49 std::vector<const G4Step*>
steps;
50 bool madeSubSteps(false);
52
53 G4int nsub_step = (
int)(StepLength /
m_config.substpsize) + 1;
54 G4double delta = 1. / ((
double)nsub_step);
55
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 =
61 G4ThreeVector subpoint2 =
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
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;
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;
103 for (const G4Step* substep : steps) {
105 G4ThreeVector stepPosition =
106 0.5 * (substep->GetPreStepPoint()->GetPosition() +
107 substep->GetPostStepPoint()->GetPosition());
108 std::vector<LArHitData> processedHits;
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 }
117 for (const auto& larhit : processedHits) {
119 et_all += larhit.energy;
120 }
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 {
131
132 G4cout << this->GetName()
133 << " WARNING ProcessHits: Call to ILArCalculatorSvc::Process "
134 "failed! Details:"
135 << G4endl << " " << "Volume: "
136 << a_step->GetPreStepPoint()->GetPhysicalVolume()->GetName()
139 << "Orig position: "
140 << substep->GetPreStepPoint()->GetPosition() << " / "
141 << substep->GetPostStepPoint()->GetPosition() << G4endl
144 << G4endl;
145 }
146 continue;
147 }
148
149
150
151
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
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) {
178 G4cout << this->GetName()
179 << " WARNING ProcessHits: numberOfProcessedHits is zero"
180 << G4endl;
181 continue;
182 }
183
184 Identifier maxEnergyIdentifier =
185 this->
ConvertID(processedHits[maxSubHitEnergyindex].
id);
186 const CaloDetDescrElement* maxEnergyCell =
188
189 Identifier invalidIdentifier;
190
191 for (const auto& larhit : processedHits) {
192 Identifier
id = this->
ConvertID(larhit.id);
193 if (larhit.id[0] == 10) {
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
225 if (numberOfProcessedHits > 1) {
227
228
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
237 if (maxSubHitEnergyindex == -1) {
238 G4cout
239 << this->GetName()
240 << " WARNING ProcessHits: no subhit index with e>-999??? "
241 << G4endl;
242 continue;
243 }
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
253 const CaloDetDescrElement* thiscell =
255 if (!maxEnergyCell) {
256
257 G4cout << this->GetName()
258 << " WARNING ProcessHits: maxEnergyCell failed: "
259 << maxEnergyIdentifier.
getString() << G4endl
260 << " "
262 << G4endl << " " << originalStepPosition.eta()
263 << G4endl << " " << originalStepPosition.phi()
264 << G4endl;
265 stepPosition = originalStepPosition;
266 } else if (maxEnergyCell == thiscell) {
267
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
281 G4ThreeVector
diff(thiscell->
x() - maxEnergyCell->
x(),
282 thiscell->
y() - maxEnergyCell->
y(),
283 thiscell->
z() - maxEnergyCell->
z());
284 stepPosition = originalStepPosition +
diff;
286 const CaloDetDescrElement* bestcell =
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
314
315
316 double time = larhit.time;
317 double energy = larhit.energy / CLHEP::MeV;
318
321 G4cout << this->GetName()
322 << " DEBUG update_map: bad identifier: " << id.getString()
323 << " skipping this hit." << G4endl;
324 }
325 continue;
326 }
328 stepPosition, id, energy, time, true,
329 numberOfProcessedHits);
330 }
331 }
332 if (madeSubSteps) {
333
334
335 while (!
steps.empty()) {
338 }
339 }
341 G4cout << "LArFCS_StepInfoSD::ProcessHits(): Etotal substeps=" << et_all
342 << G4endl << G4endl << G4endl;
343 }
344 }
345
347}
Scalar mag() const
mag method
float et(const xAOD::jFexSRJetRoI *j)
void diff(const Jet &rJet1, const Jet &rJet2, std::map< std::string, double > varDiff)
Difference between jets - Non-Class function required by trigger.
void update_map(const CLHEP::Hep3Vector &l_vec, const Identifier &l_identifier, double l_energy, double l_time, bool l_valid, int l_detector)
void getCaloDDManager()
Keep a map instead of trying to keep the full vector.
std::string getString() const
Provide a string form of the identifier - hexadecimal.
Identifier ConvertID(const LArG4Identifier &a_ident) const
Helper function for making "real" identifiers from LArG4Identifiers.
time(flags, cells_name, *args, **kw)
G4ThreeVector preStepPosition(const G4Step *theStep)
TODO.
G4ThreeVector postStepPosition(const G4Step *theStep)
TODO.