diff --git a/Pinpoint/include/AnalysisManager.hh b/Pinpoint/include/AnalysisManager.hh index 403b9a6..180553a 100644 --- a/Pinpoint/include/AnalysisManager.hh +++ b/Pinpoint/include/AnalysisManager.hh @@ -14,6 +14,8 @@ #include "AnalysisManagerMessenger.hh" #include "FPFParticle.hh" +class G4Track; + class AnalysisManager { public: @@ -54,6 +56,8 @@ class AnalysisManager { void bookHitsTrees(); void bookScintTrees(); void bookFaserTree(); + void bookTauTree(); + void bookCharmTree(); void FillEventTree(const G4Event* event); void FillPrimariesTree(const G4Event* event); @@ -62,6 +66,7 @@ class AnalysisManager { void FillHitsOutput(); void FillScintOutput(); void FillFaserOutput(); + void FillPythiaTree(const G4Event* event, TTree* tree, std::vector pdg_ids); float_t GetTotalEnergy(float_t px, float_t py, float_t pz, float_t m); @@ -87,6 +92,8 @@ class AnalysisManager { TTree* fTrk; TTree* fPrim; TTree* fGeom; + TTree* fTauTree; + TTree* fCharmTree; TDirectory* fHits; TTree* fPixelHitsTree; @@ -134,6 +141,12 @@ class AnalysisManager { double trackKinE; int trackNPoints; double trackTheta; + double trackProdX; + double trackProdY; + double trackProdZ; + double trackDecayX; + double trackDecayY; + double trackDecayZ; // std::vector trackPointX; // std::vector trackPointY; // std::vector trackPointZ; @@ -228,6 +241,24 @@ class AnalysisManager { std::vector fFaserE; // Energy std::vector fFaserEdep; // Energy deposit std::vector fFaserCharge; // Particle charge + + // + //---------------------------------------------------- + // Output variables for tau and charm trees + UInt_t fPythiaEventID; + int fPythiaTPDG; + int fPythiaPPDG; + int fPythiaTID; + int fPythiaPID; + double fPythiaProdX; + double fPythiaProdY; + double fPythiaProdZ; + double fPythiaDecayX; + double fPythiaDecayY; + double fPythiaDecayZ; + double fPythiaPx; + double fPythiaPy; + double fPythiaPz; }; #endif diff --git a/Pinpoint/src/AnalysisManager.cc b/Pinpoint/src/AnalysisManager.cc index 49a499a..14b623b 100644 --- a/Pinpoint/src/AnalysisManager.cc +++ b/Pinpoint/src/AnalysisManager.cc @@ -1,4 +1,3 @@ -#include #include #include #include @@ -67,6 +66,8 @@ AnalysisManager::AnalysisManager() fTrk = nullptr; fPrim = nullptr; fPixelHitsTree = nullptr; + fTauTree = nullptr; + fCharmTree = nullptr; // fActsParticlesTree = nullptr; fSaveTrack = false; @@ -148,9 +149,12 @@ void AnalysisManager::bookTrkTree() fTrk->Branch("trackKinE", &trackKinE, "trackKinE/D"); fTrk->Branch("trackNPoints", &trackNPoints, "trackNPoints/I"); fTrk->Branch("trackTheta", &trackTheta, "trackTheta/D"); - // fTrk->Branch("trackPointX", &trackPointX); - // fTrk->Branch("trackPointY", &trackPointY); - // fTrk->Branch("trackPointZ", &trackPointZ); + fTrk->Branch("trackProdX", &trackProdX, "trackProdX/D"); + fTrk->Branch("trackProdY", &trackProdY, "trackProdY/D"); + fTrk->Branch("trackProdZ", &trackProdZ, "trackProdZ/D"); + fTrk->Branch("trackDecayX", &trackDecayX, "trackDecayX/D"); + fTrk->Branch("trackDecayY", &trackDecayY, "trackDecayY/D"); + fTrk->Branch("trackDecayZ", &trackDecayZ, "trackDecayZ/D"); } @@ -231,7 +235,7 @@ void AnalysisManager::bookScintTrees() fScintTree->Branch("fromPrimaryLepton", &fScintFromPrimaryLepton); } - +//// --- FASER tracking spectrometer hits --- void AnalysisManager::bookFaserTree() { fFile->cd(fHits->GetName()); @@ -252,6 +256,46 @@ void AnalysisManager::bookFaserTree() fFaserHitsTree->Branch("charge", &fFaserCharge); } +//// --- Children of primary tau decay --- +void AnalysisManager::bookTauTree() +{ + fTauTree = new TTree("tau", "Children of tau decay"); + fTauTree->Branch("evtID", &fPythiaEventID, "evtID/i"); + fTauTree->Branch("PID", &fPythiaPID, "PID/I"); + fTauTree->Branch("TID", &fPythiaTID, "TID/I"); + fTauTree->Branch("PDG", &fPythiaTPDG, "PDG/I"); + fTauTree->Branch("ParentPDG", &fPythiaPPDG, "ParentPDG/I"); + fTauTree->Branch("ProdX", &fPythiaProdX, "ProdX/D"); + fTauTree->Branch("ProdY", &fPythiaProdY, "ProdY/D"); + fTauTree->Branch("ProdZ", &fPythiaProdZ, "ProdZ/D"); + fTauTree->Branch("DecayX", &fPythiaDecayX, "DecayX/D"); + fTauTree->Branch("DecayY", &fPythiaDecayY, "DecayY/D"); + fTauTree->Branch("DecayZ", &fPythiaDecayZ, "DecayZ/D"); + fTauTree->Branch("Px", &fPythiaPx, "Px/D"); + fTauTree->Branch("Py", &fPythiaPy, "Py/D"); + fTauTree->Branch("Pz", &fPythiaPz, "Pz/D"); +} + +//// --- Children of primary tau decay --- +void AnalysisManager::bookCharmTree() +{ + fCharmTree = new TTree("charm", "Children of tau decay"); + fCharmTree->Branch("evtID", &fPythiaEventID, "evtID/i"); + fCharmTree->Branch("PID", &fPythiaPID, "PID/I"); + fCharmTree->Branch("TID", &fPythiaTID, "TID/I"); + fCharmTree->Branch("PDG", &fPythiaTPDG, "PDG/I"); + fCharmTree->Branch("ParentPDG", &fPythiaPPDG, "ParentPDG/I"); + fCharmTree->Branch("ProdX", &fPythiaProdX, "ProdX/D"); + fCharmTree->Branch("ProdY", &fPythiaProdY, "ProdY/D"); + fCharmTree->Branch("ProdZ", &fPythiaProdZ, "ProdZ/D"); + fCharmTree->Branch("DecayX", &fPythiaDecayX, "DecayX/D"); + fCharmTree->Branch("DecayY", &fPythiaDecayY, "DecayY/D"); + fCharmTree->Branch("DecayZ", &fPythiaDecayZ, "DecayZ/D"); + fCharmTree->Branch("Px", &fPythiaPx, "Px/D"); + fCharmTree->Branch("Py", &fPythiaPy, "Py/D"); + fCharmTree->Branch("Pz", &fPythiaPz, "Pz/D"); +} + //--------------------------------------------------------------------- //--------------------------------------------------------------------- @@ -270,6 +314,8 @@ void AnalysisManager::BeginOfRun() bookPrimTree(); bookGeomTree(); if (fSaveTrack) bookTrkTree(); + bookTauTree(); + bookCharmTree(); bookHitsTrees(); bookScintTrees(); @@ -287,6 +333,8 @@ void AnalysisManager::EndOfRun() fFile->cd(); fEvt->Write(); fPrim->Write(); + fTauTree->Write(); + fCharmTree->Write(); FillGeomTree(); fGeom->Write(); if (fSaveTrack) fTrk->Write(); @@ -370,6 +418,10 @@ void AnalysisManager::EndOfEvent(const G4Event *event) // FILL PRIMARIES/TRAJECTORIES TREE FillPrimariesTree(event); + + fPythiaEventID = evtID; + FillPythiaTree(event, fTauTree, {15}); + FillPythiaTree(event, fCharmTree, {411, 421, 431, 4122, 4112, 4212, 4222, 4132, 4232}); if(fSaveTrack) FillTrajectoriesTree(event); //----------------------------------------------------------- @@ -490,11 +542,11 @@ void AnalysisManager::FillPrimariesTree(const G4Event *event) primVx, primVy, primVz, primVt, primPx, primPy, primPz,energy)); - G4cout << G4endl; - G4cout << "PrimaryParticleInfo: PDG code " << primPDG << G4endl - << "Particle unique ID : " << primTrackID << G4endl - << "Momentum : (" << primPx << ", " << primPy << ", " << primPz << ") MeV" << G4endl - << "Vertex : (" << primVx << ", " << primVy << ", " << primVz << ") mm" << G4endl; + // G4cout << G4endl; + // G4cout << "PrimaryParticleInfo: PDG code " << primPDG << G4endl + // << "Particle unique ID : " << primTrackID << G4endl + // << "Momentum : (" << primPx << ", " << primPy << ", " << primPz << ") MeV" << G4endl + // << "Vertex : (" << primVx << ", " << primVy << ", " << primVz << ") mm" << G4endl; fPrim->Fill(); } @@ -531,17 +583,34 @@ void AnalysisManager::FillTrajectoriesTree(const G4Event* event) 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() ); - // } + + // Find start and decay points based on z values + if (trackNPoints > 0) { + G4double minZ = std::numeric_limits::max(); + G4double maxZ = std::numeric_limits::lowest(); + G4ThreeVector prodPos, decayPos; + + for (size_t j = 0; j < trackNPoints; ++j) + { + G4ThreeVector pos = trajectory->GetPoint(j)->GetPosition(); + if (pos.z() < minZ) { + minZ = pos.z(); + prodPos = pos; + } + if (pos.z() > maxZ) { + maxZ = pos.z(); + decayPos = pos; + } + } + + trackProdX = prodPos.x(); + trackProdY = prodPos.y(); + trackProdZ = prodPos.z(); + trackDecayX = decayPos.x(); + trackDecayY = decayPos.y(); + trackDecayZ = decayPos.z(); + } fTrk->Fill(); - // trackPointX.clear(); - // trackPointY.clear(); - // trackPointZ.clear(); } G4cout << "Total number of recorded track: " << count_tracks << std::endl; } @@ -721,6 +790,68 @@ void AnalysisManager::FillFaserOutput() fFaserHitsTree->Fill(); } + +void AnalysisManager::FillPythiaTree(const G4Event* event, TTree *tree, std::vector pdg_ids) +{ + auto trajectoryContainer = event->GetTrajectoryContainer(); + if (!trajectoryContainer) + { + G4cout << "No tracks found: did you enable their storage with '/tracking/storeTrajectory 1'?" << G4endl; + return; + } + for (size_t parentIdx = 0; parentIdx < trajectoryContainer->entries(); ++parentIdx) + { + // Check for primary tau /charm + G4Trajectory *parent = static_cast((*trajectoryContainer)[parentIdx]); + int parentPDG = std::abs(parent->GetPDGEncoding()); + if ((parent->GetParentID() == 0) && (std::find(pdg_ids.begin(), pdg_ids.end(), std::abs(parentPDG)) != pdg_ids.end())) { + G4int parentID = parent->GetTrackID(); + fPythiaPPDG = parent->GetPDGEncoding(); + // std::cout << "Found primary tau/charm with pdg " << parentPDG << ", track id " << parentID << std::endl; + // Check for children of this primary tauy lepton + for (size_t trajIdx = 0; trajIdx < trajectoryContainer->entries(); ++trajIdx) { + G4Trajectory *traj = static_cast((*trajectoryContainer)[trajIdx]); + if (traj->GetParentID() == parentID) { + int trajPDG = traj->GetParticleDefinition()->GetPDGEncoding(); + // std::cout << "Found child of primary tau/charm with PDG ID " << trajPDG << std::endl; + // Write out trajecotry information + fPythiaTPDG = traj->GetParticleDefinition()->GetPDGEncoding(); + fPythiaTID = traj->GetTrackID(); + fPythiaPID = traj->GetParentID(); + G4int trackNPoints = traj->GetPointEntries(); + // Get production and decay vertex + G4double minZ = std::numeric_limits::max(); + G4double maxZ = std::numeric_limits::lowest(); + G4ThreeVector prodPos, decayPos; + for (size_t j = 0; j < trackNPoints; ++j) + { + G4ThreeVector pos = traj->GetPoint(j)->GetPosition(); + if (pos.z() < minZ) { + minZ = pos.z(); + prodPos = pos; + } + if (pos.z() > maxZ) { + maxZ = pos.z(); + decayPos = pos; + } + } + fPythiaProdX = prodPos.x(); + fPythiaProdY = prodPos.y(); + fPythiaProdZ = prodPos.z(); + fPythiaDecayX = decayPos.x(); + fPythiaDecayY = decayPos.y(); + fPythiaDecayZ = decayPos.z(); + G4ThreeVector momentum = traj->GetInitialMomentum(); + fPythiaPx = momentum.x(); + fPythiaPy = momentum.y(); + fPythiaPz = momentum.z(); + tree->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 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/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 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);