BDSIM
BDSIM is a Geant4 extension toolkit for simulation of particle transport in accelerator beamlines.
Loading...
Searching...
No Matches
BDSEventAction.cc
1/*
2Beam Delivery Simulation (BDSIM) Copyright (C) Royal Holloway,
3University of London 2001 - 2023.
4
5This file is part of BDSIM.
6
7BDSIM is free software: you can redistribute it and/or modify
8it under the terms of the GNU General Public License as published
9by the Free Software Foundation version 3 of the License.
10
11BDSIM is distributed in the hope that it will be useful, but
12WITHOUT ANY WARRANTY; without even the implied warranty of
13MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14GNU General Public License for more details.
15
16You should have received a copy of the GNU General Public License
17along with BDSIM. If not, see <http://www.gnu.org/licenses/>.
18*/
19#include "BDSAuxiliaryNavigator.hh"
20#include "BDSDebug.hh"
21#include "BDSEventAction.hh"
22#include "BDSEventInfo.hh"
23#include "BDSGlobalConstants.hh"
24#include "BDSHitEnergyDeposition.hh"
25#include "BDSHitEnergyDepositionExtra.hh"
26#include "BDSHitEnergyDepositionGlobal.hh"
27#include "BDSHitSampler.hh"
28#include "BDSHitThinThing.hh"
29#include "BDSOutput.hh"
30#include "BDSModulator.hh"
31#include "BDSNavigatorPlacements.hh"
32#include "BDSSamplerRegistry.hh"
33#include "BDSSamplerPlacementRecord.hh"
34#include "BDSSDApertureImpacts.hh"
35#include "BDSSDCollimator.hh"
36#include "BDSSDEnergyDeposition.hh"
37#include "BDSSDEnergyDepositionGlobal.hh"
38#include "BDSSDManager.hh"
39#include "BDSSDSampler.hh"
40#include "BDSSDSamplerCylinder.hh"
41#include "BDSSDSamplerSphere.hh"
42#include "BDSSDTerminator.hh"
43#include "BDSSDThinThing.hh"
44#include "BDSSDVolumeExit.hh"
45#include "BDSStackingAction.hh"
46#include "BDSTrajectoriesToStore.hh"
47#include "BDSTrajectory.hh"
48#include "BDSTrajectoryFilter.hh"
49#include "BDSTrajectoryPointHit.hh"
50#include "BDSTrajectoryPrimary.hh"
51#include "BDSUtilities.hh"
52#include "BDSWrapperMuonSplitting.hh"
53
54#include "globals.hh" // geant4 types / globals
55#include "G4Event.hh"
56#include "G4EventManager.hh"
57#include "G4HCofThisEvent.hh"
58#include "G4PrimaryVertex.hh"
59#include "G4PrimaryParticle.hh"
60#include "G4PropagatorInField.hh"
61#include "G4Run.hh"
62#include "G4SDManager.hh"
63#include "G4StackManager.hh"
64#include "G4THitsMap.hh"
65#include "G4TrajectoryContainer.hh"
66#include "G4TrajectoryPoint.hh"
67#include "G4TransportationManager.hh"
68
69#include <algorithm>
70#include <bitset>
71#include <chrono>
72#include <ctime>
73#include <map>
74#include <string>
75#include <vector>
76
77using namespace std::chrono;
78
79G4bool FireLaserCompton; // bool to ensure that Laserwire can only occur once in an event
80
81BDSEventAction::BDSEventAction(BDSOutput* outputIn):
82 output(outputIn),
83 samplerCollID_plane(-1),
84 samplerCollID_cylin(-1),
85 samplerCollID_sphere(-1),
86 eCounterID(-1),
87 eCounterFullID(-1),
88 eCounterVacuumID(-1),
89 eCounterTunnelID(-1),
90 eCounterWorldID(-1),
91 eCounterWorldContentsID(-1),
92 worldExitCollID(-1),
93 collimatorCollID(-1),
94 apertureCollID(-1),
95 thinThingCollID(-1),
96 startTime(0),
97 stopTime(0),
98 starts(0),
99 stops(0),
100 cpuStartTime(std::clock_t()),
101 primaryAbsorbedInCollimator(false),
102 currentEventIndex(0),
103 eventInfo(nullptr),
104 nTracks(0)
105{
107 verboseEventBDSIM = globals->VerboseEventBDSIM();
108 verboseEventStart = globals->VerboseEventStart();
109 verboseEventStop = BDS::VerboseEventStop(verboseEventStart, globals->VerboseEventContinueFor());
110 storeTrajectory = globals->StoreTrajectory();
111 storeTrajectoryAll = globals->StoreTrajectoryAll();
112 storeTrajectorySecondary = globals->StoreTrajectorySecondaryParticles();
113 trajectoryFilterLogicAND = globals->TrajectoryFilterLogicAND();
114 trajectoryEnergyThreshold = globals->StoreTrajectoryEnergyThreshold();
115 trajectoryCutZ = globals->TrajCutGTZ();
116 trajectoryCutR = globals->TrajCutLTR();
117 trajConnect = globals->TrajConnect();
118 trajParticleNameToStore = globals->StoreTrajectoryParticle();
119 trajParticleIDToStore = globals->StoreTrajectoryParticleID();
120 trajDepth = globals->StoreTrajectoryDepth();
121 trajSRangeToStore = globals->StoreTrajectoryELossSRange();
122 trajFiltersSet = globals->TrajectoryFiltersSet();
123 printModulo = globals->PrintModuloEvents();
124
125 // particleID to store in integer vector
126 std::stringstream iss(trajParticleIDToStore);
127 G4int i;
128 while (iss >> i)
129 {trajParticleIDIntToStore.push_back(i);}
130}
131
132BDSEventAction::~BDSEventAction()
133{;}
134
135void BDSEventAction::BeginOfEventAction(const G4Event* evt)
136{
137#ifdef BDSDEBUG
138 G4cout << __METHOD_NAME__ << "processing begin of event action" << G4endl;
139#endif
141 nTracks = 0;
143 BDSStackingAction::energyKilled = 0;
144 primaryAbsorbedInCollimator = false; // reset flag
145 currentEventIndex = evt->GetEventID();
146
147 BDSModulator::SetEventIndex(currentEventIndex);
148
149 // reset navigators to ensure no mis-navigating and that events are truly independent
150 BDSAuxiliaryNavigator::ResetNavigatorStates();
151 BDSNavigatorPlacements::ResetNavigatorStates();
152 // make further attempts to clear Geant4's tracking history between events to make them
153 // truly independent.
154 G4Navigator* trackingNavigator = G4TransportationManager::GetTransportationManager()->GetNavigatorForTracking();
155 trackingNavigator->ResetStackAndState();
156 G4TransportationManager* tm = G4TransportationManager::GetTransportationManager();
157 int i = 0;
158 for (auto it = tm->GetActiveNavigatorsIterator(); i < (int)tm->GetNoActiveNavigators(); it++)
159 {(*it)->ResetStackAndState(); i++;}
160 tm->GetPropagatorInField()->ClearPropagatorState(); // <- this one really makes a difference
161
162 // Unfortunately the interfaces to G4TransportationManager aren't great which makes this a bit
163 // pedantic. Also, the GetNavigator creates an exception if it doesn't find what it's looking
164 // for rather than just return a nullptr
165 G4bool samplerWorldExists = false;
166 std::vector<G4VPhysicalVolume*>::iterator worldIterator = tm->GetWorldsIterator();
167 for (G4int iw = 0; iw < (G4int)tm->GetNoWorlds(); iw++)
168 {
169 samplerWorldExists = samplerWorldExists || (*worldIterator)->GetName() == "SamplerWorld_main";
170 worldIterator++;
171 }
172 if (samplerWorldExists)
173 {
174 auto swtracker = tm->GetNavigator("SamplerWorld_main");
175 if (swtracker)
176 {swtracker->ResetStackAndState();}
177 }
178 fpEventManager->GetStackManager()->clear();
179
180 // update reference to event info
181 eventInfo = static_cast<BDSEventInfo*>(evt->GetUserInformation());
182
183 // number feedback
184 G4int currentEventID = evt->GetEventID();
185 BDSSDTerminator::eventNumber = currentEventID; // update static member of terminator
186 eventInfo->SetIndex(currentEventID);
187 if (currentEventID % printModulo == 0)
188 {G4cout << "---> Begin of event: " << currentEventID << G4endl;}
189 if (verboseEventBDSIM) // always print this out
190 {G4cout << __METHOD_NAME__ << "event #" << currentEventID << G4endl;}
191
192 // cache hit collection IDs for quicker access
193 if (samplerCollID_plane < 0)
194 { // if one is -1 then all need initialised.
195 G4SDManager* g4SDMan = G4SDManager::GetSDMpointer();
196 BDSSDManager* bdsSDMan = BDSSDManager::Instance();
197 samplerCollID_plane = g4SDMan->GetCollectionID(bdsSDMan->SamplerPlane()->GetName());
198 samplerCollID_cylin = g4SDMan->GetCollectionID(bdsSDMan->SamplerCylinder()->GetName());
199 samplerCollID_sphere = g4SDMan->GetCollectionID(bdsSDMan->SamplerSphere()->GetName());
200 eCounterID = g4SDMan->GetCollectionID(bdsSDMan->EnergyDeposition()->GetName());
201 eCounterFullID = g4SDMan->GetCollectionID(bdsSDMan->EnergyDepositionFull()->GetName());
202 eCounterVacuumID = g4SDMan->GetCollectionID(bdsSDMan->EnergyDepositionVacuum()->GetName());
203 eCounterTunnelID = g4SDMan->GetCollectionID(bdsSDMan->EnergyDepositionTunnel()->GetName());
204 eCounterWorldID = g4SDMan->GetCollectionID(bdsSDMan->EnergyDepositionWorld()->GetName());
205 eCounterWorldContentsID = g4SDMan->GetCollectionID(bdsSDMan->EnergyDepositionWorldContents()->GetName());
206 worldExitCollID = g4SDMan->GetCollectionID(bdsSDMan->WorldExit()->GetName());
207 collimatorCollID = g4SDMan->GetCollectionID(bdsSDMan->Collimator()->GetName());
208 apertureCollID = g4SDMan->GetCollectionID(bdsSDMan->ApertureImpacts()->GetName());
209 thinThingCollID = g4SDMan->GetCollectionID(bdsSDMan->ThinThing()->GetName());
210 const std::vector<G4String>& scorerNames = bdsSDMan->PrimitiveScorerNamesComplete();
211 for (const auto& name : scorerNames)
212 {scorerCollectionIDs[name] = g4SDMan->GetCollectionID(name);}
213 const std::vector<G4String>& extraSamplerWithFilterNames = bdsSDMan->ExtraSamplerWithFilterNamesComplete();
214 for (const auto& name : extraSamplerWithFilterNames)
215 {extraSamplerCollectionIDs[name] = g4SDMan->GetCollectionID(name);}
216 const std::vector<G4String>& extraSamplerCylinderWithFilterNames = bdsSDMan->ExtraSamplerCylinderWithFilterNamesComplete();
217 for (const auto& name : extraSamplerCylinderWithFilterNames)
218 {extraSamplerCylinderCollectionIDs[name] = g4SDMan->GetCollectionID(name);}
219 const std::vector<G4String>& extraSamplerSphereWithFilterNames = bdsSDMan->ExtraSamplerSphereWithFilterNamesComplete();
220 for (const auto& name : extraSamplerSphereWithFilterNames)
221 {extraSamplerCollectionIDs[name] = g4SDMan->GetCollectionID(name);}
222 }
223 FireLaserCompton=true;
224
225 cpuStartTime = std::clock();
226 // get the current time - last thing before we hand off to geant4
227 startTime = time(nullptr);
229 eventInfo->SetStopTime(startTime); // initialise to duration of 0
230
231 milliseconds ms = duration_cast<milliseconds>(system_clock::now().time_since_epoch());
232 starts = (G4double)ms.count()/1000.0;
233}
234
235void BDSEventAction::EndOfEventAction(const G4Event* evt)
236{
237 //G4cout << "BDSWrapperMuonSplitting::nCallsThisEvent> " << BDSWrapperMuonSplitting::nCallsThisEvent << G4endl;
238 auto flagsCache(G4cout.flags());
239 // Get event number information
240 G4int event_number = evt->GetEventID();
241 G4bool verboseThisEvent = verboseEventBDSIM && BDS::VerboseThisEvent(event_number, verboseEventStart, verboseEventStop);
242#ifdef BDSDEBUG
243 verboseThisEvent = true; // force on for debug mode
244#endif
245 const G4int nChar = 50; // for print out
246 if (verboseThisEvent)
247 {G4cout << __METHOD_NAME__ << "processing end of event"<<G4endl;}
248 eventInfo->SetIndex(event_number);
249
250 // Record if event was aborted - ie whether it's usable for analyses.
251 eventInfo->SetAborted(evt->IsAborted());
253
254 // Calculate the elapsed CPU time for the event.
255 auto cpuEndTime = std::clock();
256 G4float durationCPU = static_cast<G4float>(cpuEndTime - cpuStartTime) / CLOCKS_PER_SEC;
257 eventInfo->SetDurationCPU(durationCPU);
258
259 // Get the current wall time
260 stopTime = time(nullptr);
262
263 // Timing information (wall)
264 milliseconds ms = duration_cast<milliseconds>(system_clock::now().time_since_epoch());
265 stops = (G4double)ms.count()/1000.0;
267
268 G4double memoryUsedMb = BDS::GetMemoryUsage();
269 eventInfo->SetMemoryUsage(memoryUsedMb);
270
271 // cache if primary was absorbed in a collimator
273
274 // feedback
275 if (verboseThisEvent)
276 {eventInfo->Print();}
277
278 // Get the hits collection of this event - all hits from different SDs.
279 G4HCofThisEvent* HCE = evt->GetHCofThisEvent();
280
281 // samplers
282 typedef BDSHitsCollectionSampler shc;
283 std::vector<shc*> allSamplerHits;
284 shc* SampHC = HCE ? dynamic_cast<shc*>(HCE->GetHC(samplerCollID_plane)) : nullptr;
285 if (SampHC)
286 {allSamplerHits.push_back(SampHC);}
287 if (HCE)
288 {
289 for (const auto& nameIndex : extraSamplerCollectionIDs)
290 {allSamplerHits.push_back(dynamic_cast<shc*>(HCE->GetHC(nameIndex.second)));}
291 }
293 std::vector<shcc*> allSamplerCylinderHits;
294 shcc* hitsCylinder = HCE ? dynamic_cast<shcc*>(HCE->GetHC(samplerCollID_cylin)) : nullptr;
295 if (hitsCylinder)
296 {allSamplerCylinderHits.push_back(hitsCylinder);}
297 if (HCE)
298 {
299 for (const auto& nameIndex : extraSamplerCylinderCollectionIDs)
300 {allSamplerCylinderHits.push_back(dynamic_cast<shcc*>(HCE->GetHC(nameIndex.second)));}
301 }
303 std::vector<shcs*> allSamplerSphereHits;
304 shcs* hitsSphere = HCE ? dynamic_cast<shcs*>(HCE->GetHC(samplerCollID_sphere)) : nullptr;
305 if (hitsSphere)
306 {allSamplerSphereHits.push_back(hitsSphere);}
307 if (HCE)
308 {
309 for (const auto& nameIndex : extraSamplerSphereCollectionIDs)
310 {allSamplerSphereHits.push_back(dynamic_cast<shcs*>(HCE->GetHC(nameIndex.second)));}
311 }
312
313 // energy deposition collections - eloss, tunnel hits
315 echc* eCounterHits = HCE ? dynamic_cast<echc*>(HCE->GetHC(eCounterID)) : nullptr;
316 echc* eCounterFullHits = HCE ? dynamic_cast<echc*>(HCE->GetHC(eCounterFullID)) : nullptr;
317 echc* eCounterVacuumHits = HCE ? dynamic_cast<echc*>(HCE->GetHC(eCounterVacuumID)) : nullptr;
318 echc* eCounterTunnelHits = HCE ? dynamic_cast<echc*>(HCE->GetHC(eCounterTunnelID)) : nullptr;
319
320 // world exit hits
322 ecghc* eCounterWorldHits = HCE ? dynamic_cast<ecghc*>(HCE->GetHC(eCounterWorldID)) : nullptr;
323 ecghc* eCounterWorldContentsHits = HCE ? dynamic_cast<ecghc*>(HCE->GetHC(eCounterWorldContentsID)) : nullptr;
324 ecghc* worldExitHits = HCE ? dynamic_cast<ecghc*>(HCE->GetHC(worldExitCollID)) : nullptr;
325
326 // aperture hits
328 aihc* apertureImpactHits = HCE ? dynamic_cast<aihc*>(HCE->GetHC(apertureCollID)) : nullptr;
329
330 // thin thing hits
331 typedef BDSHitsCollectionThinThing tthc;
332 tthc* thinThingHits = HCE ? dynamic_cast<tthc*>(HCE->GetHC(thinThingCollID)) : nullptr;
333
334 std::map<G4String, G4THitsMap<G4double>*> scorerHits;
335 if (HCE)
336 {
337 for (const auto& nameIndex : scorerCollectionIDs)
338 {scorerHits[nameIndex.first] = dynamic_cast<G4THitsMap<G4double>*>(HCE->GetHC(nameIndex.second));}
339 }
340 // primary hit something? we infer this by seeing if there are any energy
341 // deposition hits at all - if there are, the primary must have 'hit' something.
342 // we don't check the world energy hits here because the hits could be from
343 // intended transport through air in part of the machine (a gap). similarly, there
344 // could be ionisation of the vacuum gas without a real impact so we don't check
345 // the vacuum energy deposition
346 if (eCounterHits)
347 {
348 if (verboseThisEvent)
349 {G4cout << std::left << std::setw(nChar) << "Energy deposition hits: " << eCounterHits->entries() << G4endl;}
350 if (eCounterHits->entries() > 0)
352 }
353 if (eCounterFullHits)
354 {
355 if (verboseThisEvent)
356 {G4cout << std::left << std::setw(nChar) << "Energy deposition full hits: " << eCounterFullHits->entries() << G4endl;}
357 if (eCounterFullHits->entries() > 0)
359 }
360 if (eCounterTunnelHits)
361 {
362 if (verboseThisEvent)
363 {G4cout << std::left << std::setw(nChar) << "Tunnel energy deposition hits: " << eCounterTunnelHits->entries() << G4endl;}
364 if (eCounterTunnelHits->entries() > 0)
366 }
367
368 // collimator hits if any
369 typedef BDSHitsCollectionCollimator chc;
370 chc* collimatorHits = HCE ? dynamic_cast<chc*>(HCE->GetHC(collimatorCollID)) : nullptr;
371
372 if (verboseThisEvent)
373 {
374 if (eCounterVacuumHits)
375 {G4cout << std::left << std::setw(nChar) << "Vacuum energy deposition hits: " << eCounterVacuumHits->entries() << G4endl;}
376 if (eCounterWorldHits)
377 {G4cout << std::left << std::setw(nChar) << "World energy deposition hits: " << eCounterWorldHits->entries() << G4endl;}
378 if (eCounterWorldContentsHits)
379 {G4cout << std::left << std::setw(nChar) << "World contents energy deposition hits: " << eCounterWorldContentsHits->entries() << G4endl;}
380 if (worldExitHits)
381 {G4cout << std::left << std::setw(nChar) << "World exit hits: " << worldExitHits->entries() << G4endl;}
382 if (collimatorHits)
383 {G4cout << std::left << std::setw(nChar) << "Collimator hits: " << collimatorHits->entries() << G4endl;}
384 }
385
386 // primary hits and losses
387 std::vector<const BDSTrajectoryPointHit*> primaryHits;
388 std::vector<const BDSTrajectoryPointHit*> primaryLosses;
389 G4TrajectoryContainer* trajCont = evt->GetTrajectoryContainer();
390 if (trajCont)
391 {
392 if (verboseThisEvent)
393 {G4cout << std::left << std::setw(nChar) << "Trajectories: " << trajCont->size() << G4endl;}
394 for (auto p : primaryTrajectoriesCache)
395 {primaryLosses.push_back(new BDSTrajectoryPointHit(p.second, p.second->LastPoint()));}
396 }
397 std::vector<const BDSTrajectoryPrimary*> primaryTrajectories;
398 for (auto kv : primaryTrajectoriesCache)
399 {primaryTrajectories.push_back(kv.second);}
400 primaryHits = BDSHitThinThing::ResolvePossibleEarlierThinHits(primaryTrajectories, thinThingHits);
401
402 BDSTrajectoriesToStore* interestingTrajectories = IdentifyTrajectoriesForStorage(evt,
403 verboseEventBDSIM || verboseThisEvent,
404 eCounterHits,
405 eCounterFullHits,
406 allSamplerHits,
407 nChar);
408
410 evt->GetPrimaryVertex(),
411 allSamplerHits,
412 allSamplerCylinderHits,
413 allSamplerSphereHits,
414 nullptr,
415 eCounterHits,
416 eCounterFullHits,
417 eCounterVacuumHits,
418 eCounterTunnelHits,
419 eCounterWorldHits,
420 eCounterWorldContentsHits,
421 worldExitHits,
422 primaryHits,
423 primaryLosses,
424 interestingTrajectories,
425 collimatorHits,
426 apertureImpactHits,
427 scorerHits,
428 BDSGlobalConstants::Instance()->TurnsTaken());
429
430 // if events per ntuples not set (default 0) - only write out at end
431 G4int evntsPerNtuple = BDSGlobalConstants::Instance()->NumberOfEventsPerNtuple();
432 if (evntsPerNtuple>0 && (event_number+1)%evntsPerNtuple == 0)
433 {
434 // note the timing information will be wrong here as the run hasn't finished but
435 // the file is bridged. There's no good way around this just now as this class
436 // can't access the timing information stored in BDSRunAction
438 }
439
440 if (verboseThisEvent)
441 {
442 G4cout << __METHOD_NAME__ << "end of event action done" << G4endl;
443 G4cout << "Energy deposition pool size: " << BDSAllocatorEnergyDeposition.GetAllocatedSize() << G4endl;
444 G4cout << "Energy deposition extra pool size: " << BDSAllocatorEnergyDepositionExtra.GetAllocatedSize() << G4endl;
445 G4cout << "Collimator hits pool size: " << BDSAllocatorCollimator.GetAllocatedSize() << G4endl;
446 G4cout << "Trajectory pool size: " << bdsTrajectoryAllocator.GetAllocatedSize() << G4endl;
447 G4cout << "Trajectory point pool size bdsim: " << bdsTrajectoryPointAllocator.GetAllocatedSize() << G4endl;
448#if G4VERSION_NUMBER > 1049
449 G4cout << "Trajectory point pool size: " << aTrajectoryPointAllocator()->GetAllocatedSize() << G4endl;
450#else
451 G4cout << "Trajectory point pool size: " << aTrajectoryPointAllocator->GetAllocatedSize() << G4endl;
452#endif
453 G4cout << "Trajectory point primary pool size: " << bdsTrajectoryPrimaryAllocator.GetAllocatedSize() << G4endl;
454 }
455
456 delete interestingTrajectories;
457 for (auto& p : primaryLosses)
458 {delete p;}
459 for (auto& p : primaryHits)
460 {delete p;}
461
462 G4cout.flags(flagsCache);
463}
464
466 G4bool verbose,
468 BDSHitsCollectionEnergyDeposition* eCounterFullHits,
469 const std::vector<BDSHitsCollectionSampler*>& allSamplerHits,
470 G4int nChar) const
471{
472 auto flagsCache(G4cout.flags());
473 G4TrajectoryContainer* trajCont = evt->GetTrajectoryContainer();
474
475 // Save interesting trajectories
476 std::map<BDSTrajectory*, bool> interestingTraj;
477 std::map<BDSTrajectory*, std::bitset<BDS::NTrajectoryFilters> > trajectoryFilters;
478
479 if (storeTrajectory && trajCont)
480 {
481 TrajectoryVector* trajVec = trajCont->GetVector();
482
483 // build trackID map, depth map
484 std::map<int, BDSTrajectory*> trackIDMap;
485 std::map<BDSTrajectory*, int> depthMap;
486 for (auto iT1 : *trajVec)
487 {
488 BDSTrajectory* traj = static_cast<BDSTrajectory*>(iT1);
489
490 // fill track ID map
491 trackIDMap[traj->GetTrackID()] = traj;
492
493 // fill depth map
494 G4int depth = traj->GetParentID() == 0 ? 0 : depthMap.at(trackIDMap.at(traj->GetParentID())) + 1;
495 traj->SetDepth(depth);
496 if (traj->GetParentID() == 0)
497 {depthMap[traj] = 0;}
498 else
499 {depthMap[traj] = depthMap.at(trackIDMap.at(traj->GetParentID())) + 1;}
500 }
501
502 // fill parent pointer - this can only be done once the map in the previous loop has been made
503 for (auto iT1 : *trajVec)
504 {
505 BDSTrajectory* traj = static_cast<BDSTrajectory*>(iT1);
506 // the parent ID may be 0 and therefore it may not be in the map but the [] operator
507 // will default-construct it giving therefore a nullptr
508 traj->SetParent(trackIDMap[iT1->GetParentID()]);
509 }
510
511 // loop over trajectories and determine if it should be stored
512 // keep tallies of how many will and won't be stored
513 G4int nYes = 0;
514 G4int nNo = 0;
515 for (auto iT1 : *trajVec)
516 {
517 std::bitset<BDS::NTrajectoryFilters> filters;
518
519 BDSTrajectory* traj = static_cast<BDSTrajectory*>(iT1);
520 G4int parentID = traj->GetParentID();
521
522 // always store primaries
523 if (parentID == 0)
524 {filters[BDSTrajectoryFilter::primary] = true;}
525 else
526 {
527 if (storeTrajectorySecondary)
528 {filters[BDSTrajectoryFilter::secondary] = true;}
529 }
530
531 // check on energy (if energy threshold is not negative)
532 if (trajectoryEnergyThreshold >= 0 &&
533 traj->GetInitialKineticEnergy() > trajectoryEnergyThreshold)
534 {filters[BDSTrajectoryFilter::energyThreshold] = true;}
535
536 // check on particle if not empty string
537 if (!trajParticleNameToStore.empty() || !trajParticleIDToStore.empty())
538 {
539 G4String particleName = traj->GetParticleName();
540 G4int particleID = traj->GetPDGEncoding();
541 G4String particleIDStr = G4String(std::to_string(particleID));
542 std::size_t found1 = trajParticleNameToStore.find(particleName);
543 bool found2 = (std::find(trajParticleIDIntToStore.begin(), trajParticleIDIntToStore.end(), particleID)
544 != trajParticleIDIntToStore.end());
545 if ((found1 != std::string::npos) || found2)
546 {filters[BDSTrajectoryFilter::particle] = true;}
547 }
548
549 // check on trajectory tree depth (trajDepth = 0 means only primaries)
550 if (depthMap[traj] <= trajDepth || storeTrajectoryAll) // all means to infinite trajDepth really
551 {filters[BDSTrajectoryFilter::depth] = true;}
552
553 // check on coordinates (and TODO momentum)
554 // clear out trajectories that don't reach point TrajCutGTZ or greater than TrajCutLTR
555 BDSTrajectoryPoint* trajEndPoint = static_cast<BDSTrajectoryPoint*>(traj->GetPoint(traj->GetPointEntries() - 1));
556
557 // end point greater than some Z
558 if (trajEndPoint->GetPosition().z() > trajectoryCutZ)
559 {filters[BDSTrajectoryFilter::minimumZ] = true;}
560
561 // less than maximum R
562 if (trajEndPoint->PostPosR() < trajectoryCutR)
563 {filters[BDSTrajectoryFilter::maximumR] = true;}
564
565 filters.any() ? nYes++ : nNo++;
566 interestingTraj.insert(std::pair<BDSTrajectory*, bool>(traj, filters.any()));
567 trajectoryFilters.insert(std::pair<BDSTrajectory*, std::bitset<BDS::NTrajectoryFilters> >(traj, filters));
568 }
569
570 // loop over energy hits to connect trajectories
571 if (!trajSRangeToStore.empty())
572 {
573 if (eCounterHits)
574 {
575 G4int nHits = (G4int)eCounterHits->entries();
577 for (G4int i = 0; i < nHits; i++)
578 {
579 hit = (*eCounterHits)[i];
580 double dS = hit->GetSHit();
581 for (const auto& v : trajSRangeToStore)
582 {
583 if ( dS >= v.first && dS <= v.second)
584 {
585 BDSTrajectory* trajToStore = trackIDMap[hit->GetTrackID()];
586 if (!interestingTraj[trajToStore])
587 {// was marked as not storing - update counters
588 nYes++;
589 nNo--;
590 }
591 interestingTraj[trajToStore] = true;
592 trajectoryFilters[trajToStore][BDSTrajectoryFilter::elossSRange] = true;
593 break;
594 }
595 }
596 }
597 }
598 if (eCounterFullHits)
599 {
600 G4int nHits = (G4int)eCounterFullHits->entries();
602 for (G4int i = 0; i < nHits; i++)
603 {
604 hit = (*eCounterFullHits)[i];
605 double dS = hit->GetSHit();
606 for (const auto& v : trajSRangeToStore)
607 {
608 if ( dS >= v.first && dS <= v.second)
609 {
610 BDSTrajectory* trajToStore = trackIDMap[hit->GetTrackID()];
611 if (!interestingTraj[trajToStore])
612 {// was marked as not storing - update counters
613 nYes++;
614 nNo--;
615 }
616 interestingTraj[trajToStore] = true;
617 trajectoryFilters[trajToStore][BDSTrajectoryFilter::elossSRange] = true;
618 break;
619 }
620 }
621 }
622 }
623 }
624
625 // loop over samplers to connect trajectories
626 for (const auto& SampHC : allSamplerHits)
627 {
628 if (!trajectorySamplerID.empty() && SampHC)
629 {
630 for (G4int i = 0; i < (G4int)SampHC->entries(); i++)
631 {
632 G4int samplerIndex = (*SampHC)[i]->samplerID;
634 if (std::find(trajectorySamplerID.begin(), trajectorySamplerID.end(), samplerIndex) !=
636 {
637 BDSTrajectory* trajToStore = trackIDMap[(*SampHC)[i]->trackID];
638 if (!interestingTraj[trajToStore])
639 {// was marked as not storing - update counters
640 nYes++;
641 nNo--;
642 }
643 interestingTraj[trajToStore] = true;
644 trajectoryFilters[trajToStore][BDSTrajectoryFilter::sampler] = true;
645 }
646 }
647 }
648 }
649
650 // If we're using AND logic (default OR) with the filters, loop over and update whether
651 // we should really store the trajectory or not. Importantly, we do this before the connect
652 // trajectory step as that flags yet more trajectories (that connect each one) back to the
653 // primary
655 {
656 for (auto& trajFlag : interestingTraj)
657 {
658 if (trajFlag.second) // if we're going to store it check the logic
659 {
660 // Use bit-wise AND ('&') on the filters matched for this trajectory with the
661 // filters set. If count of 1s the same, then trajectory should be stored,
662 // therefore if not the same, it should be set to false.
663 auto filterMatch = trajectoryFilters[trajFlag.first] & trajFiltersSet;
664 if (filterMatch.count() != trajFiltersSet.count())
665 {trajFlag.second = false;}
666 }
667 }
668 }
669
670 // Connect trajectory graphs
671 if (trajConnect && trackIDMap.size() > 1)
672 {
673 for (auto i : interestingTraj)
674 {
675 if (i.second)
676 {ConnectTrajectory(interestingTraj, i.first, trajectoryFilters);}
677 }
678 }
679 // Output interesting trajectories
680 if (verbose)
681 {G4cout << std::left << std::setw(nChar) << "Trajectories for storage: " << nYes << " out of " << nYes + nNo << G4endl;}
682 }
683 G4cout.flags(flagsCache);
684 return new BDSTrajectoriesToStore(interestingTraj, trajectoryFilters);
685}
686
687void BDSEventAction::ConnectTrajectory(std::map<BDSTrajectory*, bool>& interestingTraj,
688 BDSTrajectory* trajectoryToConnect,
689 std::map<BDSTrajectory*, std::bitset<BDS::NTrajectoryFilters> >& trajectoryFilters) const
690{
691 BDSTrajectory* parentTrajectory = trajectoryToConnect->GetParent();
692 if (parentTrajectory)
693 {
694 interestingTraj[parentTrajectory] = true;
695 trajectoryFilters[parentTrajectory][BDSTrajectoryFilter::connect] = true;
696 ConnectTrajectory(interestingTraj, parentTrajectory, trajectoryFilters);
697 }
698 else
699 {return;}
700}
701
703{
704 G4int trackID = trajectoryIn->GetTrackID();
705 if (primaryTrajectoriesCache.find(trackID) == primaryTrajectoriesCache.end())
706 {primaryTrajectoriesCache[trackID] = trajectoryIn;}
707}
G4String trajParticleIDToStore
Cache of variable from global constants.
void ConnectTrajectory(std::map< BDSTrajectory *, bool > &interestingTraj, BDSTrajectory *trajectoryToConnect, std::map< BDSTrajectory *, std::bitset< BDS::NTrajectoryFilters > > &trajectoryFilters) const
G4int eCounterVacuumID
Collection ID for the vacuum energy deposition hits.
G4bool primaryAbsorbedInCollimator
Whether primary stopped in a collimator.
G4double trajectoryCutZ
Cache of variable from global constants.
G4int collimatorCollID
Collection ID for the collimator hits.
G4int eCounterFullID
Collection ID for general energy deposition full hits.
G4int worldExitCollID
Collection ID for the world exit hits.
G4int thinThingCollID
Collection ID for the thin thing hits.
std::vector< std::pair< double, double > > trajSRangeToStore
Cache of variable from global constants.
G4bool storeTrajectory
Cache of whether to store trajectories or not.
G4bool trajectoryFilterLogicAND
Cache of variable from global constants.
G4int eCounterWorldID
Collection ID for the world energy deposition hits.
std::bitset< BDS::NTrajectoryFilters > trajFiltersSet
Cache of variable from global constants.
std::vector< int > trajParticleIDIntToStore
Cache of variable from global constants.
G4int samplerCollID_cylin
Collection ID for cylindrical sampler hits.
std::map< G4String, G4int > extraSamplerCylinderCollectionIDs
Collection IDs for extra samplers.
G4int eCounterWorldContentsID
Collection ID for the world energy deposition hits.
BDSOutput * output
Cache of output instance. Not owned by this class.
BDSTrajectoriesToStore * IdentifyTrajectoriesForStorage(const G4Event *evt, G4bool verboseThisEvent, BDSHitsCollectionEnergyDeposition *eCounterHits, BDSHitsCollectionEnergyDeposition *eCounterFullHits, const std::vector< BDSHitsCollectionSampler * > &allSamplerHits, G4int nChar=50) const
Sift through all trajectories (if any) and mark for storage.
G4int samplerCollID_plane
Collection ID for plane sampler hits.
long long int nTracks
Accumulated number of tracks for the event.
std::map< G4String, G4int > extraSamplerSphereCollectionIDs
Collection IDs for extra samplers.
time_t startTime
Time at the start of the event.
std::map< G4String, G4int > scorerCollectionIDs
Collection IDs for all scorers.
G4bool trajConnect
Cache of variable from global constants.
BDSEventInfo * eventInfo
G4int apertureCollID
Collection ID for the aperture hits.
std::vector< int > trajectorySamplerID
Cache of variable from global constants.
std::map< G4String, G4int > extraSamplerCollectionIDs
Collection IDs for extra samplers.
G4double trajectoryCutR
Cache of variable from global constants.
G4String trajParticleNameToStore
Cache of variable from global constants.
G4double stops
Precise stop time in seconds.
std::clock_t cpuStartTime
CPU time at the start of the event.
G4int trajDepth
Cache of variable from global constants.
std::map< G4int, const BDSTrajectoryPrimary * > primaryTrajectoriesCache
G4int eCounterID
Collection ID for general energy deposition hits.
void RegisterPrimaryTrajectory(const BDSTrajectoryPrimary *trajectoryIn)
Append this trajectory to vector of primaries we keep to avoid sifting at the end of event.
G4int samplerCollID_sphere
Collection ID for spherical sampler hits.
G4bool storeTrajectoryAll
Store all trajectories irrespective of filters.
time_t stopTime
Time at the end of the event.
G4double starts
Precise start time in seconds.
G4int eCounterTunnelID
Collection ID for the tunnel energy deposition hits.
G4double trajectoryEnergyThreshold
Cache of variable from global constants.
Interface to store event information use G4 hooks.
Definition: BDSEventInfo.hh:42
void SetPrimaryAbsorbedInCollimator(G4bool absorbed)
Setters.
Definition: BDSEventInfo.hh:59
void SetNTracks(long long int nTracks)
Setters.
Definition: BDSEventInfo.hh:60
void SetStopTime(const time_t &stopTimeIn)
Setters.
Definition: BDSEventInfo.hh:51
void SetMemoryUsage(G4double memoryUsageMbIn)
Setters.
Definition: BDSEventInfo.hh:58
void SetPrimaryHitMachine(G4bool hitIn)
Setters.
Definition: BDSEventInfo.hh:57
void SetStartTime(const time_t &startTimeIn)
Setters.
Definition: BDSEventInfo.hh:50
void SetIndex(G4int indexIn)
Setters.
Definition: BDSEventInfo.hh:55
void SetAborted(G4bool abortedIn)
Setters.
Definition: BDSEventInfo.hh:56
void SetDurationCPU(G4float durationCPUIn)
Setters.
Definition: BDSEventInfo.hh:53
void SetDurationWall(G4float durationWallIn)
Setters.
Definition: BDSEventInfo.hh:52
A class that holds global options and constants.
static BDSGlobalConstants * Instance()
Access method.
Information recorded for a single piece of energy deposition.
G4int GetTrackID() const
Accessor for extra piece of information.
static void SetEventIndex(G4int eventIndexIn)
Setter for the static counter so we can trace where it's used.
Definition: BDSModulator.cc:23
Output base class that defines interface for all output types.
Definition: BDSOutput.hh:73
void CloseAndOpenNewFile()
Close a file and open a new one.
Definition: BDSOutput.cc:372
void FillEvent(const BDSEventInfo *info, const G4PrimaryVertex *vertex, const std::vector< BDSHitsCollectionSampler * > &samplerHitsPlane, const std::vector< BDSHitsCollectionSamplerCylinder * > &samplerHitsCylinder, const std::vector< BDSHitsCollectionSamplerSphere * > &samplerHitsSphere, const BDSHitsCollectionSamplerLink *samplerHitsLink, const BDSHitsCollectionEnergyDeposition *energyLoss, const BDSHitsCollectionEnergyDeposition *energyLossFull, const BDSHitsCollectionEnergyDeposition *energyLossVacuum, const BDSHitsCollectionEnergyDeposition *energyLossTunnel, const BDSHitsCollectionEnergyDepositionGlobal *energyLossWorld, const BDSHitsCollectionEnergyDepositionGlobal *energyLossWorldContents, const BDSHitsCollectionEnergyDepositionGlobal *worldExitHits, const std::vector< const BDSTrajectoryPointHit * > &primaryHits, const std::vector< const BDSTrajectoryPointHit * > &primaryLosses, const BDSTrajectoriesToStore *trajectories, const BDSHitsCollectionCollimator *collimatorHits, const BDSHitsCollectionApertureImpacts *apertureImpactHits, const std::map< G4String, G4THitsMap< G4double > * > &scorerHitsMap, const G4int turnsTaken)
Copy event information from Geant4 simulation structures to output structures.
Definition: BDSOutput.cc:297
A singleton class that holds all required sensitive detector class instances.
Definition: BDSSDManager.hh:68
BDSSDEnergyDeposition * EnergyDeposition() const
SD for general energy deposition.
BDSSDVolumeExit * WorldExit() const
SD for world exit hits.
BDSSDEnergyDeposition * EnergyDepositionTunnel() const
SD for tunnel energy counter.
const std::vector< G4String > & ExtraSamplerCylinderWithFilterNamesComplete() const
Access a vector of names of extra samplers so we can identify the hits collections.
const std::vector< G4String > & ExtraSamplerSphereWithFilterNamesComplete() const
Access a vector of names of extra samplers so we can identify the hits collections.
BDSSDEnergyDeposition * EnergyDepositionFull() const
SD for general energy deposition but always include extra half of information.
BDSSDCollimator * Collimator() const
SD for collimator impact locations.
BDSSDThinThing * ThinThing() const
BDSSDApertureImpacts * ApertureImpacts() const
SD for aperture impact hits.
BDSSDSamplerCylinder * SamplerCylinder() const
SD for samplers (cylinder type).
Definition: BDSSDManager.hh:89
const std::vector< G4String > & PrimitiveScorerNamesComplete() const
Access a vector the full primitive scorer names as registered.
BDSSDSampler * SamplerPlane() const
SD for samplers (plane type). See also SamplerPlaneWithFilter below.
Definition: BDSSDManager.hh:86
BDSSDSamplerSphere * SamplerSphere() const
SD for samplers (sphere type).
Definition: BDSSDManager.hh:92
BDSSDEnergyDepositionGlobal * EnergyDepositionWorldContents() const
SD for energy deposition in things that were already placed in the externally provided world.
BDSSDEnergyDepositionGlobal * EnergyDepositionWorld() const
SD for energy deposition in the world volume.
BDSSDEnergyDeposition * EnergyDepositionVacuum() const
SD for energy deposition in vacuum volumes.
const std::vector< G4String > & ExtraSamplerWithFilterNamesComplete() const
Access a vector of names of extra samplers so we can identify the hits collections.
static G4int eventNumber
Externally accessible counter for event number. Set in BeginOfEventAction.
Information about a registered sampler.
static BDSSamplerRegistry * Instance()
Accessor for registry.
const BDSSamplerPlacementRecord & GetInfo(G4int index) const
Accessor.
Double map of trajectories to bitset of which filters matched whether to store them.
A summary trajectory object of a loss point.
A Point in a trajectory with extra information.
G4double PostPosR() const
Return the transverse local radius in x,y.
Trajectory information for only the primary.
Trajectory information from track including last scatter etc.
virtual int GetPointEntries() const
Get number of trajectory points in this trajectory.
virtual G4VTrajectoryPoint * GetPoint(G4int i) const
Access a point - use this class's container.
void SetParent(BDSTrajectory *parentIn)
Record the parent trajectory.
static G4int nCallsThisEvent
Counter for understanding occurence.
G4double GetMemoryUsage()
Get the current memory usage.
G4bool VerboseThisEvent(G4int eventIndex, G4int eventStart, G4int eventStop)
Logic of whether this event should be verbose or not. Code here so it's not duplicated.
G4int VerboseEventStop(G4int verboseEventStart, G4int verboseEventContinueFor)