diff --git a/NPSimulation/Detectors/Vamos/ChargeStateDistribution.cc b/NPSimulation/Detectors/Vamos/ChargeStateDistribution.cc new file mode 100644 index 0000000000000000000000000000000000000000..f9625aaf2165ae4e5366c86e3c2e0e8aa567c45e --- /dev/null +++ b/NPSimulation/Detectors/Vamos/ChargeStateDistribution.cc @@ -0,0 +1,163 @@ +/***************************************************************************** + * Copyright (C) 2009-2016 this file is part of the NPTool Project * + * * + * For the licensing terms see $NPTOOL/Licence/NPTool_Licence * + * For the list of contributors see $NPTOOL/Licence/Contributors * + *****************************************************************************/ + +/***************************************************************************** + * Original Author: Elia Pilotto, Omar Nasr * + * contact address: pilottoelia@gmail.com, omar.nasr@etu.unicaen.fr * + * * + * Creation Date : September 2021 * + * Last update : * + *---------------------------------------------------------------------------* + * Decription: * + * Use to kill the beam track and replace it with the reaction product * + * * + * * + * * + *---------------------------------------------------------------------------* + * Comment: * + * * + *****************************************************************************/ + +//C++ libraries +#include <iostream> +#include <string> +#include <iomanip> +#include <cmath> +#include <Randomize.hh> +#include <fstream> +//G4 libraries +#include "G4Electron.hh" +#include "G4Gamma.hh" +#include "G4IonTable.hh" +#include "G4EmCalculator.hh" +#include "G4VPhysicalVolume.hh" +#include "G4IonTable.hh" +#include "G4UserLimits.hh" +#include "G4SystemOfUnits.hh" +#include "G4PhysicalConstants.hh" +//nptool libraries +#include "NPFunction.h" +#include "NPInputParser.h" +#include "NPOptionManager.h" +#include "NPSFunction.hh" +//other +#include "ChargeStateDistribution.hh" +#include "RootOutput.h" +#include "TLorentzVector.h" + +using namespace std; + +//////////////////////////////////////////////////////////////////////////////// +NPS::ChargeStateDistribution::ChargeStateDistribution(G4String modelName, G4Region* envelope) + : G4VFastSimulationModel(modelName, envelope) { + + m_StepSize = 5*mm; +} + +//////////////////////////////////////////////////////////////////////////////// +NPS::ChargeStateDistribution::ChargeStateDistribution(G4String modelName) + : G4VFastSimulationModel(modelName) {} + +//////////////////////////////////////////////////////////////////////////////// +NPS::ChargeStateDistribution::~ChargeStateDistribution() {} + +//////////////////////////////////////////////////////////////////////////////// +G4bool NPS::ChargeStateDistribution::IsApplicable(const G4ParticleDefinition& particleType) { + if (particleType.GetPDGCharge() < 10) return false; + else return true; +} + +//////////////////////////////////////////////////////////////////////////////// +G4bool NPS::ChargeStateDistribution::ModelTrigger(const G4FastTrack& fastTrack) { + + cout << "////////////////// MODEL TRIGGER ///////////////////" << endl; + const G4Track* PrimaryTrack = fastTrack.GetPrimaryTrack(); + + G4ThreeVector V = PrimaryTrack->GetMomentum().unit(); + G4ThreeVector P = fastTrack.GetPrimaryTrackLocalPosition(); + G4VSolid* solid = fastTrack.GetPrimaryTrack()->GetVolume()->GetLogicalVolume()->GetSolid(); + double to_exit = solid->DistanceToOut(P, V); + double to_entrance = solid->DistanceToOut(P, -V); + bool is_first = (to_entrance == 0); + bool is_end = (to_exit == 0); + + if (is_first && m_shoot) { + m_shoot = false; + } + + if (is_first) { + m_shoot = true; + } + + cout << "m_StepSize= " << m_StepSize << endl; + cout << "to_entrance= " << to_entrance << endl; + cout << "to_exit= " << to_exit << endl; + cout << "is_first= " << is_first << endl; + cout << "m_shoot= " << m_shoot << endl; + + cout << "Model trigger Q = " << PrimaryTrack->GetDynamicParticle()->GetCharge() << endl; + return m_shoot; +} + +//////////////////////////////////////////////////////////////////////////////// +void NPS::ChargeStateDistribution::DoIt(const G4FastTrack& fastTrack, G4FastStep& fastStep) { + + m_shoot = false; + cout << "////////////////// DO IT ///////////////////" << endl; + G4ThreeVector localPosition = fastTrack.GetPrimaryTrackLocalPosition(); + G4ThreeVector localDirection = fastTrack.GetPrimaryTrackLocalDirection(); + const G4Track* track = fastTrack.GetPrimaryTrack(); + + G4double time = track->GetGlobalTime(); + G4ThreeVector pdirection = track->GetMomentum().unit(); + G4ThreeVector Momentum = track->GetMomentum(); + G4ThreeVector worldPosition = track->GetPosition(); + G4double originalKineticEnergy = track->GetKineticEnergy(); + + G4int Z = track->GetParticleDefinition()->GetAtomicNumber(); + G4int A = track->GetParticleDefinition()->GetAtomicMass(); + G4double Q = track->GetDynamicParticle()->GetCharge(); + + //G4DynamicParticle* dynamicParticle = const_cast<G4DynamicParticle*>(fastTrack.GetPrimaryTrack()->GetDynamicParticle()); + //dynamicParticle->SetCharge(Q + 5); // Modify the particle's charge + + //G4double newQ = track->GetDynamicParticle()->GetCharge(); + + //const G4DynamicParticle* dynamicParticle = fastTrack.GetPrimaryTrack()->GetDynamicParticle(); + //G4ParticleDefinition* particleDef = const_cast<G4ParticleDefinition*>(dynamicParticle->GetDefinition()); + + /* + fastStep.SetPrimaryTrackFinalKineticEnergy(originalKineticEnergy); + fastStep.ProposePrimaryTrackFinalMomentumDirection(localDirection,true); + fastStep.SetPrimaryTrackFinalPosition(localPosition+G4ThreeVector(0,0,1),true); + fastStep.ProposePrimaryTrackFinalTime(time); + fastStep.SetTotalEnergyDeposited(0); + fastStep.ProposePrimaryTrackPathLength(m_StepSize); + */ + + // Set the end of the step conditions + fastStep.SetPrimaryTrackFinalKineticEnergyAndDirection(0, localDirection); + fastStep.SetPrimaryTrackFinalPosition(localPosition+G4ThreeVector(0,0,1),true); + fastStep.SetTotalEnergyDeposited(0); + fastStep.SetPrimaryTrackFinalTime(time); // FIXME + fastStep.KillPrimaryTrack(); + fastStep.SetPrimaryTrackPathLength(0.0); + + G4int newCharge = Q+5; + static G4IonTable* IonTable = G4ParticleTable::GetParticleTable()->GetIonTable(); + G4ParticleDefinition* particleDef; + particleDef = IonTable->GetIon(Z,A); + + G4DynamicParticle dynamicParticle(particleDef,Momentum.unit(),originalKineticEnergy); + fastStep.CreateSecondaryTrack(dynamicParticle,localPosition+G4ThreeVector(0,0,2),time); + + cout << Q << " " << dynamicParticle.GetCharge() << endl; + + return; +} + + diff --git a/NPSimulation/Detectors/Vamos/ChargeStateDistribution.hh b/NPSimulation/Detectors/Vamos/ChargeStateDistribution.hh new file mode 100644 index 0000000000000000000000000000000000000000..e10e2f659e408663fcae93ba26237a15ca178ada --- /dev/null +++ b/NPSimulation/Detectors/Vamos/ChargeStateDistribution.hh @@ -0,0 +1,67 @@ +/***************************************************************************** + * Copyright (C) 2009-2016 this file is part of the NPTool Project * + * * + * For the licensing terms see $NPTOOL/Licence/NPTool_Licence * + * For the list of contributors see $NPTOOL/Licence/Contributors * + *****************************************************************************/ + +/***************************************************************************** + * Original Author: Pierre Morfouace * + * contact address: pierre.morfouace@cea.fr * + * * + * Creation Date : January 2023 * + * Last update : * + *---------------------------------------------------------------------------* + * Decription: Based On SamuraiFieldPropagation * + * Use to kill the beam track and replace it with the reaction product * + * * + * * + * * + *---------------------------------------------------------------------------* + * Comment: * + * * + *****************************************************************************/ + +#ifndef ChargeStateDistribution_h +#define ChargeStateDistribution_h + +#include "G4VFastSimulationModel.hh" +#include "G4Abla.hh" +#include "G4AblaInterface.hh" +#include "G4Fragment.hh" +#include "G4VPhysicalVolume.hh" +#include "G4Region.hh" + +#include "GladFieldMap.h" + +#include <string> +#include <vector> +#include <fstream> +#include <iomanip> + +namespace NPS{ + + class ChargeStateDistribution : public G4VFastSimulationModel{ + public: + ChargeStateDistribution (G4String, G4Region*); + ChargeStateDistribution (G4String); + ~ChargeStateDistribution (); + + public: + G4bool IsApplicable(const G4ParticleDefinition&); + G4bool ModelTrigger(const G4FastTrack &); + void DoIt(const G4FastTrack&, G4FastStep&); + + void RungeKuttaPropagation (const G4FastTrack& fastTrack, G4FastStep& fastStep); + + public: + void SetStepSize(double step){m_StepSize=step;}; + + private: + double m_StepSize; + bool m_shoot; + }; +} + + +#endif diff --git a/NPSimulation/Detectors/Vamos/Vamos.cc b/NPSimulation/Detectors/Vamos/Vamos.cc index 3c3855ab008769f65e103dca3e82267eed6bca8f..582eb93d19f42834628a871f6154e1ff70d8492c 100644 --- a/NPSimulation/Detectors/Vamos/Vamos.cc +++ b/NPSimulation/Detectors/Vamos/Vamos.cc @@ -38,6 +38,9 @@ #include "G4Transform3D.hh" #include "G4UnionSolid.hh" #include "G4VisAttributes.hh" +#include "G4RegionStore.hh" +#include "G4FastSimulationManager.hh" +#include "G4UserLimits.hh" // NPTool header #include "CalorimeterScorers.hh" @@ -49,6 +52,7 @@ #include "NPSHitsMap.hh" #include "RootOutput.h" #include "Vamos.hh" +#include "ChargeStateDistribution.hh" // CLHEP header #include "CLHEP/Random/RandGauss.h" @@ -108,8 +112,8 @@ namespace Vamos_NS { const G4double TMW2_Thickness = 5 * mm; // FPMW - const G4double FPMW1_Width = 200 * mm; - const G4double FPMW1_Length = 400 * mm; + const G4double FPMW1_Width = 70 * mm; + const G4double FPMW1_Length = 95 * mm; const G4double FPMW1_Thickness = 5 * mm; @@ -139,7 +143,8 @@ Vamos::Vamos() { m_DC4 = 0; m_DC1 = 0; m_DC2 = 0; - + m_ChargeStateRegion = 0; + m_StepSize = 5*mm; ICcounter = 0; // RGB Color + Transparency @@ -315,7 +320,7 @@ G4AssemblyVolume* Vamos::BuildTMW1() { G4LogicalVolume* mylar_volume = new G4LogicalVolume(box_mylar, MylarMaterial, "logic_Vamos_mylar1", 0, 0, 0); mylar_volume->SetVisAttributes(m_VisMylar); - gas_volume->SetSensitiveDetector(m_InterScorer); + //gas_volume->SetSensitiveDetector(m_InterScorer); G4ThreeVector pos_gas = G4ThreeVector(0,0,0); m_TMW1->AddPlacedVolume(gas_volume, pos_gas, 0); @@ -355,6 +360,23 @@ G4AssemblyVolume* Vamos::BuildTMW2() { pos_mylar = G4ThreeVector(0,0,0.5*TMW2_Thickness + 0.5*Mylar_Thickness); m_TMW2->AddPlacedVolume(mylar_volume, pos_mylar, 0); + + /*m_ChargeStateRegion = G4RegionStore::GetInstance()->FindOrCreateRegion("ChargeStateProcess"); + m_ChargeStateRegion->AddRootLogicalVolume(gas_volume); + m_ChargeStateRegion->SetUserLimits(new G4UserLimits(m_StepSize)); + G4FastSimulationManager* mng = m_ChargeStateRegion->GetFastSimulationManager(); + + unsigned int size = m_ChargeStateModel.size(); + for(unsigned int i=0; i<size; i++) + mng->RemoveFastSimulationModel(m_ChargeStateModel[i]); + + m_ChargeStateModel.clear(); + G4VFastSimulationModel* fsm; + fsm = new NPS::ChargeStateDistribution("ChargeStateDistribution",m_ChargeStateRegion); + //((NPS::ChargeStateDistribution*) fsm)->SetStepSize(m_StepSize); + m_ChargeStateModel.push_back(fsm);*/ + + } return m_TMW2; } @@ -368,23 +390,41 @@ G4AssemblyVolume* Vamos::BuildFPMW1() { G4Material* DetectorMaterial = MaterialManager::getInstance()->GetGasFromLibrary(m_Gas_FPMW1, m_Pressure_FPMW1, 295 * kelvin); G4Material* MylarMaterial = MaterialManager::getInstance()->GetMaterialFromLibrary("Mylar"); - - G4LogicalVolume* gas_volume = new G4LogicalVolume(box, DetectorMaterial, "logic_Vamos_FPMW1", 0, 0, 0); - gas_volume->SetVisAttributes(m_VisGasC4H10); - + G4Material* Vaccuum = MaterialManager::getInstance()->GetMaterialFromLibrary("Vaccuum"); + + G4LogicalVolume* gas_volume_FPMW1 = new G4LogicalVolume(box, Vaccuum, "logic_Vamos_FPMW1", 0, 0, 0); + gas_volume_FPMW1->SetVisAttributes(m_VisGasC4H10); + gas_volume_FPMW1->SetSensitiveDetector(m_InterScorer); + G4LogicalVolume* mylar_volume = new G4LogicalVolume(box_mylar, MylarMaterial, "logic_Vamos_mylar3", 0, 0, 0); mylar_volume->SetVisAttributes(m_VisMylar); - gas_volume->SetSensitiveDetector(m_InterScorer); G4ThreeVector pos_gas = G4ThreeVector(0,0,0); - m_FPMW1->AddPlacedVolume(gas_volume, pos_gas, 0); + m_FPMW1->AddPlacedVolume(gas_volume_FPMW1, pos_gas, 0); + + //G4ThreeVector pos_mylar = G4ThreeVector(0,0,-0.5*FPMW1_Thickness - 0.5*Mylar_Thickness); + //m_FPMW1->AddPlacedVolume(mylar_volume, pos_mylar, 0); + + //pos_mylar = G4ThreeVector(0,0,0.5*FPMW1_Thickness + 0.5*Mylar_Thickness); + //m_FPMW1->AddPlacedVolume(mylar_volume, pos_mylar, 0); + + + /*m_ChargeStateRegion = G4RegionStore::GetInstance()->FindOrCreateRegion("ChargeStateProcess"); + m_ChargeStateRegion->AddRootLogicalVolume(gas_volume_FPMW1); + m_ChargeStateRegion->SetUserLimits(new G4UserLimits(m_StepSize)); + G4FastSimulationManager* mng = m_ChargeStateRegion->GetFastSimulationManager(); + + unsigned int size = m_ChargeStateModel.size(); + for(unsigned int i=0; i<size; i++) + mng->RemoveFastSimulationModel(m_ChargeStateModel[i]); - /*G4ThreeVector pos_mylar = G4ThreeVector(0,0,-0.5*FPMW1_Thickness - 0.5*Mylar_Thickness); - m_FPMW1->AddPlacedVolume(mylar_volume, pos_mylar, 0); + m_ChargeStateModel.clear(); + G4VFastSimulationModel* fsm; + fsm = new NPS::ChargeStateDistribution("ChargeStateDistribution",m_ChargeStateRegion); + ((NPS::ChargeStateDistribution*) fsm)->SetStepSize(m_StepSize); + m_ChargeStateModel.push_back(fsm);*/ - pos_mylar = G4ThreeVector(0,0,0.5*FPMW1_Thickness + 0.5*Mylar_Thickness); - m_FPMW1->AddPlacedVolume(mylar_volume, pos_mylar, 0);*/ } return m_FPMW1; } diff --git a/NPSimulation/Detectors/Vamos/Vamos.hh b/NPSimulation/Detectors/Vamos/Vamos.hh index d191da9c253895280a4d215931e0715ac11124ed..bed3a01b94549e5b4df5a46571c15d2acde8e373 100644 --- a/NPSimulation/Detectors/Vamos/Vamos.hh +++ b/NPSimulation/Detectors/Vamos/Vamos.hh @@ -33,6 +33,7 @@ using namespace std; #include "G4AssemblyVolume.hh" #include "G4MultiFunctionalDetector.hh" #include "G4UnionSolid.hh" +#include "G4VFastSimulationModel.hh" // NPTool header #include "NPSVDetector.hh" @@ -135,6 +136,12 @@ class Vamos : public NPS::VDetector{ private: TVamosData* m_Event ; + private: + G4Region* m_ChargeStateRegion; + vector<G4VFastSimulationModel*> m_ChargeStateModel; + + int m_StepSize; + //////////////////////////////////////////////////// ///////////////Private intern Data////////////////// //////////////////////////////////////////////////// diff --git a/Projects/PISTA/PISTA.detector b/Projects/PISTA/PISTA.detector deleted file mode 100644 index 8b8aa25d53e90fb32d517840cb31e7b6a5ba2cb9..0000000000000000000000000000000000000000 --- a/Projects/PISTA/PISTA.detector +++ /dev/null @@ -1,137 +0,0 @@ -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -Target - THICKNESS= 0.44 micrometer - RADIUS= 20 mm - MATERIAL= C - ANGLE= 0 deg - X= 0 mm - Y= 0 mm - Z= 0 mm -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -PISTA - R= 100 mm - THETA= 39 deg - PHI= 315 deg -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -PISTA - R= 100 mm - THETA= 39 deg - PHI= 270 deg -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -PISTA - R= 100 mm - THETA= 39 deg - PHI= 225 deg -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -PISTA - R= 100 mm - THETA= 39 deg - PHI= 180 deg -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -PISTA - R= 100 mm - THETA= 39 deg - PHI= 135 deg -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -PISTA - R= 100 mm - THETA= 39 deg - PHI= 90 deg -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -PISTA - R= 100 mm - THETA= 39 deg - PHI= 45 deg -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -PISTA - R= 100 mm - THETA= 39 deg - PHI= 0 deg -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% Vamos Configuration -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -%Vamos -% R= 0.3 m -% Theta= 0 deg -% -%Vamos DC -% Z= 60 mm -% Gas= iC4H10 -% Pressure= 0.01 bar -% Temperature= 295 kelvin -% -%Vamos DC -% Z= 200 mm -% Gas= iC4H10 -% Pressure= 0.01 bar -% Temperature= 295 kelvin -% -%Vamos BeamCatcher -% Material= BC400 -% Width= 30 mm -% Length= 80 mm -% Thickness= 30 mm -% Pos= 0 0 600 mm -% -%Vamos MWPPAC -% Z= 750 mm -% Gas= iC4H10 -% Pressure= 0.01 bar -% Temperature= 295 kelvin -% -%Vamos DC -% Z= 940 mm -% Gas= iC4H10 -% Pressure= 0.01 bar -% Temperature= 295 kelvin -% -%Vamos DC -% Z= 1060 mm -% Gas= iC4H10 -% Pressure= 0.01 bar -% Temperature= 295 kelvin -% -%Vamos IC -% Z= 1200 mm -% Thickness= 50 mm -% Gas= CF4 -% Pressure= 0.2 bar -% Temperature= 295 kelvin -% -%Vamos IC -% Z= 1250 mm -% Thickness= 50 mm -% Gas= CF4 -% Pressure= 0.2 bar -% Temperature= 295 kelvin -% -%Vamos IC -% Z= 1300 mm -% Thickness= 50 mm -% Gas= CF4 -% Pressure= 0.2 bar -% Temperature= 295 kelvin -% -%Vamos IC -% Z= 1375 mm -% Thickness= 100 mm -% Gas= CF4 -% Pressure= 0.2 bar -% Temperature= 295 kelvin -% -%Vamos IC -% Z= 1525 mm -% Thickness= 200 mm -% Gas= CF4 -% Pressure= 0.2 bar -% Temperature= 295 kelvin -% -%Vamos IC -% Z= 1725 mm -% Thickness= 200 mm -% Gas= CF4 -% Pressure= 0.2 bar -% Temperature= 295 kelvin -% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%