Skip to content

Commit

Permalink
Merge pull request #42223 from civanch/fix_scripts_mctruth
Browse files Browse the repository at this point in the history
SIM: fixed scripts and mctruth
  • Loading branch information
cmsbuild authored Jul 11, 2023
2 parents 7d21739 + 591a0f6 commit 8e8aabe
Show file tree
Hide file tree
Showing 8 changed files with 40 additions and 60 deletions.
2 changes: 1 addition & 1 deletion SimG4CMS/HcalTestBeam/python/TB2006Analysis_cfi.py
Original file line number Diff line number Diff line change
Expand Up @@ -40,7 +40,7 @@ def testbeam2006(process):

process.common_maximum_time.MaxTrackTime = cms.double(1000.0)
process.common_maximum_time.MaxTimeNames = cms.vstring()
process.common_maximum_time.MaxTrackTimes = cms.vstring()
process.common_maximum_time.MaxTrackTimes = cms.vdouble()
process.common_maximum_time.DeadRegions = cms.vstring()

process.g4SimHits.NonBeamEvent = True
Expand Down
2 changes: 1 addition & 1 deletion SimG4Core/Application/python/g4SimHits_cfi.py
Original file line number Diff line number Diff line change
Expand Up @@ -706,7 +706,7 @@

from Configuration.Eras.Modifier_phase2_common_cff import phase2_common
phase2_common.toModify(g4SimHits,
OnlySDs = ['ZdcSensitiveDetector', 'TotemT2ScintSensitiveDetector', 'TotemSensitiveDetector', 'RomanPotSensitiveDetector', 'PLTSensitiveDetector', 'MuonSensitiveDetector', 'MtdSensitiveDetector', 'BCM1FSensitiveDetector', 'EcalSensitiveDetector', 'CTPPSSensitiveDetector', 'HGCalSensitiveDetector', 'BSCSensitiveDetector', 'CTPPSDiamondSensitiveDetector', 'FP420SensitiveDetector', 'BHMSensitiveDetector', 'HFNoseSensitiveDetector', 'HGCScintillatorSensitiveDetector', 'CastorSensitiveDetector', 'CaloTrkProcessing', 'HcalSensitiveDetector', 'TkAccumulatingSensitiveDetector'],
OnlySDs = ['ZdcSensitiveDetector', 'TotemT2ScintSensitiveDetector', 'TotemSensitiveDetector', 'RomanPotSensitiveDetector', 'PLTSensitiveDetector', 'MuonSensitiveDetector', 'MtdSensitiveDetector', 'BCM1FSensitiveDetector', 'EcalSensitiveDetector', 'CTPPSSensitiveDetector', 'HGCalSensitiveDetector', 'CTPPSDiamondSensitiveDetector', 'FP420SensitiveDetector', 'BHMSensitiveDetector', 'HFNoseSensitiveDetector', 'HGCScintillatorSensitiveDetector', 'CastorSensitiveDetector', 'CaloTrkProcessing', 'HcalSensitiveDetector', 'TkAccumulatingSensitiveDetector'],
LHCTransport = False,
MuonSD = dict(
HaveDemoChambers = False )
Expand Down
2 changes: 1 addition & 1 deletion SimG4Core/Notification/interface/TmpSimEvent.h
Original file line number Diff line number Diff line change
Expand Up @@ -14,7 +14,7 @@
class TmpSimEvent {
public:
TmpSimEvent();
virtual ~TmpSimEvent();
~TmpSimEvent();
void load(edm::SimTrackContainer& c) const;
void load(edm::SimVertexContainer& c) const;
unsigned int nTracks() const { return g4tracks_.size(); }
Expand Down
8 changes: 4 additions & 4 deletions SimG4Core/Notification/interface/TmpSimVertex.h
Original file line number Diff line number Diff line change
Expand Up @@ -7,20 +7,20 @@

class TmpSimVertex {
public:
TmpSimVertex(const math::XYZVectorD& ip, double it, int iv, unsigned int typ = 0)
: ilv_(ip), itime_(it), itrack_(iv), procType_(typ) {}
TmpSimVertex(const math::XYZVectorD& ip, double it, int iv, int typ = 0)
: ilv_(ip), itime_(it), itrack_(iv), ptype_(typ) {}
~TmpSimVertex() = default;
/// index of the parent (-1 if no parent)
const math::XYZVectorD& vertexPosition() const { return ilv_; }
double vertexGlobalTime() const { return itime_; }
int parentIndex() const { return itrack_; }
unsigned int processType() const { return procType_; }
int processType() const { return ptype_; }

private:
math::XYZVectorD ilv_;
double itime_;
int itrack_;
unsigned int procType_;
int ptype_;
};

#endif
51 changes: 24 additions & 27 deletions SimG4Core/Notification/interface/TrackWithHistory.h
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,6 @@

#include "G4Allocator.hh"

class G4VProcess;
class G4PrimaryParticle;
/** The part of the information about a SimTrack that we need from
* a G4Track
Expand All @@ -19,32 +18,35 @@ class TrackWithHistory {
* when some of the information may not available yet.
*/
TrackWithHistory(const G4Track *g4track, int pID);
TrackWithHistory(const G4PrimaryParticle *, int trackID, const math::XYZVectorD &pos, double time);
TrackWithHistory(const G4PrimaryParticle *, int trackID, const math::XYZVectorD &pos, const double time);
~TrackWithHistory() = default;

inline void *operator new(size_t);
inline void *operator new(std::size_t);
inline void operator delete(void *TrackWithHistory);

void setToBeSaved() { saved_ = true; }
int trackID() const { return trackID_; }
int particleID() const { return particleID_; }
int particleID() const { return pdgID_; }
int parentID() const { return parentID_; }
int genParticleID() const { return genParticleID_; }
int vertexID() const { return vertexID_; }
const math::XYZVectorD &momentum() const { return momentum_; }
double totalEnergy() const { return totalEnergy_; }
const math::XYZVectorD &vertexPosition() const { return vertexPosition_; }
double globalTime() const { return globalTime_; }
double localTime() const { return localTime_; }
double properTime() const { return properTime_; }
const G4VProcess *creatorProcess() const { return creatorProcess_; }
double weight() const { return weight_; }
int processType() const { return procType_; }
int getIDAtBoundary() const { return idAtBoundary_; }

void setTrackID(int i) { trackID_ = i; }
void setParentID(int i) { parentID_ = i; }
void setVertexID(int i) { vertexID_ = i; }
void setGenParticleID(int i) { genParticleID_ = i; }

double totalEnergy() const { return totalEnergy_; }
double time() const { return time_; }
double weight() const { return weight_; }
void setToBeSaved() { saved_ = true; }
bool storeTrack() const { return storeTrack_; }
bool saved() const { return saved_; }
bool crossedBoundary() const { return crossedBoundary_; }

const math::XYZVectorD &momentum() const { return momentum_; }
const math::XYZVectorD &vertexPosition() const { return vertexPosition_; }

// Boundary crossing variables
void setCrossedBoundaryPosMom(int id,
Expand All @@ -55,10 +57,8 @@ class TrackWithHistory {
positionAtBoundary_ = position;
momentumAtBoundary_ = momentum;
}
bool crossedBoundary() const { return crossedBoundary_; }
const math::XYZTLorentzVectorF &getPositionAtBoundary() const { return positionAtBoundary_; }
const math::XYZTLorentzVectorF &getMomentumAtBoundary() const { return momentumAtBoundary_; }
int getIDAtBoundary() const { return idAtBoundary_; }

// tracker surface
const math::XYZVectorD &trackerSurfacePosition() const { return tkSurfacePosition_; }
Expand All @@ -70,27 +70,24 @@ class TrackWithHistory {

private:
int trackID_;
int particleID_;
int pdgID_;
int parentID_;
int genParticleID_{-1};
int vertexID_{-1};
math::XYZVectorD momentum_;
int idAtBoundary_{-1};
int procType_{0};
double totalEnergy_;
math::XYZVectorD vertexPosition_;
double globalTime_;
double localTime_;
double properTime_;
const G4VProcess *creatorProcess_;
double time_; // lab system
double weight_;
bool storeTrack_;
bool saved_{false};

bool crossedBoundary_{false};
int idAtBoundary_{-1};
math::XYZVectorD momentum_;
math::XYZVectorD vertexPosition_;
math::XYZTLorentzVectorF positionAtBoundary_{math::XYZTLorentzVectorF(0.f, 0.f, 0.f, 0.f)};
math::XYZTLorentzVectorF momentumAtBoundary_{math::XYZTLorentzVectorF(0.f, 0.f, 0.f, 0.f)};
math::XYZVectorD tkSurfacePosition_{math::XYZVectorD(0., 0., 0.)};
math::XYZTLorentzVectorD tkSurfaceMomentum_{math::XYZTLorentzVectorD(0., 0., 0., 0.)};
bool storeTrack_{false};
bool saved_{false};
bool crossedBoundary_{false};
};

extern G4ThreadLocal G4Allocator<TrackWithHistory> *fpTrackWithHistoryAllocator;
Expand Down
7 changes: 1 addition & 6 deletions SimG4Core/Notification/src/SimTrackManager.cc
Original file line number Diff line number Diff line change
Expand Up @@ -189,12 +189,7 @@ int SimTrackManager::getOrCreateVertex(TrackWithHistory* trkH, int iParentID) {
}
}

unsigned int ptype = 0;
const G4VProcess* pr = trkH->creatorProcess();
if (nullptr != pr) {
ptype = pr->GetProcessSubType();
}
m_simEvent->add(new TmpSimVertex(trkH->vertexPosition(), trkH->globalTime(), parent, ptype));
m_simEvent->add(new TmpSimVertex(trkH->vertexPosition(), trkH->time(), parent, trkH->processType()));
m_vertexMap[parent].push_back(VertexPosition(m_nVertices, trkH->vertexPosition()));
++m_nVertices;
return (m_nVertices - 1);
Expand Down
4 changes: 2 additions & 2 deletions SimG4Core/Notification/src/TmpSimEvent.cc
Original file line number Diff line number Diff line change
Expand Up @@ -56,13 +56,13 @@ void TmpSimEvent::load(edm::SimVertexContainer& c) const {
TmpSimVertex* vtx = g4vertices_[i];
auto pos = vtx->vertexPosition();
math::XYZVectorD v3(pos.x() * invcm, pos.y() * invcm, pos.z() * invcm);
float t = vtx->vertexGlobalTime() / CLHEP::second;
float t = (float)(vtx->vertexGlobalTime() / CLHEP::second);
int iv = vtx->parentIndex();
// v3 = position in cm
// t = global time in second
// iv = index of the parent in the SimEvent SimTrack container (-1 if no parent)
SimVertex v = SimVertex(v3, t, iv, i);
v.setProcessType(vtx->processType());
v.setProcessType((unsigned int)vtx->processType());
v.setEventId(EncodedEventId(0));
c.push_back(v);
}
Expand Down
24 changes: 6 additions & 18 deletions SimG4Core/Notification/src/TrackWithHistory.cc
Original file line number Diff line number Diff line change
Expand Up @@ -14,21 +14,18 @@

G4ThreadLocal G4Allocator<TrackWithHistory>* fpTrackWithHistoryAllocator = nullptr;

//#define DEBUG

TrackWithHistory::TrackWithHistory(const G4Track* g4trk, int pID) {
trackID_ = g4trk->GetTrackID();
particleID_ = G4TrackToParticleID::particleID(g4trk);
pdgID_ = G4TrackToParticleID::particleID(g4trk);
parentID_ = pID;
auto mom = g4trk->GetMomentum();
momentum_ = math::XYZVectorD(mom.x(), mom.y(), mom.z());
totalEnergy_ = g4trk->GetTotalEnergy();
auto pos = g4trk->GetPosition();
vertexPosition_ = math::XYZVectorD(pos.x(), pos.y(), pos.z());
globalTime_ = g4trk->GetGlobalTime();
localTime_ = g4trk->GetLocalTime();
properTime_ = g4trk->GetProperTime();
creatorProcess_ = g4trk->GetCreatorProcess();
time_ = g4trk->GetGlobalTime();
auto p = g4trk->GetCreatorProcess();
procType_ = (nullptr != p) ? p->GetProcessSubType() : 0;
TrackInformation* trkinfo = static_cast<TrackInformation*>(g4trk->GetUserInformation());
storeTrack_ = trkinfo->storeTrack();
auto vgprimary = g4trk->GetDynamicParticle()->GetPrimaryParticle();
Expand All @@ -41,30 +38,21 @@ TrackWithHistory::TrackWithHistory(const G4Track* g4trk, int pID) {
// V.I. weight is computed in the same way as before
// without usage of G4Track::GetWeight()
weight_ = 10000 * genParticleID_;
#ifdef DEBUG
LogDebug("TrackInformation") << " TrackWithHistory : created history for " << trackID_ << " with mother "
<< parentID_;
#endif
}

TrackWithHistory::TrackWithHistory(const G4PrimaryParticle* ptr, int trackID, const math::XYZVectorD& pos, double time) {
trackID_ = trackID;
particleID_ = G4TrackToParticleID::particleID(ptr, trackID_);
pdgID_ = G4TrackToParticleID::particleID(ptr, trackID_);
parentID_ = trackID;
auto mom = ptr->GetMomentum();
momentum_ = math::XYZVectorD(mom.x(), mom.y(), mom.z());
totalEnergy_ = ptr->GetTotalEnergy();
vertexPosition_ = math::XYZVectorD(pos.x(), pos.y(), pos.z());
globalTime_ = localTime_ = properTime_ = time;
creatorProcess_ = nullptr;
time_ = time;
storeTrack_ = true;
auto priminfo = static_cast<GenParticleInfo*>(ptr->GetUserInformation());
if (nullptr != priminfo) {
genParticleID_ = priminfo->id();
}
weight_ = 10000. * genParticleID_;
#ifdef DEBUG
LogDebug("TrackInformation") << " TrackWithHistory : created history for " << trackID_ << " with mother "
<< parentID_;
#endif
}

0 comments on commit 8e8aabe

Please sign in to comment.