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

Geant4 Sensitive Detectors update 4: revised selection of particle type for shower libraries #22145

Merged
merged 5 commits into from Apr 22, 2018
Merged

Geant4 Sensitive Detectors update 4: revised selection of particle type for shower libraries #22145

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
3efa343
select particle type for shower libraries using utility
civanch Feb 7, 2018
4c2fe89
removed commented lines
civanch Feb 7, 2018
69a112e
added He3
civanch Feb 9, 2018
d98a2fa
optimized class
civanch Feb 12, 2018
5a60d16
fixed HF hits
civanch Feb 22, 2018
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: 0 additions & 4 deletions SimG4CMS/Calo/interface/CaloSD.h
View file
Open in desktop
Expand Up @@ -22,9 +22,6 @@

#include "FWCore/MessageLogger/interface/MessageLogger.h"

// To be replaced by something else
/* #include "Utilities/Notification/interface/TimerProxy.h" */

#include "G4VPhysicalVolume.hh"
#include "G4Track.hh"
#include "G4VGFlashSensitiveDetector.hh"
Expand Down Expand Up @@ -134,7 +131,6 @@ class CaloSD : public SensitiveCaloDetector,
double correctT;
double kmaxIon, kmaxNeutron, kmaxProton;

G4int emPDG, epPDG, gammaPDG;
bool forceSave;

