From d2f4ecb4269cddb42d773cd49be93e412d473bf6 Mon Sep 17 00:00:00 2001 From: Tobias Boeckh Date: Thu, 29 Jan 2026 13:44:53 +0100 Subject: [PATCH 1/5] Add variable to track if electron originates from primary electron --- Pinpoint/include/AnalysisManager.hh | 1 + Pinpoint/include/PixelHit.hh | 5 +++- Pinpoint/include/TrackInformation.hh | 4 +++- Pinpoint/src/AnalysisManager.cc | 3 +++ Pinpoint/src/PixelHit.cc | 4 ++-- Pinpoint/src/PixelSD.cc | 3 ++- Pinpoint/src/TrackInformation.cc | 5 +++- Pinpoint/src/TrackingAction.cc | 36 ++++++++++++++++++++++++++++ 8 files changed, 55 insertions(+), 6 deletions(-) diff --git a/Pinpoint/include/AnalysisManager.hh b/Pinpoint/include/AnalysisManager.hh index b054bed..2325215 100644 --- a/Pinpoint/include/AnalysisManager.hh +++ b/Pinpoint/include/AnalysisManager.hh @@ -189,6 +189,7 @@ class AnalysisManager { // std::vector fPixelFromPrimaryPizero; // std::vector fPixelFromFSLPizero; std::vector fPixelFromPrimaryLepton; + std::vector fPixelFromPrimaryEMShower; // Truth position of hit in x, y, z std::vector fPixelTruthX; diff --git a/Pinpoint/include/PixelHit.hh b/Pinpoint/include/PixelHit.hh index 590e250..6f84fe5 100644 --- a/Pinpoint/include/PixelHit.hh +++ b/Pinpoint/include/PixelHit.hh @@ -12,7 +12,7 @@ class PixelHit : public G4VHit { public: PixelHit() = default; - PixelHit(G4double edep, G4int rowID, G4int colID, G4int layerId, G4int trackID, G4int parentID, G4int pdgc, G4bool isPrim); + PixelHit(G4double edep, G4int rowID, G4int colID, G4int layerId, G4int trackID, G4int parentID, G4int pdgc, G4bool isPrim, G4bool isEMShower); PixelHit(const PixelHit&) = default; ~PixelHit() override = default; @@ -44,6 +44,7 @@ public: // void SetFromPrimaryPizero(G4bool fromPrimaryPizero) { fFromPrimaryPizero = fromPrimaryPizero; } // void SetFromFSLPizero(G4bool fromFSLPizero) { fFromFSLPizero = fromFSLPizero; } void SetFromPrimaryLepton(G4bool fromPrimaryLepton) { fFromPrimaryLepton = fromPrimaryLepton; } + void SetFromPrimaryEMShower(G4bool fromPrimaryEMShower) { fFromPrimaryEMShower = fromPrimaryEMShower; } // void SetTruthHitPos(G4ThreeVector pos) { fTruthHitPos = pos; } G4int GetPDGCode() const { return fPDGCode; } @@ -65,6 +66,7 @@ public: // G4bool GetFromPrimaryPizero() const { return fFromPrimaryPizero; } // G4bool GetFromFSLPizero() const { return fFromFSLPizero; } G4bool GetFromPrimaryLepton() const { return fFromPrimaryLepton; } + G4bool GetFromPrimaryEMShower() const { return fFromPrimaryEMShower; } private: G4int fTrackID = -1; @@ -76,6 +78,7 @@ private: G4int fLayerID = -1; // G4int fCharge = 0; G4bool fFromPrimaryLepton = false; + G4bool fFromPrimaryEMShower = false; // G4bool fFromPrimaryPizero = false; // G4bool fFromFSLPizero = false; // G4ThreeVector fTruthHitPos; diff --git a/Pinpoint/include/TrackInformation.hh b/Pinpoint/include/TrackInformation.hh index 59c4ae5..1a4abbc 100644 --- a/Pinpoint/include/TrackInformation.hh +++ b/Pinpoint/include/TrackInformation.hh @@ -24,9 +24,11 @@ class TrackInformation : public G4VUserTrackInformation inline void SetTrackIsFromPrimaryPizero(G4int i) {fFromPrimaryPizero = i;} inline void SetTrackIsFromFSLPizero(G4int i) {fFromFSLPizero = i;} inline void SetTrackIsFromPrimaryLepton(G4int i) {fFromPrimaryLepton = i;} + inline void SetTrackIsFromPrimaryEMShower(G4int i) {fFromPrimaryEMShower = i;} inline G4int IsTrackFromPrimaryPizero() const {return fFromPrimaryPizero;} inline G4int IsTrackFromFSLPizero() const {return fFromFSLPizero;} inline G4int IsTrackFromPrimaryLepton() const {return fFromPrimaryLepton;} + inline G4int IsTrackFromPrimaryEMShower() const {return fFromPrimaryEMShower;} inline void InsertHit(G4long hitIndex) { if (!fHitIndices) { @@ -40,7 +42,7 @@ class TrackInformation : public G4VUserTrackInformation G4int fFromPrimaryPizero; G4int fFromFSLPizero; G4int fFromPrimaryLepton; - + G4int fFromPrimaryEMShower; std::shared_ptr> fHitIndices; // std::vector* fHitIndices; //TODO: Make this a smart pointer? diff --git a/Pinpoint/src/AnalysisManager.cc b/Pinpoint/src/AnalysisManager.cc index b564339..525379a 100644 --- a/Pinpoint/src/AnalysisManager.cc +++ b/Pinpoint/src/AnalysisManager.cc @@ -198,6 +198,7 @@ void AnalysisManager::bookHitsTrees() // fPixelHitsTree->Branch("hit_fromPrimaryPizero", &fPixelFromPrimaryPizero); // fPixelHitsTree->Branch("hit_fromFSLPizero", &fPixelFromFSLPizero); fPixelHitsTree->Branch("hit_fromPrimaryLepton", &fPixelFromPrimaryLepton); + fPixelHitsTree->Branch("hit_fromPrimaryEMShower", &fPixelFromPrimaryEMShower); // if (fSaveTruthHits) // { @@ -310,6 +311,7 @@ void AnalysisManager::BeginOfEvent() // fPixelFromPrimaryPizero.clear(); // fPixelFromFSLPizero.clear(); fPixelFromPrimaryLepton.clear(); + fPixelFromPrimaryEMShower.clear(); // fPixelTruthX.clear(); // fPixelTruthY.clear(); // fPixelTruthZ.clear(); @@ -581,6 +583,7 @@ void AnalysisManager::FillHitsOutput() // fPixelFromPrimaryPizero.push_back(hit->GetFromPrimaryPizero()); // fPixelFromFSLPizero.push_back(hit->GetFromFSLPizero()); fPixelFromPrimaryLepton.push_back(hit->GetFromPrimaryLepton()); + fPixelFromPrimaryEMShower.push_back(hit->GetFromPrimaryEMShower()); // if (fSaveTruthHits) // { diff --git a/Pinpoint/src/PixelHit.cc b/Pinpoint/src/PixelHit.cc index 4d3f29e..b5d2735 100644 --- a/Pinpoint/src/PixelHit.cc +++ b/Pinpoint/src/PixelHit.cc @@ -41,7 +41,7 @@ void PixelHit::Print() << G4endl; } -PixelHit::PixelHit(G4double edep, G4int rowID, G4int colID, G4int layerID, G4int trackID, G4int parentID, G4int pdgc, G4bool isPrim) -: fEnergyDeposit{edep}, fRowID{rowID}, fColID{colID}, fLayerID{layerID}, fTrackID{trackID}, fParentID{parentID}, fPDGCode{pdgc}, fFromPrimaryLepton{isPrim} +PixelHit::PixelHit(G4double edep, G4int rowID, G4int colID, G4int layerID, G4int trackID, G4int parentID, G4int pdgc, G4bool isPrim, G4bool isEMShower) +: fEnergyDeposit{edep}, fRowID{rowID}, fColID{colID}, fLayerID{layerID}, fTrackID{trackID}, fParentID{parentID}, fPDGCode{pdgc}, fFromPrimaryLepton{isPrim}, fFromPrimaryEMShower{isEMShower} { } \ No newline at end of file diff --git a/Pinpoint/src/PixelSD.cc b/Pinpoint/src/PixelSD.cc index 6d0c536..c1013ab 100644 --- a/Pinpoint/src/PixelSD.cc +++ b/Pinpoint/src/PixelSD.cc @@ -86,6 +86,7 @@ G4bool PixelHitAccumulator::AddHit(G4Step* step) const auto* info =static_cast(track->GetUserInformation()); G4bool fromPrimaryLepton = info && info->IsTrackFromPrimaryLepton() || parentID ==0; fromPrimaryLepton = fromPrimaryLepton && (std::abs(pdgid) == 11 || std::abs(pdgid) == 13 || std::abs(pdgid) == 15); + G4bool fromPrimaryEMShower = info && info->IsTrackFromPrimaryEMShower(); assert(rowID < fNPixelsY); assert(colID < fNPixelsX); @@ -106,7 +107,7 @@ G4bool PixelHitAccumulator::AddHit(G4Step* step) } else { fPixelHits.push_back( new PixelHit(edep, rowID, colID, layerID, - trackID, parentID, pdgid, fromPrimaryLepton) + trackID, parentID, pdgid, fromPrimaryLepton, fromPrimaryEMShower) ); } diff --git a/Pinpoint/src/TrackInformation.cc b/Pinpoint/src/TrackInformation.cc index 8f4f8a2..2e2df1f 100644 --- a/Pinpoint/src/TrackInformation.cc +++ b/Pinpoint/src/TrackInformation.cc @@ -9,6 +9,7 @@ TrackInformation::TrackInformation() fFromPrimaryPizero = 0; fFromFSLPizero = 0; fFromPrimaryLepton = 0; + fFromPrimaryEMShower = 0; fHitIndices = std::make_shared>(); } @@ -18,6 +19,7 @@ TrackInformation::TrackInformation(const G4Track* aTrack) fFromPrimaryPizero = 0; fFromFSLPizero = 0; fFromPrimaryLepton = 0; + fFromPrimaryEMShower = 0; fHitIndices = std::make_shared>(); } @@ -30,7 +32,7 @@ ::operator =(const TrackInformation& aTrackInfo) fFromPrimaryPizero = aTrackInfo.fFromPrimaryPizero; fFromFSLPizero = aTrackInfo.fFromFSLPizero; fFromPrimaryLepton = aTrackInfo.fFromPrimaryLepton; - + fFromPrimaryEMShower = aTrackInfo.fFromPrimaryEMShower; return *this; } @@ -39,5 +41,6 @@ void TrackInformation::Print() const G4cout << "Is from primary pizero " << fFromPrimaryPizero << G4endl; G4cout << "Is from final state lepton decay pizero " << fFromFSLPizero << G4endl; G4cout << "Is from primary lepton (tau or muon) " << fFromPrimaryLepton << G4endl; + G4cout << "Is from primary EM shower " << fFromPrimaryEMShower << G4endl; } diff --git a/Pinpoint/src/TrackingAction.cc b/Pinpoint/src/TrackingAction.cc index 900d7fa..f27147e 100644 --- a/Pinpoint/src/TrackingAction.cc +++ b/Pinpoint/src/TrackingAction.cc @@ -71,6 +71,42 @@ void TrackingAction::PostUserTrackingAction(const G4Track* aTrack) } } + if (aTrack->GetTrackID()==1 && abs(aTrack->GetParticleDefinition()->GetPDGEncoding())==11) + { + const auto* info =static_cast(aTrack->GetUserInformation()); + TrackInformation* newInfo = new TrackInformation(); + if (info) { + newInfo->SetTrackIsFromPrimaryPizero(info->IsTrackFromPrimaryPizero()); + newInfo->SetTrackIsFromFSLPizero(info->IsTrackFromFSLPizero()); + newInfo->SetTrackIsFromPrimaryLepton(info->IsTrackFromPrimaryLepton()); + } + newInfo->SetTrackIsFromPrimaryEMShower(1); + aTrack->SetUserInformation(newInfo); + } + + const auto* info =static_cast(aTrack->GetUserInformation()); + G4bool fromPrimaryEMShower = info && info->IsTrackFromPrimaryEMShower(); + + if (fromPrimaryEMShower && + (abs(aTrack->GetParticleDefinition()->GetPDGEncoding())==11 || + abs(aTrack->GetParticleDefinition()->GetPDGEncoding())==22)) + { + G4TrackVector* secondaries = fpTrackingManager->GimmeSecondaries(); + if (secondaries) + { + size_t nSeco = secondaries->size(); + if (nSeco>0) + { + for (size_t i=0; iSetTrackIsFromPrimaryEMShower(1); + (*secondaries)[i]->SetUserInformation(info); + } + } + } + } + if (aTrack->GetParentID()==1 && aTrack->GetCreatorProcess()->GetProcessName()=="Decay") { // in case of tau decay pizero From 0c3898d96f09fb21f44ca879ca1bd244625d8399 Mon Sep 17 00:00:00 2001 From: Tobias Boeckh Date: Fri, 6 Feb 2026 16:46:57 +0100 Subject: [PATCH 2/5] Fix G4Exceptions: remove slow particles in field region which may curl --- Pinpoint/src/DetectorConstruction.cc | 51 ++++++++++--------- .../generators/GFaserGeneratorMessenger.cc | 2 +- 2 files changed, 28 insertions(+), 25 deletions(-) diff --git a/Pinpoint/src/DetectorConstruction.cc b/Pinpoint/src/DetectorConstruction.cc index ef17363..e286b36 100644 --- a/Pinpoint/src/DetectorConstruction.cc +++ b/Pinpoint/src/DetectorConstruction.cc @@ -140,8 +140,8 @@ G4VPhysicalVolume* DetectorConstruction::Construct() fNLayers = maxLayers; } auto detectorThickness = fNLayers * fLayerThickness; - auto worldSizeX = 1.2 * fDetectorWidth; - auto worldSizeY = 1.2 * fDetectorHeight; + auto worldSizeX = std::max(1.2 * fDetectorWidth, 2.4 * fOuterRadius); + auto worldSizeY = std::max(1.2 * fDetectorHeight, 2.4 * fOuterRadius); auto worldSizeZ = 1.2 * (detectorThickness + fTracker3Position); // Get materials @@ -335,35 +335,36 @@ G4VPhysicalVolume* DetectorConstruction::Construct() // solid cylinders (0 to outerRadius) and air-filled bore (0 to innerRadius) G4Material* sm2co17 = G4Material::GetMaterial("Sm2Co17"); - auto longMagnetS = new G4Tubs("Magnet0", 0., fOuterRadius, 0.5 * fLongMagnetLength, 0., 2*M_PI); - auto shortMagnetS = new G4Tubs("Magnet12", 0., fOuterRadius, 0.5 * fShortMagnetLength, 0., 2*M_PI); + // Magnets as hollow tubes — field regions placed as siblings in worldLV + // to avoid mother-daughter shared-surface navigation issues + auto longMagnetS = new G4Tubs("Magnet0", fInnerRadius, fOuterRadius, 0.5 * fLongMagnetLength, 0., 2*M_PI); + auto shortMagnetS = new G4Tubs("Magnet12", fInnerRadius, fOuterRadius, 0.5 * fShortMagnetLength, 0., 2*M_PI); - G4double fieldRadius = fInnerRadius; - auto longFieldS = new G4Tubs("FieldRegion0", 0., fieldRadius, 0.5 * fLongMagnetLength, 0., 2*M_PI); - auto shortFieldS = new G4Tubs("FieldRegion12", 0., fieldRadius, 0.5 * fShortMagnetLength, 0., 2*M_PI); + auto longFieldS = new G4Tubs("FieldRegion0", 0., fInnerRadius, 0.5 * fLongMagnetLength, 0., 2*M_PI); + auto shortFieldS = new G4Tubs("FieldRegion12", 0., fInnerRadius, 0.5 * fShortMagnetLength, 0., 2*M_PI); // Magnet 0 - auto magnet0LV = new G4LogicalVolume(longMagnetS, sm2co17, "Magnet0"); - new G4PVPlacement(nullptr, G4ThreeVector(0., 0., fMagnet0Position), magnet0LV, "Magnet0", worldLV, false, 0, fCheckOverlaps); - magnet0LV->SetVisAttributes(G4VisAttributes(G4Colour(1.0, 0.0, 1.0, 0.5))); // Magenta, semi-transparent + // auto magnet0LV = new G4LogicalVolume(longMagnetS, sm2co17, "Magnet0"); + // new G4PVPlacement(nullptr, G4ThreeVector(0., 0., fMagnet0Position), magnet0LV, "Magnet0", worldLV, false, 0, fCheckOverlaps); + // magnet0LV->SetVisAttributes(G4VisAttributes(G4Colour(1.0, 0.0, 1.0, 0.5))); // Magenta, semi-transparent auto fieldRegion0LV = new G4LogicalVolume(longFieldS, worldMaterial, "FieldRegion0"); - new G4PVPlacement(nullptr, G4ThreeVector(0., 0., 0.), fieldRegion0LV, "FieldRegion0", magnet0LV, false, 0, fCheckOverlaps); + new G4PVPlacement(nullptr, G4ThreeVector(0., 0., fMagnet0Position), fieldRegion0LV, "FieldRegion0", worldLV, false, 0, fCheckOverlaps); fieldRegion0LV->SetVisAttributes(G4VisAttributes::GetInvisible()); // Magnet 1 - auto magnet1LV = new G4LogicalVolume(shortMagnetS, sm2co17, "Magnet1"); - new G4PVPlacement(nullptr, G4ThreeVector(0., 0., fMagnet1Position), magnet1LV, "Magnet1", worldLV, false, 1, fCheckOverlaps); - magnet1LV->SetVisAttributes(G4VisAttributes(G4Colour(1.0, 0.0, 1.0, 0.5))); // Magenta, semi-transparent + // auto magnet1LV = new G4LogicalVolume(shortMagnetS, sm2co17, "Magnet1"); + // new G4PVPlacement(nullptr, G4ThreeVector(0., 0., fMagnet1Position), magnet1LV, "Magnet1", worldLV, false, 1, fCheckOverlaps); + // magnet1LV->SetVisAttributes(G4VisAttributes(G4Colour(1.0, 0.0, 1.0, 0.5))); // Magenta, semi-transparent auto fieldRegion1LV = new G4LogicalVolume(shortFieldS, worldMaterial, "FieldRegion1"); - new G4PVPlacement(nullptr, G4ThreeVector(0., 0., 0.), fieldRegion1LV, "FieldRegion1", magnet1LV, false, 1, fCheckOverlaps); + new G4PVPlacement(nullptr, G4ThreeVector(0., 0., fMagnet1Position), fieldRegion1LV, "FieldRegion1", worldLV, false, 1, fCheckOverlaps); fieldRegion1LV->SetVisAttributes(G4VisAttributes::GetInvisible()); // Magnet 2 - auto magnet2LV = new G4LogicalVolume(shortMagnetS, sm2co17, "Magnet2"); - new G4PVPlacement(nullptr, G4ThreeVector(0., 0., fMagnet2Position), magnet2LV, "Magnet2", worldLV, false, 2, fCheckOverlaps); - magnet2LV->SetVisAttributes(G4VisAttributes(G4Colour(1.0, 0.0, 1.0, 0.5))); // Magenta, semi-transparent + // auto magnet2LV = new G4LogicalVolume(shortMagnetS, sm2co17, "Magnet2"); + // new G4PVPlacement(nullptr, G4ThreeVector(0., 0., fMagnet2Position), magnet2LV, "Magnet2", worldLV, false, 2, fCheckOverlaps); + // magnet2LV->SetVisAttributes(G4VisAttributes(G4Colour(1.0, 0.0, 1.0, 0.5))); // Magenta, semi-transparent auto fieldRegion2LV = new G4LogicalVolume(shortFieldS, worldMaterial, "FieldRegion2"); - new G4PVPlacement(nullptr, G4ThreeVector(0., 0., 0.), fieldRegion2LV, "FieldRegion2", magnet2LV, false, 2, fCheckOverlaps); + new G4PVPlacement(nullptr, G4ThreeVector(0., 0., fMagnet2Position), fieldRegion2LV, "FieldRegion2", worldLV, false, 2, fCheckOverlaps); fieldRegion2LV->SetVisAttributes(G4VisAttributes::GetInvisible()); G4cout << "Created solid cylindrical magnets:" << G4endl; @@ -373,11 +374,13 @@ G4VPhysicalVolume* DetectorConstruction::Construct() G4cout << " Magnet 2: z = " << fMagnet2Position/mm << " mm, length = " << fShortMagnetLength/mm << " mm" << G4endl; G4cout << " Material: Sm2Co17" << G4endl; - // Set step limits in magnetic field regions for accurate tracking - auto userLimits = new G4UserLimits(1 * mm); // Max step size 1mm - fieldRegion0LV->SetUserLimits(userLimits); - fieldRegion1LV->SetUserLimits(userLimits); - fieldRegion2LV->SetUserLimits(userLimits); + // Set step limits and minimum energy in magnetic field regions + auto fieldRegionUserLimits = new G4UserLimits(); + fieldRegionUserLimits->SetMaxAllowedStep(0.5 * mm); + fieldRegionUserLimits->SetUserMinEkine(10.0 * MeV); + fieldRegion0LV->SetUserLimits(fieldRegionUserLimits); + fieldRegion1LV->SetUserLimits(fieldRegionUserLimits); + fieldRegion2LV->SetUserLimits(fieldRegionUserLimits); // FASER spectrometer tracking layers (use single layer per station) auto trackerS = new G4Box("Tracker", 0.5 * fTrackerSize, 0.5 * fTrackerSize, 0.5 * fSiliconThickness); diff --git a/Pinpoint/src/generators/GFaserGeneratorMessenger.cc b/Pinpoint/src/generators/GFaserGeneratorMessenger.cc index 3f1b96c..8c515f4 100644 --- a/Pinpoint/src/generators/GFaserGeneratorMessenger.cc +++ b/Pinpoint/src/generators/GFaserGeneratorMessenger.cc @@ -11,7 +11,7 @@ GFaserGeneratorMessenger::GFaserGeneratorMessenger(GFaserGenerator* action) : fGFaserAction(action) { - fGFaserGeneratorDir = new G4UIdirectory("/gen/gfaser"); + fGFaserGeneratorDir = new G4UIdirectory("/gen/gfaser/"); fGFaserGeneratorDir->SetGuidance("gfaser generator control"); fInputFileCmd = new G4UIcmdWithAString("/gen/gfaser/inputFile", this); From 289836a2933e12eac56217e3ea7eb379bc072183 Mon Sep 17 00:00:00 2001 From: Tobias Boeckh Date: Fri, 6 Feb 2026 18:06:47 +0100 Subject: [PATCH 3/5] Fix vertex z position is fUseFixedZPosition is false --- Pinpoint/src/generators/GFaserGenerator.cc | 20 ++++++++++++++++++++ 1 file changed, 20 insertions(+) diff --git a/Pinpoint/src/generators/GFaserGenerator.cc b/Pinpoint/src/generators/GFaserGenerator.cc index 2350d6e..1fd9323 100644 --- a/Pinpoint/src/generators/GFaserGenerator.cc +++ b/Pinpoint/src/generators/GFaserGenerator.cc @@ -151,6 +151,13 @@ void GFaserGenerator::GeneratePrimaries(G4Event* event) if (fUseFixedZPosition) { fVz = GenerateRandomZVertex(fLayerId); } + else { + auto *rm = G4RunManager::GetRunManager(); + auto det = (DetectorConstruction*) (rm->GetUserDetectorConstruction()); + G4int maxLayerIndex = det->GetNumberOfLayers(); + G4int randomLayer = (G4int)(G4UniformRand() * (maxLayerIndex)); + fVz = GenerateRandomZVertex(randomLayer); + } auto *runManager = G4RunManager::GetRunManager(); auto detector = (DetectorConstruction*) (runManager->GetUserDetectorConstruction()); @@ -268,6 +275,19 @@ G4double GFaserGenerator::GenerateRandomZVertex(G4int layerIndex) const { G4double layerThickness = detector->GetLayerThickness(); G4double tungstenThickness = detector->GetTungstenThickness(); G4double z = layerIndex * layerThickness + tungstenThickness * G4UniformRand(); + + if (detector->GetSimFlag() >= 0 && G4UniformRand() > 0.5) { + G4double boxThickness = detector->GetBoxThickness(); + G4double siliconThickness = detector->GetSiliconThickness(); + // If `GetSimFlag() >= 0` there are two tungsten blocks per layer: [Tungsten | Box | Silicon | Tungsten | Scintillator] + // Randomly chose in which tungsten block the neutrino interaction happens and adjust the z vertex accordingly. + z += tungstenThickness + boxThickness + siliconThickness; + } + G4cout << "LayerIndex: " << layerIndex << ", LayerThickness: " << layerThickness/mm << " mm, TungstenThickness: " << tungstenThickness/mm << " mm, " + << " SimFlag: " << detector->GetSimFlag() << ", boxThickness: " << detector->GetBoxThickness()/mm << " mm, " + << "siliconThickness: " << detector->GetSiliconThickness()/mm << " mm" + << G4endl; + G4cout << "Generated random Z vertex at: " << z << " mm = " << z/m << " m in layer " << layerIndex << G4endl; return z/m; // convert to meters From ca74f85f6ed0cc65655e3e0fdc08984fcc72ae57 Mon Sep 17 00:00:00 2001 From: Tobias Boeckh Date: Fri, 6 Feb 2026 18:37:49 +0100 Subject: [PATCH 4/5] Fix paranthesis and add more charmed hadrons --- Pinpoint/src/Py8DecayerPhysics.cc | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/Pinpoint/src/Py8DecayerPhysics.cc b/Pinpoint/src/Py8DecayerPhysics.cc index 67b264f..ab3cd05 100644 --- a/Pinpoint/src/Py8DecayerPhysics.cc +++ b/Pinpoint/src/Py8DecayerPhysics.cc @@ -94,7 +94,14 @@ void Py8DecayerPhysics::ConstructProcess() if (std::abs(particle->GetPDGEncoding()) == 15 || std::abs(particle->GetPDGEncoding()) == 521 || std::abs(particle->GetPDGEncoding()) == 411 || - std::abs(particle->GetPDGEncoding() == 421)) + std::abs(particle->GetPDGEncoding()) == 421 || + std::abs(particle->GetPDGEncoding()) == 431 || + std::abs(particle->GetPDGEncoding()) == 4122 || + std::abs(particle->GetPDGEncoding()) == 4112 || + std::abs(particle->GetPDGEncoding()) == 4212 || + std::abs(particle->GetPDGEncoding()) == 4222 || + std::abs(particle->GetPDGEncoding()) == 4132 || + std::abs(particle->GetPDGEncoding()) == 4232) { if (particle->GetDecayTable()) { delete particle->GetDecayTable(); From 4e56189344a4e64b3f12a485a3f494ca7b4eb57f Mon Sep 17 00:00:00 2001 From: Tobias Boeckh Date: Fri, 6 Feb 2026 19:38:10 +0100 Subject: [PATCH 5/5] Add PixelFromCharmedHadron var to tag charm decays (nu->X+c) --- Pinpoint/include/AnalysisManager.hh | 1 + Pinpoint/include/PixelHit.hh | 5 ++- Pinpoint/include/TrackInformation.hh | 3 ++ Pinpoint/src/AnalysisManager.cc | 3 ++ Pinpoint/src/PixelHit.cc | 4 +-- Pinpoint/src/PixelSD.cc | 3 +- Pinpoint/src/TrackInformation.cc | 4 +++ Pinpoint/src/TrackingAction.cc | 50 ++++++++++++++++++++++++++++ 8 files changed, 69 insertions(+), 4 deletions(-) diff --git a/Pinpoint/include/AnalysisManager.hh b/Pinpoint/include/AnalysisManager.hh index 403b9a6..03ff449 100644 --- a/Pinpoint/include/AnalysisManager.hh +++ b/Pinpoint/include/AnalysisManager.hh @@ -194,6 +194,7 @@ class AnalysisManager { // std::vector fPixelFromFSLPizero; std::vector fPixelFromPrimaryLepton; std::vector fPixelFromPrimaryEMShower; + std::vector fPixelFromCharmedHadron; // Truth position of hit in x, y, z std::vector fPixelTruthX; diff --git a/Pinpoint/include/PixelHit.hh b/Pinpoint/include/PixelHit.hh index 6f84fe5..ca53b2f 100644 --- a/Pinpoint/include/PixelHit.hh +++ b/Pinpoint/include/PixelHit.hh @@ -12,7 +12,7 @@ class PixelHit : public G4VHit { public: PixelHit() = default; - PixelHit(G4double edep, G4int rowID, G4int colID, G4int layerId, G4int trackID, G4int parentID, G4int pdgc, G4bool isPrim, G4bool isEMShower); + PixelHit(G4double edep, G4int rowID, G4int colID, G4int layerId, G4int trackID, G4int parentID, G4int pdgc, G4bool isPrim, G4bool isEMShower, G4bool isFromCharmed); PixelHit(const PixelHit&) = default; ~PixelHit() override = default; @@ -45,6 +45,7 @@ public: // void SetFromFSLPizero(G4bool fromFSLPizero) { fFromFSLPizero = fromFSLPizero; } void SetFromPrimaryLepton(G4bool fromPrimaryLepton) { fFromPrimaryLepton = fromPrimaryLepton; } void SetFromPrimaryEMShower(G4bool fromPrimaryEMShower) { fFromPrimaryEMShower = fromPrimaryEMShower; } + void SetFromCharmedHadron(G4bool fromCharmedHadron) { fFromCharmedHadron = fromCharmedHadron; } // void SetTruthHitPos(G4ThreeVector pos) { fTruthHitPos = pos; } G4int GetPDGCode() const { return fPDGCode; } @@ -67,6 +68,7 @@ public: // G4bool GetFromFSLPizero() const { return fFromFSLPizero; } G4bool GetFromPrimaryLepton() const { return fFromPrimaryLepton; } G4bool GetFromPrimaryEMShower() const { return fFromPrimaryEMShower; } + G4bool GetFromCharmedHadron() const { return fFromCharmedHadron; } private: G4int fTrackID = -1; @@ -79,6 +81,7 @@ private: // G4int fCharge = 0; G4bool fFromPrimaryLepton = false; G4bool fFromPrimaryEMShower = false; + G4bool fFromCharmedHadron = false; // G4bool fFromPrimaryPizero = false; // G4bool fFromFSLPizero = false; // G4ThreeVector fTruthHitPos; diff --git a/Pinpoint/include/TrackInformation.hh b/Pinpoint/include/TrackInformation.hh index 1a4abbc..d4db35d 100644 --- a/Pinpoint/include/TrackInformation.hh +++ b/Pinpoint/include/TrackInformation.hh @@ -25,10 +25,12 @@ class TrackInformation : public G4VUserTrackInformation inline void SetTrackIsFromFSLPizero(G4int i) {fFromFSLPizero = i;} inline void SetTrackIsFromPrimaryLepton(G4int i) {fFromPrimaryLepton = i;} inline void SetTrackIsFromPrimaryEMShower(G4int i) {fFromPrimaryEMShower = i;} + inline void SetTrackIsFromCharmedHadron(G4int i) {fFromCharmedHadron = i;} inline G4int IsTrackFromPrimaryPizero() const {return fFromPrimaryPizero;} inline G4int IsTrackFromFSLPizero() const {return fFromFSLPizero;} inline G4int IsTrackFromPrimaryLepton() const {return fFromPrimaryLepton;} inline G4int IsTrackFromPrimaryEMShower() const {return fFromPrimaryEMShower;} + inline G4int IsTrackFromCharmedHadron() const {return fFromCharmedHadron;} inline void InsertHit(G4long hitIndex) { if (!fHitIndices) { @@ -43,6 +45,7 @@ class TrackInformation : public G4VUserTrackInformation G4int fFromFSLPizero; G4int fFromPrimaryLepton; G4int fFromPrimaryEMShower; + G4int fFromCharmedHadron; std::shared_ptr> fHitIndices; // std::vector* fHitIndices; //TODO: Make this a smart pointer? diff --git a/Pinpoint/src/AnalysisManager.cc b/Pinpoint/src/AnalysisManager.cc index 49a499a..6d34ba2 100644 --- a/Pinpoint/src/AnalysisManager.cc +++ b/Pinpoint/src/AnalysisManager.cc @@ -201,6 +201,7 @@ void AnalysisManager::bookHitsTrees() // fPixelHitsTree->Branch("hit_fromFSLPizero", &fPixelFromFSLPizero); fPixelHitsTree->Branch("hit_fromPrimaryLepton", &fPixelFromPrimaryLepton); fPixelHitsTree->Branch("hit_fromPrimaryEMShower", &fPixelFromPrimaryEMShower); + fPixelHitsTree->Branch("hit_fromCharmedHadron", &fPixelFromCharmedHadron); // if (fSaveTruthHits) // { @@ -337,6 +338,7 @@ void AnalysisManager::BeginOfEvent() // fPixelFromFSLPizero.clear(); fPixelFromPrimaryLepton.clear(); fPixelFromPrimaryEMShower.clear(); + fPixelFromCharmedHadron.clear(); // fPixelTruthX.clear(); // fPixelTruthY.clear(); // fPixelTruthZ.clear(); @@ -612,6 +614,7 @@ void AnalysisManager::FillHitsOutput() // fPixelFromFSLPizero.push_back(hit->GetFromFSLPizero()); fPixelFromPrimaryLepton.push_back(hit->GetFromPrimaryLepton()); fPixelFromPrimaryEMShower.push_back(hit->GetFromPrimaryEMShower()); + fPixelFromCharmedHadron.push_back(hit->GetFromCharmedHadron()); // if (fSaveTruthHits) // { diff --git a/Pinpoint/src/PixelHit.cc b/Pinpoint/src/PixelHit.cc index b5d2735..2120cab 100644 --- a/Pinpoint/src/PixelHit.cc +++ b/Pinpoint/src/PixelHit.cc @@ -41,7 +41,7 @@ void PixelHit::Print() << G4endl; } -PixelHit::PixelHit(G4double edep, G4int rowID, G4int colID, G4int layerID, G4int trackID, G4int parentID, G4int pdgc, G4bool isPrim, G4bool isEMShower) -: fEnergyDeposit{edep}, fRowID{rowID}, fColID{colID}, fLayerID{layerID}, fTrackID{trackID}, fParentID{parentID}, fPDGCode{pdgc}, fFromPrimaryLepton{isPrim}, fFromPrimaryEMShower{isEMShower} +PixelHit::PixelHit(G4double edep, G4int rowID, G4int colID, G4int layerID, G4int trackID, G4int parentID, G4int pdgc, G4bool isPrim, G4bool isEMShower, G4bool isFromCharmed) +: fEnergyDeposit{edep}, fRowID{rowID}, fColID{colID}, fLayerID{layerID}, fTrackID{trackID}, fParentID{parentID}, fPDGCode{pdgc}, fFromPrimaryLepton{isPrim}, fFromPrimaryEMShower{isEMShower}, fFromCharmedHadron{isFromCharmed} { } \ No newline at end of file diff --git a/Pinpoint/src/PixelSD.cc b/Pinpoint/src/PixelSD.cc index c1013ab..2b21a14 100644 --- a/Pinpoint/src/PixelSD.cc +++ b/Pinpoint/src/PixelSD.cc @@ -87,6 +87,7 @@ G4bool PixelHitAccumulator::AddHit(G4Step* step) G4bool fromPrimaryLepton = info && info->IsTrackFromPrimaryLepton() || parentID ==0; fromPrimaryLepton = fromPrimaryLepton && (std::abs(pdgid) == 11 || std::abs(pdgid) == 13 || std::abs(pdgid) == 15); G4bool fromPrimaryEMShower = info && info->IsTrackFromPrimaryEMShower(); + G4bool fromCharmedHadron = info && info->IsTrackFromCharmedHadron(); assert(rowID < fNPixelsY); assert(colID < fNPixelsX); @@ -107,7 +108,7 @@ G4bool PixelHitAccumulator::AddHit(G4Step* step) } else { fPixelHits.push_back( new PixelHit(edep, rowID, colID, layerID, - trackID, parentID, pdgid, fromPrimaryLepton, fromPrimaryEMShower) + trackID, parentID, pdgid, fromPrimaryLepton, fromPrimaryEMShower, fromCharmedHadron) ); } diff --git a/Pinpoint/src/TrackInformation.cc b/Pinpoint/src/TrackInformation.cc index 2e2df1f..d526025 100644 --- a/Pinpoint/src/TrackInformation.cc +++ b/Pinpoint/src/TrackInformation.cc @@ -10,6 +10,7 @@ TrackInformation::TrackInformation() fFromFSLPizero = 0; fFromPrimaryLepton = 0; fFromPrimaryEMShower = 0; + fFromCharmedHadron = 0; fHitIndices = std::make_shared>(); } @@ -20,6 +21,7 @@ TrackInformation::TrackInformation(const G4Track* aTrack) fFromFSLPizero = 0; fFromPrimaryLepton = 0; fFromPrimaryEMShower = 0; + fFromCharmedHadron = 0; fHitIndices = std::make_shared>(); } @@ -33,6 +35,7 @@ ::operator =(const TrackInformation& aTrackInfo) fFromFSLPizero = aTrackInfo.fFromFSLPizero; fFromPrimaryLepton = aTrackInfo.fFromPrimaryLepton; fFromPrimaryEMShower = aTrackInfo.fFromPrimaryEMShower; + fFromCharmedHadron = aTrackInfo.fFromCharmedHadron; return *this; } @@ -42,5 +45,6 @@ void TrackInformation::Print() const G4cout << "Is from final state lepton decay pizero " << fFromFSLPizero << G4endl; G4cout << "Is from primary lepton (tau or muon) " << fFromPrimaryLepton << G4endl; G4cout << "Is from primary EM shower " << fFromPrimaryEMShower << G4endl; + G4cout << "Is from charmed hadron " << fFromCharmedHadron << G4endl; } diff --git a/Pinpoint/src/TrackingAction.cc b/Pinpoint/src/TrackingAction.cc index f27147e..d04c195 100644 --- a/Pinpoint/src/TrackingAction.cc +++ b/Pinpoint/src/TrackingAction.cc @@ -107,6 +107,56 @@ void TrackingAction::PostUserTrackingAction(const G4Track* aTrack) } } + // Track charmed hadrons and their decay products (including through neutral particles) + // Charmed hadron PDG codes: 411 (D+), 421 (D0), 431 (Ds+), 4122 (Lambda_c+), 4112 (Sigma_c0), 4212 (Sigma_c+), 4222 (Sigma_c++), 4132 (Xi_c0), 4232 (Xi_c+) + G4int absPDG = std::abs(aTrack->GetParticleDefinition()->GetPDGEncoding()); + G4bool isCharmedHadron = (absPDG == 411 || absPDG == 421 || absPDG == 431 || + absPDG == 4122 || absPDG == 4112 || absPDG == 4212 || + absPDG == 4222 || absPDG == 4132 || absPDG == 4232); + + G4bool fromCharmedHadron = info && info->IsTrackFromCharmedHadron(); + + // If this is a primary charmed hadron, its decay products are counted as primary particles + if (isCharmedHadron && (aTrack->GetParentID() == 0 || aTrack->GetParentID() == 1)) + { + G4TrackVector* secondaries = fpTrackingManager->GimmeSecondaries(); + if (secondaries) + { + size_t nSeco = secondaries->size(); + if (nSeco>0) + { + for (size_t i=0; iGetCreatorProcess()->GetProcessName()=="Decay") + { + TrackInformation* newInfo = new TrackInformation(); + newInfo->SetTrackIsFromCharmedHadron(1); + (*secondaries)[i]->SetUserInformation(newInfo); + AnalysisManager::GetInstance()->AddOnePrimaryTrack(); + } + } + } + } + } + // If this track comes from a charmed hadron, tag all its decay products + else if (isCharmedHadron || fromCharmedHadron) + { + G4TrackVector* secondaries = fpTrackingManager->GimmeSecondaries(); + if (secondaries) + { + size_t nSeco = secondaries->size(); + if (nSeco>0) + { + for (size_t i=0; iSetTrackIsFromCharmedHadron(1); + (*secondaries)[i]->SetUserInformation(newInfo); + } + } + } + } + if (aTrack->GetParentID()==1 && aTrack->GetCreatorProcess()->GetProcessName()=="Decay") { // in case of tau decay pizero