Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Implement modular biasing options #92

Draft
wants to merge 29 commits into
base: master
Choose a base branch
from
Draft
Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
29 commits
Select commit Hold shift + click to select a range
085a7cb
fail when process is not understood
lobis Dec 15, 2022
0201f6b
minor cleanup
lobis Dec 15, 2022
20efca0
add prototype for biasing example
lobis Dec 15, 2022
44f9015
Merge branch 'master' into lobis-gamma-biasing
lobis Feb 22, 2023
aef3b20
biasing templated based on example 4 (Geant4) compiles
lobis Feb 22, 2023
4d88650
working splitting in physics list
lobis Feb 22, 2023
df8707b
split only when pointing towards (0,0,0)
lobis Feb 22, 2023
0452254
attempt to integrate with TRestGeant4Metadata
lobis Feb 22, 2023
520aed1
[pre-commit.ci] auto fixes from pre-commit.com hooks
pre-commit-ci[bot] Feb 22, 2023
a5482df
attempt to fix pipeline failure
lobis Feb 22, 2023
2c66056
support for old biasing
lobis Feb 23, 2023
5165af1
update constructor order
lobis Feb 27, 2023
e61475a
update volume names
lobis Feb 27, 2023
5558ccc
[pre-commit.ci] auto fixes from pre-commit.com hooks
pre-commit-ci[bot] Feb 27, 2023
b6ac1aa
reduce number of backwards tracks
lobis Feb 27, 2023
961f86e
Merge remote-tracking branch 'origin/lobis-gamma-biasing' into lobis-…
lobis Feb 27, 2023
fb596cb
fix bad volume name
lobis Feb 27, 2023
fd1e9f8
fix bad volume name
lobis Feb 27, 2023
84a8605
[pre-commit.ci] auto fixes from pre-commit.com hooks
pre-commit-ci[bot] Feb 27, 2023
7a63d7c
grouped biasing in a single header/src file
lobis Feb 27, 2023
00e56a1
fix merge conflicts
lobis Feb 27, 2023
1a9cc40
[pre-commit.ci] auto fixes from pre-commit.com hooks
pre-commit-ci[bot] Feb 27, 2023
ae79122
add generic "cosmic" distribution (power -2.7)
lobis Feb 27, 2023
ab45454
Merge remote-tracking branch 'origin/lobis-gamma-biasing' into lobis-…
lobis Feb 27, 2023
f687f58
add option to save kill volume energy
lobis Feb 27, 2023
3b897b1
[pre-commit.ci] auto fixes from pre-commit.com hooks
pre-commit-ci[bot] Feb 27, 2023
7245b17
working biasing example
lobis Feb 28, 2023
78bc24b
Merge remote-tracking branch 'origin/lobis-gamma-biasing' into lobis-…
lobis Feb 28, 2023
9acbf0f
[pre-commit.ci] auto fixes from pre-commit.com hooks
pre-commit-ci[bot] Feb 28, 2023
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Prev Previous commit
Next Next commit
[pre-commit.ci] auto fixes from pre-commit.com hooks
for more information, see https://pre-commit.ci
  • Loading branch information
pre-commit-ci[bot] committed Feb 27, 2023
commit 3b897b11d7cb27e8c62e11b9f13b158380e13572
2 changes: 1 addition & 1 deletion examples/15.GammaBiasing/geometry/setup.gdml
Original file line number Diff line number Diff line change
Expand Up @@ -62,7 +62,7 @@
</volume>

</structure>

<setup name="Default" version="1.0">
<world ref="World"/>
</setup>
Expand Down
68 changes: 34 additions & 34 deletions src/DataModel.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -10,21 +10,21 @@

using namespace std;

