Fix z vertex position and add charm information

Pinpoint/include/AnalysisManager.hh

@@ -194,6 +194,7 @@ class AnalysisManager {
     // std::vector<G4bool> fPixelFromFSLPizero;
     std::vector<G4bool> fPixelFromPrimaryLepton;
     std::vector<G4bool> fPixelFromPrimaryEMShower;
+    std::vector<G4bool> fPixelFromCharmedHadron;
 
     // Truth position of hit in x, y, z
     std::vector<Float_t> fPixelTruthX;

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;

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

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)
           // {

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);

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}
 {
 }

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)
     );
   }
 

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();

Pinpoint/src/TrackInformation.cc

@@ -10,6 +10,7 @@ TrackInformation::TrackInformation()
   fFromFSLPizero = 0;
   fFromPrimaryLepton = 0;
   fFromPrimaryEMShower = 0;
+  fFromCharmedHadron = 0;
   fHitIndices = std::make_shared<std::vector<G4long>>();
 }
 
@@ -20,6 +21,7 @@ TrackInformation::TrackInformation(const G4Track* aTrack)
   fFromFSLPizero = 0;
   fFromPrimaryLepton = 0;
   fFromPrimaryEMShower = 0;
+  fFromCharmedHadron = 0;
   fHitIndices = std::make_shared<std::vector<G4long>>();
 }
 
@@ -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;
 }
 

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; i<nSeco; ++i)
+        {
+          if ((*secondaries)[i]->GetCreatorProcess()->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; i<nSeco; ++i)
+        {
+          TrackInformation* newInfo = new TrackInformation();
+          newInfo->SetTrackIsFromCharmedHadron(1);
+          (*secondaries)[i]->SetUserInformation(newInfo);
+        }
+      }
+    }
+  }
+
   if (aTrack->GetParentID()==1 && aTrack->GetCreatorProcess()->GetProcessName()=="Decay") 
   {
     // in case of tau decay pizero

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

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);