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.

Sign up for GitHub

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

Jump to bottom

Fixed fast sim sampling of ECAL hits #34485

Merged
merged 2 commits into from Jul 15, 2021
Merged

Fixed fast sim sampling of ECAL hits #34485

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
f7a9186
added fixes
civanch Jul 14, 2021
e8e72ea
code checks
civanch Jul 14, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
4 changes: 2 additions & 2 deletions SimG4CMS/Calo/src/CaloSD.cc
View file
Open in desktop
Expand Up @@ -598,7 +598,7 @@ CaloG4Hit* CaloSD::createNewHit(const G4Step* aStep, const G4Track* theTrack) {
storeHit(aHit);
TrackInformation* trkInfo = cmsTrackInformation(theTrack);

bool currentlyInsideFineVolume = isItFineCalo(aStep->GetPostStepPoint()->GetTouchable());
bool currentlyInsideFineVolume = !doFineCalo_ ? false : isItFineCalo(aStep->GetPostStepPoint()->GetTouchable());

#ifdef EDM_ML_DEBUG
edm::LogVerbatim("DoFineCalo") << "Creating new hit " << aHit->getUnitID() << " using "
Expand All @@ -617,7 +617,7 @@ CaloG4Hit* CaloSD::createNewHit(const G4Step* aStep, const G4Track* theTrack) {

// If fine calo is activated for the current volume, perform track/hit
// saving logic for fineCalo
if (doFineCalo_ && currentlyInsideFineVolume) {
if (currentlyInsideFineVolume) {
hitBookkeepingFineCalo(aStep, theTrack, aHit);
}
// 'Traditional', non-fine history bookkeeping
Expand Down
5 changes: 4 additions & 1 deletion SimG4Core/Application/interface/LowEnergyFastSimModel.h
View file
Open in desktop
Expand Up @@ -9,6 +9,7 @@
#include "GFlashHitMaker.hh"
#include "G4Region.hh"
#include "G4Types.hh"
#include "G4ThreeVector.hh"

class TrackingAction;

Expand All @@ -24,8 +25,10 @@ class LowEnergyFastSimModel : public G4VFastSimulationModel {
G4double fEmax;
const G4Envelope* fRegion;
const TrackingAction* fTrackingAction;
G4bool fCheck;
G4ThreeVector fTailPos;
GFlashHitMaker fHitMaker;
LowEnergyFastSimParam param;
LowEnergyFastSimParam fParam;
};

#endif
40 changes: 23 additions & 17 deletions SimG4Core/Application/src/LowEnergyFastSimModel.cc
View file
Open in desktop
Expand Up @@ -11,10 +11,14 @@
#include "G4Region.hh"
#include "G4PhysicalConstants.hh"

constexpr double twomass = 2 * CLHEP::electron_mass_c2;
constexpr G4double twomass = 2 * CLHEP::electron_mass_c2;

LowEnergyFastSimModel::LowEnergyFastSimModel(const G4String& name, G4Region* region, const edm::ParameterSet& parSet)
: G4VFastSimulationModel(name, region), fRegion(region), fTrackingAction(nullptr) {
: G4VFastSimulationModel(name, region),
fRegion(region),
fTrackingAction(nullptr),
fCheck(false),
fTailPos(0., 0., 0.) {
fEmax = parSet.getParameter<double>("LowEnergyGflashEcalEmax") * CLHEP::GeV;
}

Expand All @@ -24,32 +28,32 @@ G4bool LowEnergyFastSimModel::IsApplicable(const G4ParticleDefinition& particle)

G4bool LowEnergyFastSimModel::ModelTrigger(const G4FastTrack& fastTrack) {
const G4Track* track = fastTrack.GetPrimaryTrack();
if (nullptr == fTrackingAction) {
fTrackingAction = static_cast<const TrackingAction*>(G4EventManager::GetEventManager()->GetUserTrackingAction());
if (fCheck) {
if (nullptr == fTrackingAction) {
fTrackingAction = static_cast<const TrackingAction*>(G4EventManager::GetEventManager()->GetUserTrackingAction());
}
int pdgMother = std::abs(fTrackingAction->geant4Track()->GetDefinition()->GetPDGEncoding());
if (pdgMother == 11 || pdgMother == 22)
return false;
}
int pdgMother = std::abs(fTrackingAction->geant4Track()->GetDefinition()->GetPDGEncoding());
if (pdgMother == 11 || pdgMother == 22)
return false;
G4double energy = track->GetKineticEnergy();
return energy < fEmax && fRegion == fastTrack.GetEnvelope();
return (energy < fEmax && fRegion == fastTrack.GetEnvelope());
}

void LowEnergyFastSimModel::DoIt(const G4FastTrack& fastTrack, G4FastStep& fastStep) {
fastStep.KillPrimaryTrack();
fastStep.SetPrimaryTrackPathLength(0.0);
double energy = fastTrack.GetPrimaryTrack()->GetKineticEnergy();
G4double energy = fastTrack.GetPrimaryTrack()->GetKineticEnergy();

const G4ThreeVector& pos = fastTrack.GetPrimaryTrack()->GetPosition();

G4double inPointEnergy = param.GetInPointEnergyFraction(energy) * energy;
G4double inPointEnergy = fParam.GetInPointEnergyFraction(energy) * energy;

// take into account positron annihilation (not included in in-point)
if (-11 == fastTrack.GetPrimaryTrack()->GetDefinition()->GetPDGEncoding())
energy += twomass;

const G4ThreeVector& momDir = fastTrack.GetPrimaryTrack()->GetMomentumDirection();
const G4ThreeVector& ortho = momDir.orthogonal();
const G4ThreeVector& cross = momDir.cross(ortho);

// in point energy deposition
GFlashEnergySpot spot;
Expand All @@ -59,17 +63,19 @@ void LowEnergyFastSimModel::DoIt(const G4FastTrack& fastTrack, G4FastStep& fastS

// tail energy deposition
G4double etail = energy - inPointEnergy;
const G4int nspots = int(etail) + 1;
const G4int nspots = (G4int)(etail) + 1;
const G4double tailEnergy = etail / (G4double)nspots;
for (G4int i = 0; i < nspots; ++i) {
const G4double radius = param.GetRadius(energy);
const G4double z = param.GetZ();
const G4double r = fParam.GetRadius(energy);
const G4double z = fParam.GetZ();

const G4double phi = CLHEP::twopi * G4UniformRand();
const G4ThreeVector tailPos = pos + z * momDir + radius * std::cos(phi) * ortho + radius * std::sin(phi) * cross;
fTailPos.set(r * std::cos(phi), r * std::sin(phi), z);
fTailPos.rotateUz(momDir);
fTailPos += pos;

spot.SetEnergy(tailEnergy);
spot.SetPosition(tailPos);
spot.SetPosition(fTailPos);
fHitMaker.make(&spot, &fastTrack);
}
}