TRestGeant4Event::TRestGeant4Event(const G4Event *event) : TRestGeant4Event() {
TRestGeant4Event::TRestGeant4Event(const G4Event* event) : TRestGeant4Event() {
SetID(event->GetEventID());
SetOK(true);
time_t system_time = time(nullptr);

SetTime((Double_t) system_time);
SetTime((Double_t)system_time);

auto primaryVertex = event->GetPrimaryVertex();
const auto &position = primaryVertex->GetPosition();
const auto& position = primaryVertex->GetPosition();
fPrimaryPosition = {position.x() / CLHEP::mm, position.y() / CLHEP::mm, position.z() / CLHEP::mm};
for (int i = 0; i < primaryVertex->GetNumberOfParticle(); i++) {
const auto &primaryParticle = primaryVertex->GetPrimary(i);
const auto& primaryParticle = primaryVertex->GetPrimary(i);
fPrimaryParticleNames.emplace_back(primaryParticle->GetParticleDefinition()->GetParticleName());
fPrimaryEnergies.emplace_back(primaryParticle->GetKineticEnergy() / CLHEP::keV);
const auto &momentum = primaryParticle->GetMomentumDirection();
const auto& momentum = primaryParticle->GetMomentumDirection();
fPrimaryDirections.emplace_back(momentum.x(), momentum.y(), momentum.z());
}

Expand All @@ -49,7 +49,7 @@ TRestGeant4Event::TRestGeant4Event(const G4Event *event) : TRestGeant4Event() {
*/
}

bool TRestGeant4Event::InsertTrack(const G4Track *track) {
bool TRestGeant4Event::InsertTrack(const G4Track* track) {
if (fInitialStep.GetNumberOfHits() != 1) {
cout << "fInitialStep does not have exactly one step! Problem with stepping verbose" << endl;
exit(1);
Expand All @@ -60,33 +60,33 @@ bool TRestGeant4Event::InsertTrack(const G4Track *track) {
// First track of sub-event (primary)
fSubEventPrimaryParticleName = track->GetParticleDefinition()->GetParticleName();
fSubEventPrimaryEnergy = track->GetKineticEnergy() / CLHEP::keV;
const auto &position = track->GetPosition();
const auto& position = track->GetPosition();
fSubEventPrimaryPosition = {position.x() / CLHEP::mm, position.y() / CLHEP::mm,
position.z() / CLHEP::mm};
const auto &momentum = track->GetMomentumDirection();
const auto& momentum = track->GetMomentumDirection();
fSubEventPrimaryDirection = {momentum.x(), momentum.y(), momentum.z()};
}

fTrackIDToTrackIndex[track->GetTrackID()] = int(fTracks.size()); // before insertion

fTracks.emplace_back(track);

auto &insertedTrack = fTracks.back();
auto& insertedTrack = fTracks.back();

insertedTrack.SetHits(fInitialStep);
insertedTrack.SetEvent(this);

TRestGeant4Track *parentTrack = GetTrackByID(track->GetParentID());
TRestGeant4Track* parentTrack = GetTrackByID(track->GetParentID());
if (parentTrack) {
parentTrack->AddSecondaryTrackID(track->GetTrackID());
}

return true;
}

void TRestGeant4Event::UpdateTrack(const G4Track *track) { fTracks.back().UpdateTrack(track); }
void TRestGeant4Event::UpdateTrack(const G4Track* track) { fTracks.back().UpdateTrack(track); }

void TRestGeant4Event::InsertStep(const G4Step *step) {
void TRestGeant4Event::InsertStep(const G4Step* step) {
if (step->GetTrack()->GetCurrentStepNumber() == 0) {
// initial step (from SteppingVerbose) is generated before TrackingAction can insert the first track
fInitialStep = TRestGeant4Hits();
Expand All @@ -97,11 +97,11 @@ void TRestGeant4Event::InsertStep(const G4Step *step) {
}
}

bool OutputManager::IsValidTrack(const G4Track *) const { return true; }
bool OutputManager::IsValidTrack(const G4Track*) const { return true; }

bool OutputManager::IsValidStep(const G4Step *) const { return true; }
bool OutputManager::IsValidStep(const G4Step*) const { return true; }

TRestGeant4Track::TRestGeant4Track(const G4Track *track) : TRestGeant4Track() {
TRestGeant4Track::TRestGeant4Track(const G4Track* track) : TRestGeant4Track() {
fTrackID = track->GetTrackID();
fParentID = track->GetParentID();

Expand All @@ -124,13 +124,13 @@ TRestGeant4Track::TRestGeant4Track(const G4Track *track) : TRestGeant4Track() {

fGlobalTimestamp = track->GetGlobalTime() / CLHEP::microsecond;

const G4ThreeVector &trackOrigin = track->GetPosition();
const G4ThreeVector& trackOrigin = track->GetPosition();
fInitialPosition = {trackOrigin.x(), trackOrigin.y(), trackOrigin.z()};
}

void TRestGeant4Track::InsertStep(const G4Step *step) { fHits.InsertStep(step); }
void TRestGeant4Track::InsertStep(const G4Step* step) { fHits.InsertStep(step); }

void TRestGeant4Track::UpdateTrack(const G4Track *track) {
void TRestGeant4Track::UpdateTrack(const G4Track* track) {
if (track->GetTrackID() != fTrackID) {
G4cout << "Geant4Track::UpdateTrack - mismatch of trackID!" << endl;
exit(1);
Expand All @@ -140,19 +140,19 @@ void TRestGeant4Track::UpdateTrack(const G4Track *track) {
fTimeLength = track->GetGlobalTime() / CLHEP::microsecond - fGlobalTimestamp;
}

Int_t TRestGeant4PhysicsInfo::GetProcessIDFromGeant4Process(const G4VProcess *process) {
Int_t TRestGeant4PhysicsInfo::GetProcessIDFromGeant4Process(const G4VProcess* process) {
return process->GetProcessType() * 1000 + process->GetProcessSubType();
}

void TRestGeant4Hits::InsertStep(const G4Step *step) {
const G4Track *track = step->GetTrack();
void TRestGeant4Hits::InsertStep(const G4Step* step) {
const G4Track* track = step->GetTrack();

TRestGeant4Metadata *metadata = GetGeant4Metadata();
TRestGeant4Metadata* metadata = GetGeant4Metadata();

const auto &geometryInfo = metadata->GetGeant4GeometryInfo();
const auto& geometryInfo = metadata->GetGeant4GeometryInfo();

const auto &volumeNameGeant4 = step->GetPreStepPoint()->GetPhysicalVolume()->GetName();
const auto &volumeName = geometryInfo.GetAlternativeNameFromGeant4PhysicalName(volumeNameGeant4);
const auto& volumeNameGeant4 = step->GetPreStepPoint()->GetPhysicalVolume()->GetName();
const auto& volumeName = geometryInfo.GetAlternativeNameFromGeant4PhysicalName(volumeNameGeant4);

if (!metadata->IsActiveVolume(volumeName) && step->GetTrack()->GetCurrentStepNumber() != 0) {
// we always store the first step
Expand All @@ -161,9 +161,9 @@ void TRestGeant4Hits::InsertStep(const G4Step *step) {

const bool kill = metadata->IsKillVolume(volumeName);

const auto &particle = step->GetTrack()->GetDefinition();
const auto &particleID = particle->GetPDGEncoding();
const auto &particleName = particle->GetParticleName();
const auto& particle = step->GetTrack()->GetDefinition();
const auto& particleID = particle->GetPDGEncoding();
const auto& particleName = particle->GetParticleName();

auto energy = step->GetTotalEnergyDeposit() / CLHEP::keV;

Expand Down Expand Up @@ -204,17 +204,17 @@ void TRestGeant4Hits::InsertStep(const G4Step *step) {
metadata->fGeant4PhysicsInfo.InsertProcessName(processID, processName, processTypeName);

auto sensitiveVolumeName =
geometryInfo.GetAlternativeNameFromGeant4PhysicalName(metadata->GetSensitiveVolume());
geometryInfo.GetAlternativeNameFromGeant4PhysicalName(metadata->GetSensitiveVolume());

G4Track *aTrack = step->GetTrack();
G4Track* aTrack = step->GetTrack();

Double_t x = aTrack->GetPosition().x() / CLHEP::mm;
Double_t y = aTrack->GetPosition().y() / CLHEP::mm;
Double_t z = aTrack->GetPosition().z() / CLHEP::mm;

const TVector3 hitPosition(x, y, z);
const Double_t hitGlobalTime = track->GetGlobalTime() / CLHEP::microsecond;
const G4ThreeVector &momentum = track->GetMomentumDirection();
const G4ThreeVector& momentum = track->GetMomentumDirection();

AddHit(hitPosition, energy, hitGlobalTime); // this increases fNHits

Expand All @@ -228,9 +228,9 @@ void TRestGeant4Hits::InsertStep(const G4Step *step) {
}

void OutputManager::RemoveUnwantedTracks() {
const auto &metadata = fSimulationManager->GetRestMetadata();
const auto& metadata = fSimulationManager->GetRestMetadata();
set<int> trackIDsToKeep; // We populate this container with the tracks we want to keep
for (const auto &track: fEvent->fTracks) {
for (const auto& track : fEvent->fTracks) {
// If one children track is kept, we keep all the parents
if (trackIDsToKeep.count(track.GetTrackID()) > 0) {
continue;
Expand All @@ -256,7 +256,7 @@ void OutputManager::RemoveUnwantedTracks() {
// const size_t numberOfTracksBefore = fEvent->fTracks.size();

vector<TRestGeant4Track> tracksAfterRemoval;
for (const auto &track: fEvent->fTracks) {
for (const auto& track : fEvent->fTracks) {
// we do this to preserve original order
if (trackIDsToKeep.count(track.GetTrackID()) > 0) {
tracksAfterRemoval.push_back(*(fEvent->GetTrackByID(track.GetTrackID())));
Expand Down