diff --git a/FullSimLight/fullSimLight.cc b/FullSimLight/fullSimLight.cc index db9a07a9dff5ab6b2f3f30d3cf868a095b4b78dc..e5210d9fbbe976fbf15853b5d023ca83b8106b6a 100644 --- a/FullSimLight/fullSimLight.cc +++ b/FullSimLight/fullSimLight.cc @@ -53,7 +53,6 @@ static G4double parDeltaChord; //DELTACHORD LENGTH void GetInputArguments(int argc, char** argv); void Help(); - int main(int argc, char** argv) { // JFB if the G4 environment does not already set path to these variables, look for diff --git a/FullSimLight/include/FSLDetectorConstruction.hh b/FullSimLight/include/FSLDetectorConstruction.hh index c272c78334e4392c7657dcdf3a3eddafa7c76c62..e580ee316642923317feb35f55ce583fb4455d0b 100644 --- a/FullSimLight/include/FSLDetectorConstruction.hh +++ b/FullSimLight/include/FSLDetectorConstruction.hh @@ -19,12 +19,9 @@ // #include "G4UIcmdWithADoubleAndUnit.hh" #include "G4FieldManager.hh" -#include "G4ChordFinder.hh" // - #include <iomanip> -#include "G4ChordFinder.hh" #include "G4SystemOfUnits.hh" #include "G4Field.hh" #include "G4Mag_UsualEqRhs.hh" @@ -56,6 +53,8 @@ #include <cassert> #include "HistoManager.hh" #include "HistoMessenger.hh" +#include "G4Track.hh" +#include "G4SteppingManager.hh" // // Units @@ -101,6 +100,8 @@ public: void SetDumpGDML(const bool dumpGDML) {fDumpGDML=dumpGDML;} void SetDeltaChord(const G4double &parDeltaChord) {fDeltaChord=parDeltaChord;} //////////added void SetStepMax(G4double theStepMax) { fStepMax = theStepMax; }//////////added + void TrackAndPrintParticleDetails(const G4Track* track);//////added + void UserSteppingAction(const G4Step* step); /////////Changes/////////// //void SetDeltaChord(const std::string DeltaChord); /// Common method to construct a driver with a stepper of requested type. diff --git a/FullSimLight/src/FSLDetectorConstruction.cc b/FullSimLight/src/FSLDetectorConstruction.cc index 9aaff2e9772013fdbc34e244f590bfedbf4d3103..a84863a41266299e38370fbf8b1dfa912db17f9d 100644 --- a/FullSimLight/src/FSLDetectorConstruction.cc +++ b/FullSimLight/src/FSLDetectorConstruction.cc @@ -145,7 +145,13 @@ #include "G4RegionStore.hh" #include "G4Region.hh" - +#include "G4Track.hh" +#include "G4ParticleDefinition.hh" +#include "G4DynamicParticle.hh" +#include "G4SteppingManager.hh" +#include "G4Track.hh" +#include "G4Step.hh" +#include "G4ThreeVector.hh" // @@ -177,19 +183,7 @@ FSLDetectorConstruction::~FSLDetectorConstruction() { delete fDetectorMessenger; } -/*///////////////////// -FSLDetectorConstruction::FSLDetectorConstruction() : fWorld(nullptr), fDetectorMessenger(nullptr) -{ - // ... other member variable initializations - // Create the HistoManager and configure it - HistoManager* histoManager = HistoManager::GetPointer(); - histoManager->SetFileName("y.root"); // Set the output file name - histoManager->SetVerbose(1); // Set verbosity level as needed - -} - -*///////////////////// GeoPhysVol* FSLDetectorConstruction::CreateTheWorld(GeoPhysVol* world) { @@ -217,28 +211,7 @@ GeoPhysVol* FSLDetectorConstruction::CreateTheWorld(GeoPhysVol* world) } return world; } -/*/////////////////////////////////////////////////////////////////////////////// -G4VPhysicalVolume* FSLDetectorConstruction::Construct() -{ - // ... existing code - - // Your code for creating the detector geometry - - // Example: Start of the event - HistoManager* histoManager = HistoManager::GetPointer(); - histoManager->BeginOfEvent(); - - // ... existing code for building the detector geometry - - // Example: End of the event - histoManager->EndOfEvent(); - - // ... remaining code - - return fWorld; -} -*////////////////////////////////////////////////////////////////////////////// G4VPhysicalVolume *FSLDetectorConstruction::Construct() { fTimer.Start(); @@ -341,68 +314,6 @@ G4VPhysicalVolume *FSLDetectorConstruction::Construct() fWorld = physWorld; fWorld->GetLogicalVolume()->SetVisAttributes(G4VisAttributes::GetInvisible()); -//G4LogicalVolumeStore* logVolStore = G4LogicalVolumeStore::GetInstance(); - -/* G4RegionStore* regionStore = G4RegionStore::GetInstance(); - G4out <<"list of region if it worked:" << G4endl; - G4RegionStore* rv = (*regionStore)[i]; - G4String regionName = rv->GetName(); - - // Print the logical volume's name and information - G4cout << "Logical Volume " << i << ": " << regionName << G4endl; */ -/* -G4cout << "List of all logical volumes:" << G4endl; -for (size_t i = 0; i < logVolStore->size(); ++i) { - G4LogicalVolume* lv = (*logVolStore)[i]; - G4String volumeName = lv->GetName(); - - // Print the logical volume's name and information - G4cout << "Logical Volume " << i << ": " << volumeName << G4endl; - // You can print additional information about the logical volume here - - // If you need to access other properties of the logical volume, you can do so here - // For example, you can get the material or other properties: - // G4Material* material = lv->GetMaterial(); - // G4cout << " Material: " << material->GetName() << G4endl; - - // If you need to access the daughter volumes, you can iterate through them as well: - for (size_t j = 0; j < lv->GetNoDaughters(); ++j) { - G4VPhysicalVolume* daughterVolume = lv->GetDaughter(j); - G4String daughterName = daughterVolume->GetName(); - G4cout << " Daughter Volume " << j << ": " << daughterName << G4endl; - // Print additional information about the daughter volume if needed - } -} -*/ -// ... (rest of your code) -// Include necessary Geant4 headers -// ... - -/* - if (regionStore) { - G4cout << "List of defined regions:" << G4endl; - G4RegionStore::iterator it; - for (it = regionStore->begin(); it != regionStore->end(); ++it) { - G4Region* region = *it; - G4cout << "Region Name: " << region->GetName() << G4endl; - - // Access region properties - G4bool isEnabled = region->IsEnabled(); - G4bool isProductionCutSet = region->IsProductionCutSet(); - G4double productionCut = region->GetProductionCut(); - - G4cout << " Is Enabled: " << isEnabled << G4endl; - G4cout << " Is Production Cut Set: " << isProductionCutSet << G4endl; - G4cout << " Production Cut: " << productionCut << G4endl; - } - } else { - G4cout << "No regions defined." << G4endl; - } -} - -// ... -*/ /////////////////////// - if (fWorld == 0) { G4ExceptionDescription ed; ed << "World volume not set properly check your setup selection criteria or GDML input!" << G4endl; @@ -583,14 +494,9 @@ if (fFieldConstant && std::abs(fFieldValue) > 0.0) { G4FieldManager *fieldMgr = G4TransportationManager::GetTransportationManager()->GetFieldManager(); fieldMgr->SetDetectorField(uniformMagField); fieldMgr->CreateChordFinder(uniformMagField); -////////STEPMAX///// -//fieldMgr->SetMaximumTrackLength(fStepMax); // Set the maximum step length for particle tracking -// Set the maxStep parameter for particle tracking - // fieldMgr->SetMaximumTrackLength(5 * cm); // Adjust maxStep as needed - // fieldMgr->GetChordFinder()->SetDeltaChord( fDeltaChord); ///////////////////////////////////////////DeltaChord(mm) Parameter Set//////////////////// -//fieldMgr->GetChordFinder()->SetDeltaChord( 0.001); //DeltaChord(mm) Parameter Set 1 +//fieldMgr->GetChordFinder()->SetDeltaChord( 0.001); //DeltaChord(mm) Parametet 1 //fieldMgr->GetChordFinder()->SetDeltaChord( 0.1); //DeltaChord(mm) Parameter Set 2 //fieldMgr->GetChordFinder()->SetDeltaChord( 0.002); //DeltaChord(mm) Parameter Set 3/ generaliseremove value and put name //fDeltaChord = std::stod(fDeltaChord); @@ -614,14 +520,157 @@ else { fField.Put(g4Field); } } +///////////ChordFinder////////// +/* G4FieldManager *fieldMgr = G4TransportationManager::GetTransportationManager()->GetFieldManager(); + fieldMgr->CreateChordFinder(); + fieldMgr->GetChordFinder()->SetDeltaChord( fDeltaChord); -// -// Define the maximum step length -//G4double theStepMax = 1.0 * mm; // Set your desired maximum step length (e.g., 1 mm) -G4double theStepMax; -// Get the logical volume store -G4LogicalVolumeStore* logVolStore = G4LogicalVolumeStore::GetInstance(); +// Getting the chord finder and setting delta chord +G4ChordFinder* chordFinder = fieldMgr->GetChordFinder(); +chordFinder->SetDeltaChord(fDeltaChord); + +// Assuming m_currChordFinder is a pointer to your chord finder class +fDeltaChord->SetDeltaChord(); +fDeltaChord->SetDeltaOneStep(); +fDeltaChord->SetDeltaIntersection(); + + +////Tracker from CMS//////// +void G4FieldManager::ConfigureForTrack(const G4Track *ptrack) { + if (ptrack->GetKineticEnergy() // G4doc. < 2.5 * MeV )// { + TrackAndPrintParticleDetails(ptrack); + + if (!Tracker) { + setChordFinderForTracker(); + } + } else { + // restore defaults + setDefaultChordFinder(); + } +} +*/ + +/*void G4UserSteppingAction::UserSteppingAction(const G4Step* step) { + G4Track* track = step->GetTrack(); + G4ParticleDefinition* particle = track->GetDefinition(); + G4ThreeVector position = track->GetPosition(); + G4ThreeVector momentum = track->GetMomentum(); + + G4cout << "Particle: " << particle->GetParticleName() << G4endl; + G4cout << "Position: " << position << G4endl; + G4cout << "Momentum: " << momentum << G4endl; + // Add more information as needed + + // Call base class method for processing further + G4UserSteppingAction::UserSteppingAction(step); +*/ +/*void G4FieldManager::TrackAndPrintParticleDetails(const G4Track* track) { + const G4ParticleDefinition* particleDefinition = track->GetParticleDefinition(); + const G4DynamicParticle* dynamicParticle = track->GetDynamicParticle(); + + G4cout << "Particle Details:" << G4endl; + G4cout << " Track ID: " << track->GetTrackID() << G4endl; + G4cout << " Particle Name: " << particleDefinition->GetParticleName() << G4endl; + G4cout << " Particle Mass: " << particleDefinition->GetPDGMass() / GeV << " GeV/c^2" << G4endl; + G4cout << " Kinetic Energy: " << track->GetKineticEnergy() / GeV << " GeV" << G4endl; + G4cout << " Momentum: " << track->GetMomentum() / GeV << " GeV/c" << G4endl; + G4cout << " Position: " << track->GetPosition() / mm << " mm" << G4endl; + G4cout << " Direction: " << track->GetMomentumDirection() << G4endl; + + G4cout << G4endl; +}*/ + + ///////////////////////////////////////////////////////////// +///////////// MuonSys ///////////// +/*G4LogicalVolume* muonSysLogicalVolume = nullptr; + +// Loop through the logical volumes to find "MuonSys" by name +for (auto lv : *logVolStore) { + if (lv->GetName() == "MuonSys") { + muonSysLogicalVolume = lv; + break; + } +} +////////////////////DeltaIntersection-MuonSys///////////// +// Check for "MuonSys" +if (muonSysLogicalVolume) { + // Create a G4Region and associate it with "MuonSys" logical volume + G4Region* muonSysRegion = new G4Region("MuonSysRegion"); + muonSysRegion->AddRootLogicalVolume(muonSysLogicalVolume); + + // Set the step size limit for "MuonSys" logical volume + G4double theStepMaxMuonSys = 1.0 * mm; // StepMax value for MuonSys + G4UserLimits* muonSysUserLimits = muonSysLogicalVolume->GetUserLimits(); + if (!muonSysUserLimits) { + muonSysUserLimits = new G4UserLimits(theStepMaxMuonSys); + muonSysLogicalVolume->SetUserLimits(muonSysUserLimits); + } else { + // If user limits already exist, update the maximum step length + muonSysUserLimits->SetMaxAllowedStep(theStepMaxMuonSys); + } + + // Get the field manager for the "MuonSysRegion" + G4FieldManager* muonSysFieldMgr = muonSysRegion->GetFieldManager(); + + // Set the deltaIntersection for the "MuonSysRegion" + G4double theDeltaIntersectionMuonSys = 0.1 * mm; // DeltaIntersection value for MuonSys + muonSysFieldMgr->SetDeltaIntersection(theDeltaIntersectionMuonSys); + + // Print the step size limit for "MuonSys" logical volume + G4cout << "Step size limit on MuonSys: " << theStepMaxMuonSys / mm << " mm" << G4endl; + // Print the deltaIntersection for "MuonSys" region + G4cout << "DeltaIntersection for MuonSysRegion: " << theDeltaIntersectionMuonSys / mm << " mm" << G4endl; + } else { + G4cerr << "Error: Logical volume 'MuonSys' not found in the geometry description!" << G4endl; + }*/ +//} +///////////////////////////////////////////////////////// +///////////// SCT ///////////// +/* +// Create a list of SCT logical volume names +std::vector<std::string> sctLogicalVolumeNames = {"SCT_Barrel", "SCT_ForwardA", "SCT_ForwardC"}; + +for (const std::string& sctName : sctLogicalVolumeNames) { + G4LogicalVolume* sctLogicalVolume = nullptr; + + // Loop through the logical volumes to find the SCT logical volume by name + for (auto lv : *logVolStore) { + if (lv->GetName() == G4String(sctName.c_str())) { // Convert std::string to G4String + sctLogicalVolume = lv; + break; + } + } + + // Check if the SCT logical volume was found + if (sctLogicalVolume) { + // Create a G4Region and associate it with the SCT logical volume + G4Region* sctRegion = new G4Region((sctName + "Region").c_str()); + sctRegion->AddRootLogicalVolume(sctLogicalVolume); + + // Set the step size limit for the SCT logical volume + G4double theStepMaxSCT = 2.0 * mm; // StepMax value for SCT + G4UserLimits* sctUserLimits = sctLogicalVolume->GetUserLimits(); + if (!sctUserLimits) { + sctUserLimits = new G4UserLimits(theStepMaxSCT); + sctLogicalVolume->SetUserLimits(sctUserLimits); + } else { + // If user limits already exist, update the maximum step length + sctUserLimits->SetMaxAllowedStep(theStepMaxSCT); + } + + // Print the step size limit for the SCT logical volume + G4cout << "Step size limit on " << sctName << ": " << theStepMaxSCT / mm << " mm" << G4endl; + } else { + G4cerr << "Error: Logical volume '" << sctName << "' not found in the geometry description!" << G4endl; + } +} +*/ + +///////////////////// + +/* +//G4double theStepMax = 1.0 * mm; // Loop through logical volumes and set the step limits for (auto lv : *logVolStore) { G4UserLimits* userLimits = lv->GetUserLimits(); @@ -635,9 +684,8 @@ for (auto lv : *logVolStore) { } // Print the value of maxStepLength (theStepMax) -G4cout << "StepMax (theStepMax) set to: " << theStepMax << std::endl; +G4cout << "StepMax (theStepMax) set to: " << theStepMax << std::endl; */ - // This is thread-local //To print the values of parameters G4FieldManager *fieldMgr = G4TransportationManager::GetTransportationManager()->GetFieldManager(); @@ -648,7 +696,8 @@ G4cout << "EpsilonMax: " << fieldMgr->GetMaximumEpsilonStep() / mm << " mm" << G G4cout << "Delta Chord Length: " << fieldMgr->GetChordFinder()->GetDeltaChord() / mm << " mm" << G4endl; } - +// } +//} //#if G4VERSION_NUMBER < 1040 // // auto stepper = getStepper(m_integratorStepper, field);