private:
Expand Down
3 changes: 2 additions & 1 deletion SimG4CMS/Calo/interface/HCalSD.h
View file
Open in desktop
Expand Up @@ -97,9 +97,10 @@ class HCalSD : public CaloSD, public Observer<const BeginOfJob *> {
bool testNumber, neutralDensity, testNS_;
double birk1, birk2, birk3, betaThr;
bool useHF, useShowerLibrary, useParam, applyFidCut;
bool isEM;
double eminHitHB, eminHitHE, eminHitHO, eminHitHF;
double deliveredLumi;
G4int mumPDG, mupPDG, depth_;
G4int depth_;
std::vector<double> gpar;
std::vector<int> hfLevels;
std::vector<G4String> hfNames, fibreNames, matNames;
Expand Down
4 changes: 0 additions & 4 deletions SimG4CMS/Calo/interface/HFShowerLibrary.h
View file
Open in desktop
Expand Up @@ -76,10 +76,6 @@ class HFShowerLibrary {
double dphi, rMin, rMax;
std::vector<double> gpar;

int emPDG, epPDG, gammaPDG;
int pi0PDG, etaPDG, nuePDG, numuPDG, nutauPDG;
int anuePDG, anumuPDG, anutauPDG, geantinoPDG;

int npe;
HFShowerPhotonCollection pe;
HFShowerPhotonCollection* photo;
Expand Down
20 changes: 3 additions & 17 deletions SimG4CMS/Calo/src/CaloSD.cc
View file
Open in desktop
Expand Up @@ -6,6 +6,7 @@
#include "SimG4CMS/Calo/interface/CaloSD.h"
#include "SimDataFormats/SimHitMaker/interface/CaloSlaveSD.h"
#include "SimG4Core/Notification/interface/TrackInformation.h"
#include "SimG4Core/Notification/interface/G4TrackToParticleID.h"
#include "SimG4Core/Application/interface/EventAction.h"

#include "G4EventManager.hh"
Expand Down Expand Up @@ -109,11 +110,8 @@ bool CaloSD::ProcessHits(G4GFlashSpot* aSpot, G4TouchableHistory*) {

if (aSpot != nullptr) {
theTrack = const_cast<G4Track *>(aSpot->GetOriginatorTrack()->GetPrimaryTrack());
G4int particleCode = theTrack->GetDefinition()->GetPDGEncoding();

if (particleCode == emPDG ||
particleCode == epPDG ||
particleCode == gammaPDG ) {
if (G4TrackToParticleID::isGammaElectronPositron(theTrack)) {
edepositEM = aSpot->GetEnergySpot()->GetEnergy();
edepositHAD = 0.;
} else {
Expand Down Expand Up @@ -253,10 +251,7 @@ bool CaloSD::getStepInfo(G4Step* aStep) {
#endif
}

G4int particleCode = theTrack->GetDefinition()->GetPDGEncoding();
if (particleCode == emPDG ||
particleCode == epPDG ||
particleCode == gammaPDG ) {
if (G4TrackToParticleID::isGammaElectronPositron(theTrack)) {
edepositEM = getEnergyDeposit(aStep);
edepositHAD = 0.;
} else {
Expand Down Expand Up @@ -459,15 +454,6 @@ double CaloSD::getAttenuation(const G4Step* aStep, double birk1, double birk2, d
}

void CaloSD::update(const BeginOfRun *) {
G4ParticleTable * theParticleTable = G4ParticleTable::GetParticleTable();
G4String particleName;
emPDG = theParticleTable->FindParticle(particleName="e-")->GetPDGEncoding();
epPDG = theParticleTable->FindParticle(particleName="e+")->GetPDGEncoding();
gammaPDG = theParticleTable->FindParticle(particleName="gamma")->GetPDGEncoding();
#ifdef DebugLog
edm::LogInfo("CaloSim") << "CaloSD: Particle code for e- = " << emPDG
<< " for e+ = " << epPDG << " for gamma = " << gammaPDG;
#endif
initRun();
runInit = true;
}
Expand Down
109 changes: 49 additions & 60 deletions SimG4CMS/Calo/src/HCalSD.cc
View file
Open in desktop
Expand Up @@ -8,6 +8,7 @@
#include "SimG4CMS/Calo/interface/HFFibreFiducial.h"
#include "SimG4CMS/Calo/interface/HcalTestNS.h"
#include "SimG4Core/Notification/interface/TrackInformation.h"
#include "SimG4Core/Notification/interface/G4TrackToParticleID.h"
#include "DataFormats/HcalDetId/interface/HcalDetId.h"
#include "DetectorDescription/Core/interface/DDFilter.h"
#include "DetectorDescription/Core/interface/DDFilteredView.h"
Expand Down Expand Up @@ -296,7 +297,7 @@ HCalSD::HCalSD(const std::string& name, const DDCompactView & cpv,
edm::LogInfo("HcalSim") << "HCalSD: (" << i << ") " << matNames[i]
<< " pointer " << materials[i];

mumPDG = mupPDG = 0;
isEM = false;

if (useLayerWt) readWeightFromFile(file);

Expand Down Expand Up @@ -365,20 +366,20 @@ bool HCalSD::ProcessHits(G4Step * aStep, G4TouchableHistory * ) {
depth_ = (aStep->GetPreStepPoint()->GetTouchable()->GetReplicaNumber(0))%10;
const G4LogicalVolume* lv =
aStep->GetPreStepPoint()->GetPhysicalVolume()->GetLogicalVolume();
isEM = (G4TrackToParticleID::isGammaElectronPositron(aStep->GetTrack())) ? true : false;
G4String nameVolume = lv->GetName();
if (isItHF(aStep)) {
G4int parCode = aStep->GetTrack()->GetDefinition()->GetPDGEncoding();
double weight(1.0);
if (m_HFDarkening) {
G4ThreeVector hitPoint = aStep->GetPreStepPoint()->GetPosition();
double r = hitPoint.perp()/CLHEP::cm;
double z = std::abs(hitPoint.z())/CLHEP::cm;
double dose_acquired = 0.;
if (z>=HFDarkening::lowZLimit && z <= HFDarkening::upperZLimit) {
unsigned int hfZLayer = (int)((z - HFDarkening::lowZLimit)/5);
if (hfZLayer >= HFDarkening::upperZLimit) hfZLayer = (HFDarkening::upperZLimit-1);
if (z>=HFDarkening::lowZLimit && z <= HFDarkening::upperZLimit) {
unsigned int hfZLayer = (int)((z - HFDarkening::lowZLimit)/5);
if (hfZLayer >= HFDarkening::upperZLimit) hfZLayer = (HFDarkening::upperZLimit-1);
float normalized_lumi = m_HFDarkening->int_lumi(deliveredLumi);
for (int i = hfZLayer; i != HFDarkening::numberOfZLayers; ++i) {
for (int i = hfZLayer; i != HFDarkening::numberOfZLayers; ++i) {
dose_acquired = m_HFDarkening->dose(i,r);
weight *= m_HFDarkening->degradation(normalized_lumi*dose_acquired);
}
Expand All @@ -403,16 +404,18 @@ bool HCalSD::ProcessHits(G4Step * aStep, G4TouchableHistory * ) {
<< " hits afterParamS*";
#endif
} else {
bool notaMuon = true;
if (parCode == mupPDG || parCode == mumPDG ) notaMuon = false;
if (useShowerLibrary && notaMuon) {
// shower library is applied only for gamma, e+- or stable hadrons
// muons are tracked via HF, hyperons does not interact
if (useShowerLibrary && !G4TrackToParticleID::isMuon(aStep->GetTrack())) {
if(isEM || G4TrackToParticleID::isStableHadronIon(aStep->GetTrack())) {
getFromLibrary(aStep, weight);
}
#ifdef EDM_ML_DEBUG
LogDebug("HcalSim") << "HCalSD: Starts shower library from "
<< nameVolume << " for Track "
<< aStep->GetTrack()->GetTrackID() << " ("
<< aStep->GetTrack()->GetDefinition()->GetParticleName() << ")";
#endif
getFromLibrary(aStep, weight);
} else if (isItFibre(lv)) {
#ifdef EDM_ML_DEBUG
LogDebug("HcalSim") << "HCalSD: Hit at Fibre in " << nameVolume
Expand Down Expand Up @@ -608,13 +611,6 @@ void HCalSD::update(const BeginOfJob * job) {

void HCalSD::initRun() {
G4ParticleTable * theParticleTable = G4ParticleTable::GetParticleTable();
G4String particleName;
mumPDG = theParticleTable->FindParticle(particleName="mu-")->GetPDGEncoding();
mupPDG = theParticleTable->FindParticle(particleName="mu+")->GetPDGEncoding();
#ifdef EDM_ML_DEBUG
LogDebug("HcalSim") << "HCalSD: Particle code for mu- = " << mumPDG
<< " for mu+ = " << mupPDG;
#endif
if (showerLibrary) showerLibrary->initRun(theParticleTable,hcalConstants);
if (showerParam) showerParam->initRun(theParticleTable,hcalConstants);
if (hfshower) hfshower->initRun(theParticleTable,hcalConstants);
Expand Down Expand Up @@ -785,52 +781,50 @@ bool HCalSD::isItinFidVolume (const G4ThreeVector& hitPoint) {
void HCalSD::getFromLibrary (G4Step* aStep, double weight) {
preStepPoint = aStep->GetPreStepPoint();
theTrack = aStep->GetTrack();
int det = 5;
bool ok;
int det = 5;
bool ok(false);

std::vector<HFShowerLibrary::Hit> hits = showerLibrary->getHits(aStep, ok, weight, false);

double etrack = preStepPoint->GetKineticEnergy();
int primaryID = setTrackID(aStep);
if(!hits.empty()) {
double etrack = preStepPoint->GetKineticEnergy();
int primaryID = setTrackID(aStep);

// Reset entry point for new primary
posGlobal = preStepPoint->GetPosition();
resetForNewPrimary(posGlobal, etrack);
// Reset entry point for new primary
posGlobal = preStepPoint->GetPosition();
resetForNewPrimary(posGlobal, etrack);

G4int particleCode = theTrack->GetDefinition()->GetPDGEncoding();
if (particleCode==emPDG || particleCode==epPDG || particleCode==gammaPDG) {
edepositEM = 1.*GeV;
edepositHAD = 0.;
} else {
edepositEM = 0.;
edepositHAD = 1.*GeV;
}
if (isEM) {
edepositEM = 1.*GeV;
edepositHAD = 0.;
} else {
edepositEM = 0.;
edepositHAD = 1.*GeV;
}
#ifdef EDM_ML_DEBUG
edm::LogInfo("HcalSim") << "HCalSD::getFromLibrary " <<hits.size()
<< " hits for " << GetName() << " of " << primaryID
<< " with " << theTrack->GetDefinition()->GetParticleName()
<< " of " << preStepPoint->GetKineticEnergy()/GeV << " GeV";
edm::LogInfo("HcalSim") << "HCalSD::getFromLibrary " <<hits.size()
<< " hits for " << GetName() << " of " << primaryID
<< " with " << theTrack->GetDefinition()->GetParticleName()
<< " of " << preStepPoint->GetKineticEnergy()/GeV << " GeV";
#endif
for (unsigned int i=0; i<hits.size(); ++i) {
G4ThreeVector hitPoint = hits[i].position;
if (isItinFidVolume (hitPoint)) {
int depth = hits[i].depth;
double time = hits[i].time;
unsigned int unitID = setDetUnitId(det, hitPoint, depth);
currentID.setID(unitID, time, primaryID, 0);
for (unsigned int i=0; i<hits.size(); ++i) {
G4ThreeVector hitPoint = hits[i].position;
if (isItinFidVolume (hitPoint)) {
int depth = hits[i].depth;
double time = hits[i].time;
unsigned int unitID = setDetUnitId(det, hitPoint, depth);
currentID.setID(unitID, time, primaryID, 0);
#ifdef plotDebug
plotProfile(aStep, hitPoint, 1.0*GeV, time, depth);
bool emType = false;
if (particleCode==emPDG || particleCode==epPDG || particleCode==gammaPDG)
emType = true;
plotHF(hitPoint,emType);
plotProfile(aStep, hitPoint, 1.0*GeV, time, depth);
plotHF(hitPoint,isEM);
#endif

// check if it is in the same unit and timeslice as the previous one
if (currentID == previousID) {
updateHit(currentHit);
} else {
if (!checkHit()) currentHit = createNewHit();
// check if it is in the same unit and timeslice as the previous one
if (currentID == previousID) {
updateHit(currentHit);
} else {
if (!checkHit()) currentHit = createNewHit();
}
}
}
}
Expand All @@ -848,14 +842,12 @@ void HCalSD::getFromLibrary (G4Step* aStep, double weight) {
void HCalSD::hitForFibre (const G4Step* aStep, double weight) { // if not ParamShower

const G4StepPoint* preStepPoint = aStep->GetPreStepPoint();
const G4Track* theTrack = aStep->GetTrack();
int primaryID = setTrackID(aStep);

int det = 5;
std::vector<HFShower::Hit> hits = hfshower->getHits(aStep, weight);

G4int particleCode = theTrack->GetDefinition()->GetPDGEncoding();
if (particleCode==emPDG || particleCode==epPDG || particleCode==gammaPDG) {
if (isEM) {
edepositEM = 1.*GeV;
edepositHAD = 0.;
} else {
Expand All @@ -880,10 +872,7 @@ void HCalSD::hitForFibre (const G4Step* aStep, double weight) { // if not ParamS
currentID.setID(unitID, time, primaryID, 0);
#ifdef plotDebug
plotProfile(aStep, hitPoint, edepositEM, time, depth);
bool emType = false;
if (particleCode==emPDG || particleCode==epPDG || particleCode==gammaPDG)
emType = true;
plotHF(hitPoint,emType);
plotHF(hitPoint,isEM);
#endif
// check if it is in the same unit and timeslice as the previous one
if (currentID == previousID) {
Expand Down
51 changes: 19 additions & 32 deletions SimG4CMS/Calo/src/HFShowerLibrary.cc
View file
Open in desktop
Expand Up @@ -8,6 +8,7 @@
#include "DetectorDescription/Core/interface/DDFilter.h"
#include "DetectorDescription/Core/interface/DDFilteredView.h"
#include "DetectorDescription/Core/interface/DDValue.h"
#include "SimG4Core/Notification/interface/G4TrackToParticleID.h"

#include "FWCore/Utilities/interface/Exception.h"

Expand Down Expand Up @@ -119,9 +120,6 @@ HFShowerLibrary::HFShowerLibrary(const std::string & name, const DDCompactView &

fibre = new HFFibre(name, cpv, p);
photo = new HFShowerPhotonCollection;
emPDG = epPDG = gammaPDG = 0;
pi0PDG = etaPDG = nuePDG = numuPDG = nutauPDG= 0;
anuePDG= anumuPDG = anutauPDG = geantinoPDG = 0;
}

HFShowerLibrary::~HFShowerLibrary() {
Expand All @@ -135,30 +133,6 @@ void HFShowerLibrary::initRun(G4ParticleTable * theParticleTable,
HcalDDDSimConstants* hcons) {

if (fibre) fibre->initRun(hcons);

G4String parName;
emPDG = theParticleTable->FindParticle(parName="e-")->GetPDGEncoding();
epPDG = theParticleTable->FindParticle(parName="e+")->GetPDGEncoding();
gammaPDG = theParticleTable->FindParticle(parName="gamma")->GetPDGEncoding();
pi0PDG = theParticleTable->FindParticle(parName="pi0")->GetPDGEncoding();
etaPDG = theParticleTable->FindParticle(parName="eta")->GetPDGEncoding();
nuePDG = theParticleTable->FindParticle(parName="nu_e")->GetPDGEncoding();
numuPDG = theParticleTable->FindParticle(parName="nu_mu")->GetPDGEncoding();
nutauPDG= theParticleTable->FindParticle(parName="nu_tau")->GetPDGEncoding();
anuePDG = theParticleTable->FindParticle(parName="anti_nu_e")->GetPDGEncoding();
anumuPDG= theParticleTable->FindParticle(parName="anti_nu_mu")->GetPDGEncoding();
anutauPDG= theParticleTable->FindParticle(parName="anti_nu_tau")->GetPDGEncoding();
geantinoPDG= theParticleTable->FindParticle(parName="geantino")->GetPDGEncoding();
#ifdef DebugLog
edm::LogInfo("HFShower") << "HFShowerLibrary: Particle codes for e- = "
<< emPDG << ", e+ = " << epPDG << ", gamma = "
<< gammaPDG << ", pi0 = " << pi0PDG << ", eta = "
<< etaPDG << ", geantino = " << geantinoPDG
<< "\n nu_e = " << nuePDG << ", nu_mu = "
<< numuPDG << ", nu_tau = " << nutauPDG
<< ", anti_nu_e = " << anuePDG << ", anti_nu_mu = "
<< anumuPDG << ", anti_nu_tau = " << anutauPDG;
#endif

//Radius (minimum and maximum)
std::vector<double> rTable = hcons->getRTableHF();
Expand Down Expand Up @@ -198,6 +172,11 @@ std::vector<HFShowerLibrary::Hit> HFShowerLibrary::getHits(G4Step * aStep,
G4String partType = track->GetDefinition()->GetParticleName();
int parCode = track->GetDefinition()->GetPDGEncoding();

// VI: for ions use internally pdg code of alpha in order to keep
// consistency with previous simulation
if(track->GetDefinition()->IsGeneralIon()) { parCode = 1000020040; }


#ifdef DebugLog
G4ThreeVector localPos = preStepPoint->GetTouchable()->GetHistory()->GetTopTransform().TransformPoint(hitPoint);
double zoff = localPos.z() + 0.5*gpar[1];
Expand All @@ -214,7 +193,10 @@ std::vector<HFShowerLibrary::Hit> HFShowerLibrary::getHits(G4Step * aStep,
#endif

double tSlice = (postStepPoint->GetGlobalTime())/nanosecond;
double pin = preStepPoint->GetTotalEnergy();

// use kinetic energy for protons and ions
double pin = (track->GetDefinition()->GetBaryonNumber() > 0)
? preStepPoint->GetKineticEnergy() : preStepPoint->GetTotalEnergy();

return fillHits(hitPoint,momDir,parCode,pin,ok,weight,tSlice,onlyLong);
}
Expand All @@ -226,12 +208,17 @@ std::vector<HFShowerLibrary::Hit> HFShowerLibrary::fillHits(const G4ThreeVector

std::vector<HFShowerLibrary::Hit> hit;
ok = false;
if (parCode == pi0PDG || parCode == etaPDG || parCode == nuePDG ||
parCode == numuPDG || parCode == nutauPDG || parCode == anuePDG ||
parCode == anumuPDG || parCode == anutauPDG || parCode == geantinoPDG)
bool isEM = G4TrackToParticleID::isGammaElectronPositron(parCode);
// shower is built only for gamma, e+- and stable hadrons
if (!isEM && !G4TrackToParticleID::isStableHadron(parCode)) {
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

hi @civanch @fabiocos -following my questions from the release meeting - I would suggest to add some debug printout here to see in which cases these two stanzas of code are not equivalent in practice for some generic sample. It looks like a few events are enough to know the common examples.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hi @davidlange6 , I practically did this when was trying to clean up calorimeter sensitive detectors. Now this PR #21964 is closed, because we decided first to fix a weak place - selection of particles for SL. In #21984 I was doing exactly as you propose and achieve at the end all WF except 2 to be identical. Most part of problems come from Castor - secondary particles hit it even for typical ttbar events (primary are filtered out by eta cut) and for a hit random numbers are used. Also association between hits and MCtruth is sensitive to minor details of HF and Castor. What I also observed : even if random sequence is identical there are differences in results if hits are not absolutely identical.

For this PR I do not do the same but it is clear, that depending of WF it is various hyperons, which are now do not considered for shower libraries. I do not see a need to include them in the list, because they can safely decay and the products of their decay will be treated by SL. Unfortunately, for sampling of a hyperon decay random numbers are also spent.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Right - I don't mean to suggest a change - rather just an explicit check that its just hyperons on a handful of events... @civanch

return hit;
}
ok = true;

// remove low-energy component
const double threshold = 50*MeV;
if(pin < threshold) { return hit; }

double pz = momDir.z();
double zint = hitPoint.z();

Expand All @@ -244,7 +231,7 @@ std::vector<HFShowerLibrary::Hit> HFShowerLibrary::fillHits(const G4ThreeVector
double ctheta = cos(momDir.theta());
double stheta = sin(momDir.theta());

if (parCode == emPDG || parCode == epPDG || parCode == gammaPDG ) {
if (isEM) {
if (pin<pmom[nMomBin-1]) {
interpolate(0, pin);
} else {
Expand Down
4 changes: 0 additions & 4 deletions SimG4CMS/Forward/interface/CastorShowerLibrary.h
View file
Open in desktop
Expand Up @@ -70,10 +70,6 @@ class CastorShowerLibrary {

std::vector<double> pmom;

int emPDG, epPDG, gammaPDG;
int pi0PDG, etaPDG, nuePDG, numuPDG, nutauPDG;
int anuePDG, anumuPDG, anutauPDG, geantinoPDG;
int mumPDG, mupPDG;
// new variables (bins in eta and phi)
unsigned int nBinsE, nBinsEta, nBinsPhi;
unsigned int nEvtPerBinE, nEvtPerBinEta, nEvtPerBinPhi;
Expand Down