From d2f4ecb4269cddb42d773cd49be93e412d473bf6 Mon Sep 17 00:00:00 2001 From: Tobias Boeckh Date: Thu, 29 Jan 2026 13:44:53 +0100 Subject: [PATCH 1/4] 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 22a6ddd9486162c460bfdfa5147d01a0c115062d Mon Sep 17 00:00:00 2001 From: Tobias Boeckh Date: Mon, 2 Feb 2026 00:28:48 +0100 Subject: [PATCH 2/4] Only write out tracks with 10GeV kin. energy, disable position of all hits on track, add theta --- Pinpoint/include/AnalysisManager.hh | 7 +++--- Pinpoint/src/AnalysisManager.cc | 37 ++++++++++++++++------------- 2 files changed, 24 insertions(+), 20 deletions(-) diff --git a/Pinpoint/include/AnalysisManager.hh b/Pinpoint/include/AnalysisManager.hh index 2325215..50729cc 100644 --- a/Pinpoint/include/AnalysisManager.hh +++ b/Pinpoint/include/AnalysisManager.hh @@ -130,9 +130,10 @@ class AnalysisManager { int trackPDG; double trackKinE; int trackNPoints; - std::vector trackPointX; - std::vector trackPointY; - std::vector trackPointZ; + double trackTheta; + // std::vector trackPointX; + // std::vector trackPointY; + // std::vector trackPointZ; //--------------------------------------------------- // Output variables for PRIMARIES tree diff --git a/Pinpoint/src/AnalysisManager.cc b/Pinpoint/src/AnalysisManager.cc index 525379a..21d95eb 100644 --- a/Pinpoint/src/AnalysisManager.cc +++ b/Pinpoint/src/AnalysisManager.cc @@ -146,9 +146,10 @@ void AnalysisManager::bookTrkTree() fTrk->Branch("trackPDG", &trackPDG, "trackPDG/I"); fTrk->Branch("trackKinE", &trackKinE, "trackKinE/D"); fTrk->Branch("trackNPoints", &trackNPoints, "trackNPoints/I"); - fTrk->Branch("trackPointX", &trackPointX); - fTrk->Branch("trackPointY", &trackPointY); - fTrk->Branch("trackPointZ", &trackPointZ); + fTrk->Branch("trackTheta", &trackTheta, "trackTheta/D"); + // fTrk->Branch("trackPointX", &trackPointX); + // fTrk->Branch("trackPointY", &trackPointY); + // fTrk->Branch("trackPointZ", &trackPointZ); } @@ -292,9 +293,9 @@ void AnalysisManager::BeginOfEvent() // track ID to primary ancestor association trackToPrimaryAncestor.clear(); - trackPointX.clear(); - trackPointY.clear(); - trackPointZ.clear(); + // trackPointX.clear(); + // trackPointY.clear(); + // trackPointZ.clear(); fPixelRowIDs.clear(); fPixelColIDs.clear(); @@ -497,23 +498,25 @@ void AnalysisManager::FillTrajectoriesTree(const G4Event* event) for (size_t i = 0; i < trajectoryContainer->entries(); ++i) { auto trajectory = static_cast((*trajectoryContainer)[i]); + trackKinE = trajectory->GetInitialKineticEnergy(); + if (trackKinE < 10*GeV) continue; trackTID = trajectory->GetTrackID(); trackPID = trajectory->GetParentID(); trackPDG = trajectory->GetPDGEncoding(); - trackKinE = trajectory->GetInitialKineticEnergy(); trackNPoints = trajectory->GetPointEntries(); + trackTheta = trajectory->GetInitialMomentum().theta(); count_tracks++; - for (size_t j = 0; j < trackNPoints; ++j) - { - G4ThreeVector pos = trajectory->GetPoint(j)->GetPosition(); - trackPointX.push_back( pos.x() ); - trackPointY.push_back( pos.y() ); - trackPointZ.push_back( pos.z() ); - } + // for (size_t j = 0; j < trackNPoints; ++j) + // { + // G4ThreeVector pos = trajectory->GetPoint(j)->GetPosition(); + // trackPointX.push_back( pos.x() ); + // trackPointY.push_back( pos.y() ); + // trackPointZ.push_back( pos.z() ); + // } fTrk->Fill(); - trackPointX.clear(); - trackPointY.clear(); - trackPointZ.clear(); + // trackPointX.clear(); + // trackPointY.clear(); + // trackPointZ.clear(); } G4cout << "Total number of recorded track: " << count_tracks << std::endl; } From 2033e0fe1dcea8ba0b3e5424725ffaa8960e6b54 Mon Sep 17 00:00:00 2001 From: Tobias Boeckh Date: Mon, 2 Feb 2026 11:37:58 +0100 Subject: [PATCH 3/4] Add FASER spectrometer --- Pinpoint/include/AnalysisManager.hh | 20 +++- Pinpoint/include/DetectorConstruction.hh | 15 +++ Pinpoint/include/FaserHit.hh | 83 +++++++++++++++++ Pinpoint/include/FaserSD.hh | 25 +++++ Pinpoint/src/AnalysisManager.cc | 77 ++++++++++++++++ Pinpoint/src/DetectorConstruction.cc | 112 ++++++++++++++++++++++- Pinpoint/src/FaserHit.cc | 52 +++++++++++ Pinpoint/src/FaserSD.cc | 66 +++++++++++++ 8 files changed, 447 insertions(+), 3 deletions(-) create mode 100644 Pinpoint/include/FaserHit.hh create mode 100644 Pinpoint/include/FaserSD.hh create mode 100644 Pinpoint/src/FaserHit.cc create mode 100644 Pinpoint/src/FaserSD.cc diff --git a/Pinpoint/include/AnalysisManager.hh b/Pinpoint/include/AnalysisManager.hh index 50729cc..403b9a6 100644 --- a/Pinpoint/include/AnalysisManager.hh +++ b/Pinpoint/include/AnalysisManager.hh @@ -53,6 +53,7 @@ class AnalysisManager { void bookGeomTree(); void bookHitsTrees(); void bookScintTrees(); + void bookFaserTree(); void FillEventTree(const G4Event* event); void FillPrimariesTree(const G4Event* event); @@ -60,6 +61,7 @@ class AnalysisManager { void FillGeomTree(); void FillHitsOutput(); void FillScintOutput(); + void FillFaserOutput(); float_t GetTotalEnergy(float_t px, float_t py, float_t pz, float_t m); @@ -90,6 +92,7 @@ class AnalysisManager { TTree* fPixelHitsTree; TTree* fScintTree = nullptr; + TTree* fFaserHitsTree = nullptr; // track to primary ancestor std::map trackToPrimaryAncestor; @@ -209,7 +212,22 @@ class AnalysisManager { std::vector fScintEdep; std::vector fScintFromMuon; std::vector fScintFromPrimaryLepton; - + + //---------------------------------------------------- + //OUTPUT VARIABLES FOR FASER SPECTROMETER (fFaserHitsTree) + UInt_t fFaserEventID; + std::vector fFaserTrackerID; // Tracker layer: 1, 2, or 3 + std::vector fFaserTrackID; + std::vector fFaserPDG; + std::vector fFaserX; // Hit x position + std::vector fFaserY; // Hit y position + std::vector fFaserZ; // Hit z position + std::vector fFaserPx; // Momentum x + std::vector fFaserPy; // Momentum y + std::vector fFaserPz; // Momentum z + std::vector fFaserE; // Energy + std::vector fFaserEdep; // Energy deposit + std::vector fFaserCharge; // Particle charge }; #endif diff --git a/Pinpoint/include/DetectorConstruction.hh b/Pinpoint/include/DetectorConstruction.hh index 4940b70..39d1690 100644 --- a/Pinpoint/include/DetectorConstruction.hh +++ b/Pinpoint/include/DetectorConstruction.hh @@ -151,6 +151,21 @@ class DetectorConstruction : public G4VUserDetectorConstruction G4int sim_flag = 0; G4bool scint_bar_flag = true; G4int scintContainer_CopyNumber = 0; + G4double fMagneticField = 0.57 * tesla; + G4double fLongMagnetLength = 1500.0 * mm; + G4double fShortMagnetLength = 1000.0 * mm; + G4double fInnerRadius = 100.0 * mm; + G4double fOuterRadius = 215.0 * mm; + // Position of FASER magnets and tracking stations relative to VetoNu scintillator + G4double fVetoNuPosition = -3112 * mm; + G4double fMagnet0Position = -815.3 * mm - fVetoNuPosition; + G4double fMagnet1Position = 637.4 * mm - fVetoNuPosition; + G4double fMagnet2Position = 1837.4 * mm - fVetoNuPosition; + + G4double fTrackerSize = 250.0 * mm; + G4double fTracker1Position = 47.4 * mm - fVetoNuPosition; + G4double fTracker2Position = 1237.4 * mm - fVetoNuPosition; + G4double fTracker3Position = 2427.4 * mm - fVetoNuPosition; G4bool fCheckOverlaps = true; diff --git a/Pinpoint/include/FaserHit.hh b/Pinpoint/include/FaserHit.hh new file mode 100644 index 0000000..1479bde --- /dev/null +++ b/Pinpoint/include/FaserHit.hh @@ -0,0 +1,83 @@ +#ifndef pinpoint_FaserHit_hh +#define pinpoint_FaserHit_hh + +#include "G4VHit.hh" +#include "G4THitsCollection.hh" +#include "G4Allocator.hh" +#include "G4ThreeVector.hh" +#include "G4Threading.hh" +#include "G4LorentzVector.hh" + +class FaserHit : public G4VHit +{ +public: + FaserHit() = default; + FaserHit(G4int trackerID, G4int trackID, G4int parentID, G4int pdgc, + const G4ThreeVector& pos, const G4ThreeVector& mom, + G4double energy, G4double edep, G4double charge); + FaserHit(const FaserHit&) = default; + ~FaserHit() override = default; + + FaserHit& operator=(const FaserHit&) = default; + G4bool operator==(const FaserHit&) const; + + inline void* operator new(size_t); + inline void operator delete(void*); + + void Draw() override; + void Print() override; + + void SetTrackerID(G4int id) { fTrackerID = id; } + void SetTrackID(G4int trackID) { fTrackID = trackID; } + void SetParentID(G4int parentID) { fParentID = parentID; } + void SetPDGCode(G4int pdg) { fPDGCode = pdg; } + void SetPosition(const G4ThreeVector& pos) { fPosition = pos; } + void SetMomentum(const G4ThreeVector& mom) { fMomentum = mom; } + void SetEnergy(G4double energy) { fEnergy = energy; } + void SetEnergyDeposit(G4double edep) { fEnergyDeposit = edep; } + void SetCharge(G4double charge) { fCharge = charge; } + + G4int GetTrackerID() const { return fTrackerID; } + G4int GetTrackID() const { return fTrackID; } + G4int GetParentID() const { return fParentID; } + G4int GetPDGCode() const { return fPDGCode; } + G4ThreeVector GetPosition() const { return fPosition; } + G4ThreeVector GetMomentum() const { return fMomentum; } + G4double GetEnergy() const { return fEnergy; } + G4double GetEnergyDeposit() const { return fEnergyDeposit; } + G4double GetCharge() const { return fCharge; } + G4double GetPx() const { return fMomentum.x(); } + G4double GetPy() const { return fMomentum.y(); } + G4double GetPz() const { return fMomentum.z(); } + G4double GetX() const { return fPosition.x(); } + G4double GetY() const { return fPosition.y(); } + G4double GetZ() const { return fPosition.z(); } + +private: + G4int fTrackerID = -1; // FASER tracker layer: 0, 1, or 2 + G4int fTrackID = -1; + G4int fParentID = -1; + G4int fPDGCode = 0; + G4ThreeVector fPosition; + G4ThreeVector fMomentum; + G4double fEnergy = 0.0; + G4double fEnergyDeposit = 0.0; + G4double fCharge = 0.0; +}; + +using FaserHitsCollection = G4THitsCollection; + +extern G4ThreadLocal G4Allocator* FaserHitAllocator; + +inline void* FaserHit::operator new(size_t) +{ + if (!FaserHitAllocator) FaserHitAllocator = new G4Allocator; + return (void*)FaserHitAllocator->MallocSingle(); +} + +inline void FaserHit::operator delete(void* hit) +{ + FaserHitAllocator->FreeSingle((FaserHit*)hit); +} + +#endif diff --git a/Pinpoint/include/FaserSD.hh b/Pinpoint/include/FaserSD.hh new file mode 100644 index 0000000..1575efe --- /dev/null +++ b/Pinpoint/include/FaserSD.hh @@ -0,0 +1,25 @@ +#ifndef pinpoint_FaserSD_hh +#define pinpoint_FaserSD_hh + +#include "FaserHit.hh" +#include "G4VSensitiveDetector.hh" +#include "G4SystemOfUnits.hh" + +class G4Step; +class G4HCofThisEvent; + +class FaserSD : public G4VSensitiveDetector +{ +public: + FaserSD(const G4String& name, const G4String& hitsCollectionName); + ~FaserSD() override = default; + + void Initialize(G4HCofThisEvent* hitCollection) override; + G4bool ProcessHits(G4Step* step, G4TouchableHistory* history) override; + void EndOfEvent(G4HCofThisEvent* hitCollection) override; + +private: + FaserHitsCollection* fHitsCollection = nullptr; +}; + +#endif diff --git a/Pinpoint/src/AnalysisManager.cc b/Pinpoint/src/AnalysisManager.cc index 21d95eb..49a499a 100644 --- a/Pinpoint/src/AnalysisManager.cc +++ b/Pinpoint/src/AnalysisManager.cc @@ -34,6 +34,7 @@ #include "FPFParticle.hh" #include "PixelHit.hh" #include "ScintHit.hh" +#include "FaserHit.hh" //--------------------------------------------------------------------- @@ -230,6 +231,27 @@ void AnalysisManager::bookScintTrees() fScintTree->Branch("fromPrimaryLepton", &fScintFromPrimaryLepton); } + +void AnalysisManager::bookFaserTree() +{ + fFile->cd(fHits->GetName()); + + fFaserHitsTree = new TTree("faserHits", "FASER spectrometer hits"); + fFaserHitsTree->Branch("event_id", &fFaserEventID, "event_id/i"); + fFaserHitsTree->Branch("trackerID", &fFaserTrackerID); + fFaserHitsTree->Branch("trackID", &fFaserTrackID); + fFaserHitsTree->Branch("pdg", &fFaserPDG); + fFaserHitsTree->Branch("x", &fFaserX); + fFaserHitsTree->Branch("y", &fFaserY); + fFaserHitsTree->Branch("z", &fFaserZ); + fFaserHitsTree->Branch("px", &fFaserPx); + fFaserHitsTree->Branch("py", &fFaserPy); + fFaserHitsTree->Branch("pz", &fFaserPz); + fFaserHitsTree->Branch("energy", &fFaserE); + fFaserHitsTree->Branch("edep", &fFaserEdep); + fFaserHitsTree->Branch("charge", &fFaserCharge); +} + //--------------------------------------------------------------------- //--------------------------------------------------------------------- @@ -251,6 +273,7 @@ void AnalysisManager::BeginOfRun() bookHitsTrees(); bookScintTrees(); + bookFaserTree(); } @@ -271,6 +294,7 @@ void AnalysisManager::EndOfRun() fFile->cd(fHits->GetName()); fPixelHitsTree->Write(); fScintTree->Write(); + fFaserHitsTree->Write(); // fActsParticlesTree->Write(); fFile->cd(); // go back to top @@ -361,6 +385,7 @@ void AnalysisManager::EndOfEvent(const G4Event *event) FillHitsOutput(); FillScintOutput(); + FillFaserOutput(); } //--------------------------------------------------------------------- @@ -644,6 +669,58 @@ void AnalysisManager::FillScintOutput() fScintTree->Fill(); } +//// FASER SPECTROMETER HITS --- +void AnalysisManager::FillFaserOutput() +{ + fFaserEventID = evtID; + + // Clear vectors + fFaserTrackerID.clear(); + fFaserTrackID.clear(); + fFaserPDG.clear(); + fFaserX.clear(); + fFaserY.clear(); + fFaserZ.clear(); + fFaserPx.clear(); + fFaserPy.clear(); + fFaserPz.clear(); + fFaserE.clear(); + fFaserEdep.clear(); + fFaserCharge.clear(); + + G4int nHC = fHCofEvent->GetNumberOfCollections(); + + for(G4int i = 0; i < nHC; ++i) + { + auto* hc = fHCofEvent->GetHC(i); + auto* faserHC = dynamic_cast(hc); + if(!faserHC) continue; + if(faserHC->GetName() != "FaserHitsCollection") continue; + + G4cout << "Found FASER tracker hit collection with " << faserHC->GetSize() << " hits" << G4endl; + + for(size_t h = 0; h < faserHC->entries(); ++h) + { + auto* hit = (*faserHC)[h]; + + fFaserTrackerID.push_back(hit->GetTrackerID()); + fFaserTrackID.push_back(hit->GetTrackID()); + fFaserPDG.push_back(hit->GetPDGCode()); + fFaserX.push_back(hit->GetX()); + fFaserY.push_back(hit->GetY()); + fFaserZ.push_back(hit->GetZ()); + fFaserPx.push_back(hit->GetPx()); + fFaserPy.push_back(hit->GetPy()); + fFaserPz.push_back(hit->GetPz()); + fFaserE.push_back(hit->GetEnergy()); + fFaserEdep.push_back(hit->GetEnergyDeposit()); + fFaserCharge.push_back(hit->GetCharge()); + } + } + + fFaserHitsTree->Fill(); +} + float_t AnalysisManager::GetTotalEnergy(float_t px, float_t py, float_t pz, float_t m) { return TMath::Sqrt(px * px + py * py + pz * pz + m * m); diff --git a/Pinpoint/src/DetectorConstruction.cc b/Pinpoint/src/DetectorConstruction.cc index fe2f4be..6c12d17 100644 --- a/Pinpoint/src/DetectorConstruction.cc +++ b/Pinpoint/src/DetectorConstruction.cc @@ -3,17 +3,26 @@ #include "DetectorConstruction.hh" #include "PixelSD.hh" #include "ScintSD.hh" +#include "FaserSD.hh" #include "G4LogicalVolume.hh" #include "G4PVPlacement.hh" #include "G4Box.hh" +#include "G4Tubs.hh" #include "G4Cons.hh" #include "G4PVReplica.hh" #include "G4PVParameterised.hh" #include "G4SDManager.hh" +#include "G4UniformMagField.hh" +#include "G4FieldManager.hh" +#include "G4TransportationManager.hh" +#include "G4Mag_UsualEqRhs.hh" +#include "G4ClassicalRK4.hh" +#include "G4ChordFinder.hh" #include #include "G4VisAttributes.hh" #include "G4LogicalSkinSurface.hh" #include "G4PhysicalVolumeStore.hh" +#include "G4LogicalVolumeStore.hh" #include #include @@ -34,7 +43,17 @@ void DetectorConstruction::DefineMaterial() //Scintillator Material and Properties G4NistManager* nist = G4NistManager::Instance(); scintillator = nist->FindOrBuildMaterial("G4_PLASTIC_SC_VINYLTOLUENE"); - scintillator->GetIonisation()->SetBirksConstant(0.126 * mm / MeV); + scintillator->GetIonisation()->SetBirksConstant(0.126 * mm / MeV); + + // Define Samarium Cobalt (Sm2Co17) material + // Density: 8.4 g/cm³ + G4double density = 8.4 * g/cm3; + G4int nComponents = 2; + G4Material* sm2co17 = new G4Material("Sm2Co17", density, nComponents); + G4Element* Sm = nist->FindOrBuildElement("Sm"); + G4Element* Co = nist->FindOrBuildElement("Co"); + sm2co17->AddElement(Sm, 2); + sm2co17->AddElement(Co, 17); // std::vector refractiveIndexScint = { 1.58, 1.58 }; // std::vector absorptionScint = {0.1*cm, 0.1*cm}; @@ -119,7 +138,7 @@ G4VPhysicalVolume* DetectorConstruction::Construct() auto detectorThickness = fNLayers * fLayerThickness; auto worldSizeX = 1.2 * fDetectorWidth; auto worldSizeY = 1.2 * fDetectorHeight; - auto worldSizeZ = 1.2 * detectorThickness; + auto worldSizeZ = 1.2 * (detectorThickness + fTracker3Position); // Get materials G4NistManager* nist = G4NistManager::Instance(); @@ -307,6 +326,60 @@ G4VPhysicalVolume* DetectorConstruction::Construct() // G4cout << "Detector consists of " << fNLayers << " layers of: [ " << fTungstenThickness / mm << "mm of " << tungstenMaterial->GetName() << " + " << fSiliconThickness / mm << "mm of " // << siliconMaterial->GetName() << " + " << fBoxThickness / mm << "mm of " << worldMaterial->GetName() << " ] " << G4endl; + + // FASER spectrometer magnets + G4Material* sm2co17 = G4Material::GetMaterial("Sm2Co17"); + + auto longMangetS = 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); + + // Magnet 0 + auto magnet0LV = new G4LogicalVolume(longMangetS, 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 + + // 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 + + // 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 + + G4cout << "Created cylindrical magnet regions:" << G4endl; + G4cout << " Inner radius: " << fInnerRadius/mm << " mm, Outer radius: " << fOuterRadius/mm << " mm" << G4endl; + G4cout << " Magnet 0: z = " << fMagnet0Position/mm << " mm, length = " << fLongMagnetLength/mm << " mm" << G4endl; + G4cout << " Magnet 1: z = " << fMagnet1Position/mm << " mm, length = " << fShortMagnetLength/mm << " mm" << G4endl; + G4cout << " Magnet 2: z = " << fMagnet2Position/mm << " mm, length = " << fShortMagnetLength/mm << " mm" << G4endl; + G4cout << " Material: Sm2Co17" << G4endl; + + // FASER spectrometer tracking layers (use single layer per station) + auto trackerS = new G4Box("Tracker", 0.5 * fTrackerSize, 0.5 * fTrackerSize, 0.5 * fSiliconThickness); + + // Tracker 1 + auto tracker1LV = new G4LogicalVolume(trackerS, siliconMaterial, "Tracker1"); + new G4PVPlacement(nullptr, G4ThreeVector(0., 0., fTracker1Position), tracker1LV, "Tracker1", worldLV, false, 0, fCheckOverlaps); + tracker1LV->SetVisAttributes(G4VisAttributes(G4Colour(0.0, 1.0, 0.0, 0.7))); // Green, semi-transparent + + // Tracker 2 + auto tracker2LV = new G4LogicalVolume(trackerS, siliconMaterial, "Tracker2"); + new G4PVPlacement(nullptr, G4ThreeVector(0., 0., fTracker2Position), tracker2LV, "Tracker2", worldLV, false, 1, fCheckOverlaps); + tracker2LV->SetVisAttributes(G4VisAttributes(G4Colour(0.0, 1.0, 0.0, 0.7))); // Green, semi-transparent + + // Tracker 3 + auto tracker3LV = new G4LogicalVolume(trackerS, siliconMaterial, "Tracker3"); + new G4PVPlacement(nullptr, G4ThreeVector(0., 0., fTracker3Position), tracker3LV, "Tracker3", worldLV, false, 2, fCheckOverlaps); + tracker3LV->SetVisAttributes(G4VisAttributes(G4Colour(0.0, 1.0, 0.0, 0.7))); // Green, semi-transparent + + G4cout << "Created tracking layers:" << G4endl; + G4cout << " Size: " << fTrackerSize/mm << " x " << fTrackerSize/mm << " mm, Thickness: " << fSiliconThickness/um << " um" << G4endl; + G4cout << " Tracker 1: z = " << fTracker1Position/mm << " mm" << G4endl; + G4cout << " Tracker 2: z = " << fTracker2Position/mm << " mm" << G4endl; + G4cout << " Tracker 3: z = " << fTracker3Position/mm << " mm" << G4endl; + G4cout << " Material: " << siliconMaterial->GetName() << G4endl; + // Write GDML file if it doesn't exist std::ifstream file(fWriteFile); if (!file.good()) { @@ -354,6 +427,41 @@ void DetectorConstruction::ConstructSDandField() ComputeSiliconZPositions(); ComputePixelCentersXY(); } + + // Tracker SD + G4cout << "Adding tracker SD" << G4endl; + auto faserSD = new FaserSD("FaserSpectrometer", "FaserHitsCollection"); + G4SDManager::GetSDMpointer()->AddNewDetector(faserSD); + + G4LogicalVolume* tracker1LV = G4LogicalVolumeStore::GetInstance()->GetVolume("Tracker1"); + G4LogicalVolume* tracker2LV = G4LogicalVolumeStore::GetInstance()->GetVolume("Tracker2"); + G4LogicalVolume* tracker3LV = G4LogicalVolumeStore::GetInstance()->GetVolume("Tracker3"); + + if(tracker1LV) tracker1LV->SetSensitiveDetector(faserSD); + if(tracker2LV) tracker2LV->SetSensitiveDetector(faserSD); + if(tracker3LV) tracker3LV->SetSensitiveDetector(faserSD); + + // Create uniform magnetic field (pointing in X direction) + G4ThreeVector fieldVector(fMagneticField, 0., 0.); + auto magField = new G4UniformMagField(fieldVector); + + // Create field manager and chord finder + auto fieldMgr = new G4FieldManager(); + fieldMgr->SetDetectorField(magField); + fieldMgr->CreateChordFinder(magField); + + // Get the magnet logical volumes and assign field manager + G4LogicalVolume* magnet0LV = G4LogicalVolumeStore::GetInstance()->GetVolume("Magnet0"); + G4LogicalVolume* magnet1LV = G4LogicalVolumeStore::GetInstance()->GetVolume("Magnet1"); + G4LogicalVolume* magnet2LV = G4LogicalVolumeStore::GetInstance()->GetVolume("Magnet2"); + + if(magnet0LV) magnet0LV->SetFieldManager(fieldMgr, true); + if(magnet1LV) magnet1LV->SetFieldManager(fieldMgr, true); + if(magnet2LV) magnet2LV->SetFieldManager(fieldMgr, true); + + G4cout << "Configured magnetic field:" << G4endl; + G4cout << " Field strength: " << fMagneticField/tesla << " T (X-direction)" << G4endl; + G4cout << " Applied to Magnet0, Magnet1, and Magnet2" << G4endl; } diff --git a/Pinpoint/src/FaserHit.cc b/Pinpoint/src/FaserHit.cc new file mode 100644 index 0000000..560b3d4 --- /dev/null +++ b/Pinpoint/src/FaserHit.cc @@ -0,0 +1,52 @@ +#include "FaserHit.hh" +#include "G4UnitsTable.hh" +#include "G4VVisManager.hh" +#include "G4Circle.hh" +#include "G4Colour.hh" +#include "G4VisAttributes.hh" +#include "G4SystemOfUnits.hh" + +#include + +G4ThreadLocal G4Allocator* FaserHitAllocator = nullptr; + +G4bool FaserHit::operator==(const FaserHit& right) const +{ + return ( this == &right ) ? true : false; +} + +void FaserHit::Draw() +{ + G4VVisManager* pVVisManager = G4VVisManager::GetConcreteInstance(); + if(pVVisManager) + { + G4Circle circle(fPosition); + circle.SetScreenSize(4.); + circle.SetFillStyle(G4Circle::filled); + G4Colour colour(0.,1.,0.); + G4VisAttributes attribs(colour); + circle.SetVisAttributes(attribs); + pVVisManager->Draw(circle); + } +} + +void FaserHit::Print() +{ + G4cout + << " TrackerID: " << fTrackerID + << " TrackID: " << fTrackID + << " PDG: " << fPDGCode + << " Position: (" << fPosition.x()/mm << ", " << fPosition.y()/mm << ", " << fPosition.z()/mm << ") mm" + << " Momentum: (" << fMomentum.x()/GeV << ", " << fMomentum.y()/GeV << ", " << fMomentum.z()/GeV << ") GeV" + << " Energy: " << fEnergy/GeV << " GeV" + << " Edep: " << G4BestUnit(fEnergyDeposit,"Energy") + << G4endl; +} + +FaserHit::FaserHit(G4int trackerID, G4int trackID, G4int parentID, G4int pdgc, + const G4ThreeVector& pos, const G4ThreeVector& mom, + G4double energy, G4double edep, G4double charge) +: fTrackerID{trackerID}, fTrackID{trackID}, fParentID{parentID}, fPDGCode{pdgc}, + fPosition{pos}, fMomentum{mom}, fEnergy{energy}, fEnergyDeposit{edep}, fCharge{charge} +{ +} diff --git a/Pinpoint/src/FaserSD.cc b/Pinpoint/src/FaserSD.cc new file mode 100644 index 0000000..c5877e3 --- /dev/null +++ b/Pinpoint/src/FaserSD.cc @@ -0,0 +1,66 @@ +#include "FaserSD.hh" +#include "G4HCofThisEvent.hh" +#include "G4Step.hh" +#include "G4ThreeVector.hh" +#include "G4SDManager.hh" +#include "G4Track.hh" +#include "G4VTouchable.hh" +#include "G4TouchableHistory.hh" +#include "G4ParticleDefinition.hh" +#include "G4ios.hh" + +FaserSD::FaserSD(const G4String& name, const G4String& hitsCollectionName) + : G4VSensitiveDetector(name) +{ + collectionName.insert(hitsCollectionName); +} + +void FaserSD::Initialize(G4HCofThisEvent* hce) +{ + fHitsCollection = new FaserHitsCollection(SensitiveDetectorName, collectionName[0]); + + G4int hcID = G4SDManager::GetSDMpointer()->GetCollectionID(collectionName[0]); + hce->AddHitsCollection(hcID, fHitsCollection); +} + +G4bool FaserSD::ProcessHits(G4Step* step, G4TouchableHistory* /*history*/) +{ + const auto* track = step->GetTrack(); + const auto* preStepPoint = step->GetPreStepPoint(); + const auto& touchable = preStepPoint->GetTouchableHandle(); + + G4double edep = step->GetTotalEnergyDeposit(); + + G4double charge = track->GetDefinition()->GetPDGCharge(); + if (charge == 0) { + return false; + } + + // Get tracker ID from copy number (0, 1, or 2) + G4int trackerID = touchable->GetCopyNumber(); + + G4int trackID = track->GetTrackID(); + G4int parentID = track->GetParentID(); + G4int pdgCode = track->GetParticleDefinition()->GetPDGEncoding(); + + G4ThreeVector position = preStepPoint->GetPosition(); + G4ThreeVector momentum = preStepPoint->GetMomentum(); + G4double energy = preStepPoint->GetTotalEnergy(); + + auto* hit = new FaserHit(trackerID, trackID, parentID, pdgCode, + position, momentum, energy, edep, charge); + fHitsCollection->insert(hit); + + return true; +} + +void FaserSD::EndOfEvent(G4HCofThisEvent* /*hce*/) +{ + if (verboseLevel > 1) { + std::size_t nofHits = fHitsCollection->entries(); + G4cout << G4endl << "-------->Hits Collection: in this event there are " << nofHits + << " hits in the Faser tracker detector: " << G4endl; + for (std::size_t i = 0; i < nofHits; i++) + (*fHitsCollection)[i]->Print(); + } +} From fd4e450dd8f08762aac3f07dabc1d2944321da4b Mon Sep 17 00:00:00 2001 From: Tobias Boeckh Date: Mon, 2 Feb 2026 16:04:07 +0100 Subject: [PATCH 4/4] Fix magnetic field (make magnets solid with air-filled bore --- Pinpoint/include/DetectorConstruction.hh | 6 +- Pinpoint/include/MagneticField.hh | 50 ++++++++++++++ Pinpoint/src/DetectorConstruction.cc | 83 ++++++++++++++++-------- 3 files changed, 112 insertions(+), 27 deletions(-) create mode 100644 Pinpoint/include/MagneticField.hh diff --git a/Pinpoint/include/DetectorConstruction.hh b/Pinpoint/include/DetectorConstruction.hh index 39d1690..5ba1ce1 100644 --- a/Pinpoint/include/DetectorConstruction.hh +++ b/Pinpoint/include/DetectorConstruction.hh @@ -9,6 +9,8 @@ #include "G4OpticalSurface.hh" class G4VPhysicalVolume; +class G4FieldManager; +class MagneticField; class DetectorConstruction : public G4VUserDetectorConstruction { @@ -151,7 +153,6 @@ class DetectorConstruction : public G4VUserDetectorConstruction G4int sim_flag = 0; G4bool scint_bar_flag = true; G4int scintContainer_CopyNumber = 0; - G4double fMagneticField = 0.57 * tesla; G4double fLongMagnetLength = 1500.0 * mm; G4double fShortMagnetLength = 1000.0 * mm; G4double fInnerRadius = 100.0 * mm; @@ -179,6 +180,9 @@ class DetectorConstruction : public G4VUserDetectorConstruction G4OpticalSurface* scintWrap; G4Material* scintillator = nullptr; + + static G4ThreadLocal MagneticField* fMagneticField; + static G4ThreadLocal G4FieldManager* fFieldMgr; }; diff --git a/Pinpoint/include/MagneticField.hh b/Pinpoint/include/MagneticField.hh new file mode 100644 index 0000000..524f473 --- /dev/null +++ b/Pinpoint/include/MagneticField.hh @@ -0,0 +1,50 @@ +#ifndef MagneticField_h +#define MagneticField_h 1 + +#include "G4MagneticField.hh" +#include "MagneticField.hh" +#include "G4GenericMessenger.hh" +#include "globals.hh" + +#include + +class G4GenericMessenger; + +/// Uniform magnetic field in x direction with 0.57 T strength + +class MagneticField : public G4MagneticField +{ + public: + MagneticField() { DefineCommands(); } + ~MagneticField() override { delete fMessenger; } + + void GetFieldValue(const G4double point[4], double* bField) const override { + bField[0] = fBx; + bField[1] = 0.; + bField[2] = 0.; + } + + void SetField(G4double val) { + fBx = val; + G4cout << "MagneticField::SetField called - new field strength: " << fBx/CLHEP::tesla << " T" << G4endl; + } + G4double GetField() const { return fBx; } + + private: + void DefineCommands() + { + // Define /Pinpoint/field command directory using generic messenger class + fMessenger = new G4GenericMessenger(this, "/Pinpoint/field/", "Field control"); + + // fieldValue command + auto& valueCmd = fMessenger->DeclareMethodWithUnit("value", "tesla", &MagneticField::SetField, + "Set field strength."); + valueCmd.SetParameterName("field", true); + valueCmd.SetDefaultValue("0.57"); + } + + G4GenericMessenger* fMessenger = nullptr; + G4double fBx = 0.57 * CLHEP::tesla; +}; + +#endif diff --git a/Pinpoint/src/DetectorConstruction.cc b/Pinpoint/src/DetectorConstruction.cc index 6c12d17..ef17363 100644 --- a/Pinpoint/src/DetectorConstruction.cc +++ b/Pinpoint/src/DetectorConstruction.cc @@ -4,6 +4,7 @@ #include "PixelSD.hh" #include "ScintSD.hh" #include "FaserSD.hh" +#include "MagneticField.hh" #include "G4LogicalVolume.hh" #include "G4PVPlacement.hh" #include "G4Box.hh" @@ -12,9 +13,8 @@ #include "G4PVReplica.hh" #include "G4PVParameterised.hh" #include "G4SDManager.hh" -#include "G4UniformMagField.hh" #include "G4FieldManager.hh" -#include "G4TransportationManager.hh" +#include "G4UserLimits.hh" #include "G4Mag_UsualEqRhs.hh" #include "G4ClassicalRK4.hh" #include "G4ChordFinder.hh" @@ -26,11 +26,15 @@ #include #include +G4ThreadLocal MagneticField* DetectorConstruction::fMagneticField = nullptr; +G4ThreadLocal G4FieldManager* DetectorConstruction::fFieldMgr = nullptr; DetectorConstruction::DetectorConstruction() : G4VUserDetectorConstruction() { messenger = new DetectorConstructionMessenger(this); + + fMagneticField = new MagneticField(); } DetectorConstruction::~DetectorConstruction() @@ -327,34 +331,54 @@ G4VPhysicalVolume* DetectorConstruction::Construct() // << siliconMaterial->GetName() << " + " << fBoxThickness / mm << "mm of " << worldMaterial->GetName() << " ] " << G4endl; - // FASER spectrometer magnets + // FASER spectrometer magnets: + // solid cylinders (0 to outerRadius) and air-filled bore (0 to innerRadius) G4Material* sm2co17 = G4Material::GetMaterial("Sm2Co17"); - auto longMangetS = 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); + 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); + + 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); // Magnet 0 - auto magnet0LV = new G4LogicalVolume(longMangetS, sm2co17, "Magnet0"); + 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); + 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 fieldRegion1LV = new G4LogicalVolume(shortFieldS, worldMaterial, "FieldRegion1"); + new G4PVPlacement(nullptr, G4ThreeVector(0., 0., 0.), fieldRegion1LV, "FieldRegion1", magnet1LV, 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 fieldRegion2LV = new G4LogicalVolume(shortFieldS, worldMaterial, "FieldRegion2"); + new G4PVPlacement(nullptr, G4ThreeVector(0., 0., 0.), fieldRegion2LV, "FieldRegion2", magnet2LV, false, 2, fCheckOverlaps); + fieldRegion2LV->SetVisAttributes(G4VisAttributes::GetInvisible()); - G4cout << "Created cylindrical magnet regions:" << G4endl; - G4cout << " Inner radius: " << fInnerRadius/mm << " mm, Outer radius: " << fOuterRadius/mm << " mm" << G4endl; + G4cout << "Created solid cylindrical magnets:" << G4endl; + G4cout << " Outer radius: " << fOuterRadius/mm << " mm (solid Sm2Co17)" << G4endl; G4cout << " Magnet 0: z = " << fMagnet0Position/mm << " mm, length = " << fLongMagnetLength/mm << " mm" << G4endl; G4cout << " Magnet 1: z = " << fMagnet1Position/mm << " mm, length = " << fShortMagnetLength/mm << " mm" << G4endl; 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); + // FASER spectrometer tracking layers (use single layer per station) auto trackerS = new G4Box("Tracker", 0.5 * fTrackerSize, 0.5 * fTrackerSize, 0.5 * fSiliconThickness); @@ -428,7 +452,7 @@ void DetectorConstruction::ConstructSDandField() ComputePixelCentersXY(); } - // Tracker SD + // FASER Tracking spectrometer SD G4cout << "Adding tracker SD" << G4endl; auto faserSD = new FaserSD("FaserSpectrometer", "FaserHitsCollection"); G4SDManager::GetSDMpointer()->AddNewDetector(faserSD); @@ -441,27 +465,34 @@ void DetectorConstruction::ConstructSDandField() if(tracker2LV) tracker2LV->SetSensitiveDetector(faserSD); if(tracker3LV) tracker3LV->SetSensitiveDetector(faserSD); - // Create uniform magnetic field (pointing in X direction) - G4ThreeVector fieldVector(fMagneticField, 0., 0.); - auto magField = new G4UniformMagField(fieldVector); + // Setup magnetic field manager + fFieldMgr = new G4FieldManager(); + fFieldMgr->SetDetectorField(fMagneticField); + fFieldMgr->CreateChordFinder(fMagneticField); - // Create field manager and chord finder - auto fieldMgr = new G4FieldManager(); - fieldMgr->SetDetectorField(magField); - fieldMgr->CreateChordFinder(magField); + // Get the air-filled field region logical volumes and assign field manager + G4LogicalVolume* fieldRegion0LV = G4LogicalVolumeStore::GetInstance()->GetVolume("FieldRegion0"); + G4LogicalVolume* fieldRegion1LV = G4LogicalVolumeStore::GetInstance()->GetVolume("FieldRegion1"); + G4LogicalVolume* fieldRegion2LV = G4LogicalVolumeStore::GetInstance()->GetVolume("FieldRegion2"); - // Get the magnet logical volumes and assign field manager - G4LogicalVolume* magnet0LV = G4LogicalVolumeStore::GetInstance()->GetVolume("Magnet0"); - G4LogicalVolume* magnet1LV = G4LogicalVolumeStore::GetInstance()->GetVolume("Magnet1"); - G4LogicalVolume* magnet2LV = G4LogicalVolumeStore::GetInstance()->GetVolume("Magnet2"); + G4bool forceToAllDaughters = true; + if(fieldRegion0LV) { + fieldRegion0LV->SetFieldManager(fFieldMgr, forceToAllDaughters); + G4cout << "Assigned magnetic field to FieldRegion0" << G4endl; + } + if(fieldRegion1LV) { + fieldRegion1LV->SetFieldManager(fFieldMgr, forceToAllDaughters); + G4cout << "Assigned magnetic field to FieldRegion1" << G4endl; + } + if(fieldRegion2LV) { + fieldRegion2LV->SetFieldManager(fFieldMgr, forceToAllDaughters); + G4cout << "Assigned magnetic field to FieldRegion2" << G4endl; + } - if(magnet0LV) magnet0LV->SetFieldManager(fieldMgr, true); - if(magnet1LV) magnet1LV->SetFieldManager(fieldMgr, true); - if(magnet2LV) magnet2LV->SetFieldManager(fieldMgr, true); + G4cout << "Configured magnetic field using custom MagneticField class:" << G4endl; + G4cout << " Field strength: " << fMagneticField->GetField()/tesla << " T (X-direction)" << G4endl; + G4cout << " Applied to air-filled field regions inside magnets (r < " << fInnerRadius/mm << " mm)" << G4endl; - G4cout << "Configured magnetic field:" << G4endl; - G4cout << " Field strength: " << fMagneticField/tesla << " T (X-direction)" << G4endl; - G4cout << " Applied to Magnet0, Magnet1, and Magnet2" << G4endl; }