Faser spectrometer

Pinpoint/include/AnalysisManager.hh

@@ -53,13 +53,15 @@ class AnalysisManager {
     void bookGeomTree();
     void bookHitsTrees();
     void bookScintTrees();
+    void bookFaserTree();
 
     void FillEventTree(const G4Event* event);
     void FillPrimariesTree(const G4Event* event);
     void FillTrajectoriesTree(const G4Event* event);
     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<G4int, G4int> trackToPrimaryAncestor;
@@ -130,9 +133,10 @@ class AnalysisManager {
     int trackPDG;
     double trackKinE;
     int trackNPoints;
-    std::vector<double> trackPointX;
-    std::vector<double> trackPointY;
-    std::vector<double> trackPointZ;
+    double trackTheta;
+    // std::vector<double> trackPointX;
+    // std::vector<double> trackPointY;
+    // std::vector<double> trackPointZ;
 
     //---------------------------------------------------
     // Output variables for PRIMARIES tree
@@ -189,6 +193,7 @@ class AnalysisManager {
     // std::vector<G4bool> fPixelFromPrimaryPizero;
     // std::vector<G4bool> fPixelFromFSLPizero;
     std::vector<G4bool> fPixelFromPrimaryLepton;
+    std::vector<G4bool> fPixelFromPrimaryEMShower;
 
     // Truth position of hit in x, y, z
     std::vector<Float_t> fPixelTruthX;
@@ -207,7 +212,22 @@ class AnalysisManager {
     std::vector<float> fScintEdep;
     std::vector<int> fScintFromMuon;
     std::vector<int> fScintFromPrimaryLepton;
-
+
+    //----------------------------------------------------
+    //OUTPUT VARIABLES FOR FASER SPECTROMETER (fFaserHitsTree)
+    UInt_t fFaserEventID;
+    std::vector<int> fFaserTrackerID;     // Tracker layer: 1, 2, or 3
+    std::vector<int> fFaserTrackID;
+    std::vector<int> fFaserPDG;
+    std::vector<float> fFaserX;           // Hit x position
+    std::vector<float> fFaserY;           // Hit y position
+    std::vector<float> fFaserZ;           // Hit z position
+    std::vector<float> fFaserPx;          // Momentum x
+    std::vector<float> fFaserPy;          // Momentum y
+    std::vector<float> fFaserPz;          // Momentum z
+    std::vector<float> fFaserE;           // Energy
+    std::vector<float> fFaserEdep;        // Energy deposit
+    std::vector<float> fFaserCharge;      // Particle charge
   };
 
 #endif

Pinpoint/include/DetectorConstruction.hh

@@ -9,6 +9,8 @@
 #include "G4OpticalSurface.hh"
 
 class G4VPhysicalVolume;
+class G4FieldManager;
+class MagneticField;
 
 class DetectorConstruction : public G4VUserDetectorConstruction
 {
@@ -151,6 +153,20 @@ class DetectorConstruction : public G4VUserDetectorConstruction
     G4int sim_flag = 0;
     G4bool scint_bar_flag = true;
     G4int scintContainer_CopyNumber = 0;
+    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;
 
@@ -164,6 +180,9 @@ class DetectorConstruction : public G4VUserDetectorConstruction
     G4OpticalSurface* scintWrap;
 
     G4Material* scintillator = nullptr;
+
+    static G4ThreadLocal MagneticField* fMagneticField;
+    static G4ThreadLocal G4FieldManager* fFieldMgr;
 };
 
 

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<FaserHit>;
+
+extern G4ThreadLocal G4Allocator<FaserHit>* FaserHitAllocator;
+
+inline void* FaserHit::operator new(size_t)
+{
+  if (!FaserHitAllocator) FaserHitAllocator = new G4Allocator<FaserHit>;
+  return (void*)FaserHitAllocator->MallocSingle();
+}
+
+inline void FaserHit::operator delete(void* hit)
+{
+  FaserHitAllocator->FreeSingle((FaserHit*)hit);
+}
+
+#endif

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

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 <CLHEP/Units/SystemOfUnits.h>
+
+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

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;

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<std::vector<G4long>> fHitIndices;
 
     // std::vector<G4long>* fHitIndices; //TODO: Make this a smart pointer?