From 3a863ef94352da66ea1a49e411c8416e8fb4052e Mon Sep 17 00:00:00 2001 From: adrien matta <matta@lpccaen.in2p3.fr> Date: Wed, 21 Dec 2016 17:34:23 +0100 Subject: [PATCH] * Adding new physics facility to simulate gazeous detector (chio,tpc,...) - Follow the optical photon model - Add a property table to the material, with drift velocity and diffusivity - add new particle drift electron - add new transport process for drift electon - Still need some work (drift direction is hard coded, small bug in transport cause G4 to issue exeception) --- Inputs/DetectorConfiguration/csi.detector | 2 +- NPSimulation/Detectors/Chio/CMakeLists.txt | 2 + NPSimulation/Detectors/Chio/Chio.cc | 339 ++++++++++++++++++ NPSimulation/Detectors/Chio/Chio.hh | 117 ++++++ NPSimulation/Process/CMakeLists.txt | 2 +- NPSimulation/Process/G4DEAbsorption.cc | 133 +++++++ NPSimulation/Process/G4DEAbsorption.hh | 127 +++++++ NPSimulation/Process/G4DEAmplification.cc | 200 +++++++++++ NPSimulation/Process/G4DEAmplification.hh | 130 +++++++ NPSimulation/Process/G4DEProcessSubType.hh | 52 +++ NPSimulation/Process/G4DETransport.cc | 165 +++++++++ NPSimulation/Process/G4DETransport.hh | 127 +++++++ NPSimulation/Process/G4DriftElectron.cc | 36 +- .../Process/G4DriftElectronPhysics.cc | 253 +++++++++++++ .../Process/G4DriftElectronPhysics.hh | 116 ++++++ .../Process/G4DriftElectronProcessIndex.hh | 76 ++++ NPSimulation/Process/G4IonizationWithDE.cc | 248 +++++++++++++ NPSimulation/Process/G4IonizationWithDE.hh | 254 +++++++++++++ NPSimulation/Process/PhysicsList.cc | 21 +- NPSimulation/Process/PhysicsList.hh | 5 + Projects/BeLise/34Si.beam | 19 + Projects/BeLise/chio.detector | 15 + 22 files changed, 2416 insertions(+), 23 deletions(-) create mode 100644 NPSimulation/Detectors/Chio/CMakeLists.txt create mode 100644 NPSimulation/Detectors/Chio/Chio.cc create mode 100644 NPSimulation/Detectors/Chio/Chio.hh create mode 100644 NPSimulation/Process/G4DEAbsorption.cc create mode 100644 NPSimulation/Process/G4DEAbsorption.hh create mode 100644 NPSimulation/Process/G4DEAmplification.cc create mode 100644 NPSimulation/Process/G4DEAmplification.hh create mode 100644 NPSimulation/Process/G4DEProcessSubType.hh create mode 100644 NPSimulation/Process/G4DETransport.cc create mode 100644 NPSimulation/Process/G4DETransport.hh create mode 100644 NPSimulation/Process/G4DriftElectronPhysics.cc create mode 100644 NPSimulation/Process/G4DriftElectronPhysics.hh create mode 100644 NPSimulation/Process/G4DriftElectronProcessIndex.hh create mode 100644 NPSimulation/Process/G4IonizationWithDE.cc create mode 100644 NPSimulation/Process/G4IonizationWithDE.hh create mode 100644 Projects/BeLise/34Si.beam create mode 100644 Projects/BeLise/chio.detector diff --git a/Inputs/DetectorConfiguration/csi.detector b/Inputs/DetectorConfiguration/csi.detector index 771f5f7a1..6a0aed081 100644 --- a/Inputs/DetectorConfiguration/csi.detector +++ b/Inputs/DetectorConfiguration/csi.detector @@ -21,7 +21,7 @@ CsI FaceFront= 100 mm FaceBack= 100 mm Radius= 0 mm - Scintillator= CsI + Scintillator= CsI_Scintillator LeadThickness= 0 mm %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% diff --git a/NPSimulation/Detectors/Chio/CMakeLists.txt b/NPSimulation/Detectors/Chio/CMakeLists.txt new file mode 100644 index 000000000..afbcd0e5f --- /dev/null +++ b/NPSimulation/Detectors/Chio/CMakeLists.txt @@ -0,0 +1,2 @@ +add_library(NPSChio SHARED Chio.cc) +target_link_libraries(NPSChio NPSCore ${ROOT_LIBRARIES} ${Geant4_LIBRARIES} ${NPLib_LIBRARIES} -lNPChio) diff --git a/NPSimulation/Detectors/Chio/Chio.cc b/NPSimulation/Detectors/Chio/Chio.cc new file mode 100644 index 000000000..443f6ffd1 --- /dev/null +++ b/NPSimulation/Detectors/Chio/Chio.cc @@ -0,0 +1,339 @@ +/***************************************************************************** + * Copyright (C) 2009-2106 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: Adrien Matta contact address: matta@lpccaen.in2p3.fr * + * * + * Creation Date : December 2016 * + * Last update : * + *---------------------------------------------------------------------------* + * Decription: * + * This class describe Chio simulation * + * * + *---------------------------------------------------------------------------* + * Comment: * + * * + *****************************************************************************/ + +// C++ headers +#include <sstream> +#include <cmath> +#include <limits> +//G4 Geometry object +#include "G4Tubs.hh" +#include "G4Box.hh" + +//G4 sensitive +#include "G4SDManager.hh" +#include "G4MultiFunctionalDetector.hh" + +//G4 various object +#include "G4Material.hh" +#include "G4Transform3D.hh" +#include "G4PVPlacement.hh" +#include "G4VisAttributes.hh" +#include "G4Colour.hh" + +// G4 Field +#include "G4FieldManager.hh" +#include "G4ElectricField.hh" +#include "G4UniformElectricField.hh" +#include "G4TransportationManager.hh" +#include "G4EqMagElectricField.hh" +#include "G4MagIntegratorStepper.hh" +#include "G4ClassicalRK4.hh" +#include "G4MagIntegratorDriver.hh" +#include "G4ChordFinder.hh" +#include "G4MaterialPropertiesTable.hh" + +// NPTool header +#include "Chio.hh" +#include "CalorimeterScorers.hh" +#include "RootOutput.h" +#include "MaterialManager.hh" +#include "NPSDetectorFactory.hh" +#include "NPOptionManager.h" +#include "FastDriftElectron.hh" +// CLHEP header +#include "CLHEP/Random/RandGauss.h" + +using namespace std; +using namespace CLHEP; + + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +namespace Chio_NS{ + // Energy and time Resolution + const double EnergyThreshold = 0.1*MeV; + const double ResoTime = 4.5*ns ; + const double ResoEnergy = 1.0*MeV ; + const double Radius = 50*mm ; + const double Width = 100*mm ; + const double Thickness = 300*mm ; + const string Material = "BC400"; +} +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +// Chio Specific Method +Chio::Chio(){ + m_Event = new TChio_anData() ; + m_ChioScorer = 0; + m_SquareDetector = 0; + m_CylindricalDetector = 0; + + + // RGB Color + Transparency + m_VisChamber = new G4VisAttributes(G4Colour(0.5, 0.5, 0.5, 0.25)); + m_VisWindows = new G4VisAttributes(G4Colour(1, 0, 0, 0.25)); + m_VisGas = new G4VisAttributes(G4Colour(0, 0.5, 0.5, 0.5)); + +} + +Chio::~Chio(){ +} +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +void Chio::AddDetector(G4ThreeVector POS, string Shape){ + // Convert the POS value to R theta Phi as Spherical coordinate is easier in G4 + m_R.push_back(POS.mag()); + m_Theta.push_back(POS.theta()); + m_Phi.push_back(POS.phi()); + m_Shape.push_back(Shape); +} + + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +void Chio::AddDetector(double R, double Theta, double Phi, string Shape){ + m_R.push_back(R); + m_Theta.push_back(Theta); + m_Phi.push_back(Phi); + m_Shape.push_back(Shape); +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +G4LogicalVolume* Chio::BuildDetector(){ + if(!m_SquareDetector){ + // Main volume + G4Box* sChamber = new G4Box("Chio_Box",7*cm*0.5, + 7*cm*0.5,12*cm*0.5+2*12*micrometer*0.5); + + // Gaz volume + G4Box* sGas = new G4Box("Chio_Gas",6*cm*0.5, + 6*cm*0.5,12*cm*0.5-2*12*micrometer*0.5); + + // Entrance/Exit windows + G4Box* sWindows = new G4Box("Chio_Windows",7*cm*0.5, + 7*cm*0.5,12*micrometer*0.5); + + + G4Material* Fe= MaterialManager::getInstance()->GetMaterialFromLibrary("Fe"); + G4Material* CF4= MaterialManager::getInstance()->GetGasFromLibrary("CF4",0.0693276*bar,273.15*kelvin); + G4Material* Mylar= MaterialManager::getInstance()->GetMaterialFromLibrary("Mylar"); + + G4MaterialPropertiesTable* MPT = new G4MaterialPropertiesTable(); + MPT->AddConstProperty("DE_PAIRENERGY",30*eV); + MPT->AddConstProperty("DE_YIELD",1e-3); +// MPT->AddConstProperty("DE_AMPLIFICATION",1e4); + MPT->AddConstProperty("DE_ABSLENGTH",1*km); + MPT->AddConstProperty("DE_DRIFTSPEED",8e-3*mm/ns); + MPT->AddConstProperty("DE_TRANSVERSALSPREAD",6e-5*mm2/ns); + MPT->AddConstProperty("DE_LONGITUDINALSPREAD",4e-5*mm2/ns); + + CF4->SetMaterialPropertiesTable(MPT); + m_SquareDetector = new G4LogicalVolume(sChamber,Fe,"logic_Chio_Box",0,0,0); + G4LogicalVolume* logicGas = new G4LogicalVolume(sGas,CF4,"logic_Gas",0,0,0); + G4LogicalVolume* logicWindows = new G4LogicalVolume(sWindows,Mylar,"logic_Windows",0,0,0); + + G4RotationMatrix* Rot = new G4RotationMatrix(); + + new G4PVPlacement(G4Transform3D(*Rot,G4ThreeVector(0,0,0)), + logicGas, + "ChioGas",m_SquareDetector,false,0); + + new G4PVPlacement(G4Transform3D(*Rot,G4ThreeVector(0,0,6*cm-6*micrometer)), + logicWindows, + "ChioExitWindows",m_SquareDetector,false,0); + + new G4PVPlacement(G4Transform3D(*Rot,G4ThreeVector(0,0,-6*cm+6*micrometer)), + logicWindows, + "ChioEntranceWindows",m_SquareDetector,false,0); + + +/* G4Region* DriftRegion = new G4Region("DriftRegion"); + DriftRegion->AddRootLogicalVolume(logicGas); + G4ProductionCuts* cuts = new G4ProductionCuts; + cuts->SetProductionCut(1e-9*micrometer); // same cuts for gamma, e- and e+ + DriftRegion->SetProductionCuts(cuts); + + new DriftElectron("DriftElectron",DriftRegion);*/ + m_SquareDetector->SetVisAttributes(m_VisChamber); + logicGas->SetVisAttributes(m_VisGas); + logicWindows->SetVisAttributes(m_VisWindows); + m_SquareDetector->SetSensitiveDetector(m_ChioScorer); + } + return m_SquareDetector; +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +// Virtual Method of NPS::VDetector class + +// Read stream at Configfile to pick-up parameters of detector (Position,...) +// Called in DetecorConstruction::ReadDetextorConfiguration Method +void Chio::ReadConfiguration(NPL::InputParser parser){ + vector<NPL::InputBlock*> blocks = parser.GetAllBlocksWithToken("Chio"); + if(NPOptionManager::getInstance()->GetVerboseLevel()) + cout << "//// " << blocks.size() << " detectors found " << endl; + + vector<string> cart = {"POS","Shape"}; + vector<string> sphe = {"R","Theta","Phi","Shape"}; + + for(unsigned int i = 0 ; i < blocks.size() ; i++){ + if(blocks[i]->HasTokenList(cart)){ + if(NPOptionManager::getInstance()->GetVerboseLevel()) + cout << endl << "//// Chio " << i+1 << endl; + + G4ThreeVector Pos = NPS::ConvertVector(blocks[i]->GetTVector3("POS","mm")); + string Shape = blocks[i]->GetString("Shape"); + AddDetector(Pos,Shape); + } + else if(blocks[i]->HasTokenList(sphe)){ + if(NPOptionManager::getInstance()->GetVerboseLevel()) + cout << endl << "//// Chio " << i+1 << endl; + double R = blocks[i]->GetDouble("R","mm"); + double Theta = blocks[i]->GetDouble("Theta","deg"); + double Phi = blocks[i]->GetDouble("Phi","deg"); + string Shape = blocks[i]->GetString("Shape"); + AddDetector(R,Theta,Phi,Shape); + } + else{ + cout << "ERROR: check your input file formatting " << endl; + exit(1); + } + } +} + + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +// Construct detector and inialise sensitive part. +// Called After DetecorConstruction::AddDetector Method +void Chio::ConstructDetector(G4LogicalVolume* world){ + for (unsigned short i = 0 ; i < m_R.size() ; i++) { + + G4double wX = m_R[i] * sin(m_Theta[i] ) * cos(m_Phi[i] ) ; + G4double wY = m_R[i] * sin(m_Theta[i] ) * sin(m_Phi[i] ) ; + G4double wZ = m_R[i] * cos(m_Theta[i] ) ; + G4ThreeVector Det_pos = G4ThreeVector(wX, wY, wZ) ; + // So the face of the detector is at R instead of the middle + Det_pos+=Det_pos.unit()*Chio_NS::Thickness*0.5; + // Building Detector reference frame + G4double ii = cos(m_Theta[i]) * cos(m_Phi[i]); + G4double jj = cos(m_Theta[i]) * sin(m_Phi[i]); + G4double kk = -sin(m_Theta[i]); + G4ThreeVector Y(ii,jj,kk); + G4ThreeVector w = Det_pos.unit(); + G4ThreeVector u = w.cross(Y); + G4ThreeVector v = w.cross(u); + v = v.unit(); + u = u.unit(); + + G4RotationMatrix* Rot = new G4RotationMatrix(u,v,w); + + if(m_Shape[i] == "Square"){ + new G4PVPlacement(G4Transform3D(*Rot,Det_pos), + BuildDetector(), + "Chio",world,false,i+1); + } + } +} +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +// Add Detector branch to the EventTree. +// Called After DetecorConstruction::AddDetector Method +void Chio::InitializeRootOutput(){ + RootOutput *pAnalysis = RootOutput::getInstance(); + TTree *pTree = pAnalysis->GetTree(); + if(!pTree->FindBranch("Chio")){ + pTree->Branch("Chio", "TChioData", &m_Event) ; + } + pTree->SetBranchAddress("Chio", &m_Event) ; +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +// Read sensitive part and fill the Root tree. +// Called at in the EventAction::EndOfEventAvtion +void Chio::ReadSensitive(const G4Event* event){ + m_Event->Clear(); + + /////////// + // Calorimeter scorer + NPS::HitsMap<G4double*>* CaloHitMap; + std::map<G4int, G4double**>::iterator Calo_itr; + + G4int CaloCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("ChioScorer/Calorimeter"); + CaloHitMap = (NPS::HitsMap<G4double*>*)(event->GetHCofThisEvent()->GetHC(CaloCollectionID)); + /* + // Loop on the Calo map + for (Calo_itr = CaloHitMap->GetMap()->begin() ; Calo_itr != CaloHitMap->GetMap()->end() ; Calo_itr++){ + + G4double* Info = *(Calo_itr->second); + double Energy = RandGauss::shoot(Info[0],Chio_NS::ResoEnergy); + if(Energy>Chio_NS::EnergyThreshold){ + double Time = RandGauss::shoot(Info[1],Chio_NS::ResoTime); + int DetectorNbr = (int) Info[2]; + //m_Event->SetEnergy(DetectorNbr,Energy); + //m_Event->SetTime(DetectorNbr,Time); + } + }*/ + // clear map for next event + CaloHitMap->clear(); +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +//////////////////////////////////////////////////////////////// +void Chio::InitializeScorers() { + // This check is necessary in case the geometry is reloaded + bool already_exist = false; + m_ChioScorer = CheckScorer("ChioScorer",already_exist) ; + + if(already_exist) + return ; + + // Otherwise the scorer is initialised + vector<int> level; level.push_back(0); + G4VPrimitiveScorer* Calorimeter= new CALORIMETERSCORERS::PS_Calorimeter("Calorimeter",level, 0) ; + //and register it to the multifunctionnal detector + m_ChioScorer->RegisterPrimitive(Calorimeter); + G4SDManager::GetSDMpointer()->AddNewDetector(m_ChioScorer) ; +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +//////////////////////////////////////////////////////////////////////////////// +// Construct Method to be pass to the DetectorFactory // +//////////////////////////////////////////////////////////////////////////////// +NPS::VDetector* Chio::Construct(){ + return (NPS::VDetector*) new Chio(); +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +//////////////////////////////////////////////////////////////////////////////// +// Registering the construct method to the factory // +//////////////////////////////////////////////////////////////////////////////// +extern"C" { + class proxy_nps_Chio{ + public: + proxy_nps_Chio(){ + NPS::DetectorFactory::getInstance()->AddToken("Chio","Chio"); + NPS::DetectorFactory::getInstance()->AddDetector("Chio",Chio::Construct); + } + }; + + proxy_nps_Chio p_nps_Chio; +} diff --git a/NPSimulation/Detectors/Chio/Chio.hh b/NPSimulation/Detectors/Chio/Chio.hh new file mode 100644 index 000000000..65ab3a9cf --- /dev/null +++ b/NPSimulation/Detectors/Chio/Chio.hh @@ -0,0 +1,117 @@ +#ifndef Chio_h +#define Chio_h 1 +/***************************************************************************** + * Copyright (C) 2009-XYEARX 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: Adrien Matta contact address: matta@lpccaen.in2p3.fr * + * * + * Creation Date : December 2016 * + * Last update : * + *---------------------------------------------------------------------------* + * Decription: * + * This class describe Chio simulation * + * * + *---------------------------------------------------------------------------* + * Comment: * + * * + *****************************************************************************/ + +// C++ header +#include <string> +#include <vector> +using namespace std; + +// G4 headers +#include "G4ThreeVector.hh" +#include "G4RotationMatrix.hh" +#include "G4LogicalVolume.hh" +#include "G4MultiFunctionalDetector.hh" + +// NPTool header +#include "NPSVDetector.hh" +#include "TChio_anData.h" +#include "NPInputParser.h" + +class Chio : public NPS::VDetector{ + //////////////////////////////////////////////////// + /////// Default Constructor and Destructor ///////// + //////////////////////////////////////////////////// + public: + Chio() ; + virtual ~Chio() ; + + //////////////////////////////////////////////////// + /////// Specific Function of this Class /////////// + //////////////////////////////////////////////////// + public: + // Cartesian + void AddDetector(G4ThreeVector POS, string Shape); + // Spherical + void AddDetector(double R,double Theta,double Phi,string Shape); + + + G4LogicalVolume* BuildDetector(); + + private: + G4LogicalVolume* m_SquareDetector; + G4LogicalVolume* m_CylindricalDetector; + + //////////////////////////////////////////////////// + ////// Inherite from NPS::VDetector class ///////// + //////////////////////////////////////////////////// + public: + // Read stream at Configfile to pick-up parameters of detector (Position,...) + // Called in DetecorConstruction::ReadDetextorConfiguration Method + void ReadConfiguration(NPL::InputParser) ; + + // Construct detector and inialise sensitive part. + // Called After DetecorConstruction::AddDetector Method + void ConstructDetector(G4LogicalVolume* world) ; + + // Add Detector branch to the EventTree. + // Called After DetecorConstruction::AddDetector Method + void InitializeRootOutput() ; + + // Read sensitive part and fill the Root tree. + // Called at in the EventAction::EndOfEventAvtion + void ReadSensitive(const G4Event* event) ; + + public: // Scorer + // Initialize all Scorer used by the MUST2Array + void InitializeScorers() ; + + // Associated Scorer + G4MultiFunctionalDetector* m_ChioScorer ; + //////////////////////////////////////////////////// + ///////////Event class to store Data//////////////// + //////////////////////////////////////////////////// + private: + TChio_anData* m_Event ; + + //////////////////////////////////////////////////// + ///////////////Private intern Data////////////////// + //////////////////////////////////////////////////// + private: // Geometry + // Detector Coordinate + vector<double> m_R; + vector<double> m_Theta; + vector<double> m_Phi; + + // Shape type + vector<string> m_Shape ; + + // Visualisation Attribute + G4VisAttributes* m_VisChamber; + G4VisAttributes* m_VisWindows; + G4VisAttributes* m_VisGas; + + // Needed for dynamic loading of the library + public: + static NPS::VDetector* Construct(); +}; +#endif diff --git a/NPSimulation/Process/CMakeLists.txt b/NPSimulation/Process/CMakeLists.txt index 7235036b3..2c6eeb037 100644 --- a/NPSimulation/Process/CMakeLists.txt +++ b/NPSimulation/Process/CMakeLists.txt @@ -1,2 +1,2 @@ -add_library(NPSProcess SHARED FastDriftElectron.cc NPIonIonInelasticPhysic.cc PhysicsList.cc G4DriftElectron.cc ) +add_library(NPSProcess SHARED FastDriftElectron.cc NPIonIonInelasticPhysic.cc PhysicsList.cc G4DriftElectron.cc G4IonizationWithDE.cc G4DriftElectronPhysics.cc G4DEAbsorption.cc G4DEAmplification.cc G4DETransport.cc) target_link_libraries(NPSProcess ${ROOT_LIBRARIES} ${Geant4_LIBRARIES} ${NPLib_LIBRARIES} -lNPInitialConditions -lNPInteractionCoordinates) diff --git a/NPSimulation/Process/G4DEAbsorption.cc b/NPSimulation/Process/G4DEAbsorption.cc new file mode 100644 index 000000000..1fa070121 --- /dev/null +++ b/NPSimulation/Process/G4DEAbsorption.cc @@ -0,0 +1,133 @@ +// +// ******************************************************************** +// * License and Disclaimer * +// * * +// * The Geant4 software is copyright of the Copyright Holders of * +// * the Geant4 Collaboration. It is provided under the terms and * +// * conditions of the Geant4 Software License, included in the file * +// * LICENSE and available at http://cern.ch/geant4/license . These * +// * include a list of copyright holders. * +// * * +// * Neither the authors of this software system, nor their employing * +// * institutes,nor the agencies providing financial support for this * +// * work make any representation or warranty, express or implied, * +// * regarding this software system or assume any liability for its * +// * use. Please see the license in the file LICENSE and URL above * +// * for the full disclaimer and the limitation of liability. * +// * * +// * This code implementation is the result of the scientific and * +// * technical work of the GEANT4 collaboration. * +// * By using, copying, modifying or distributing the software (or * +// * any work based on the software) you agree to acknowledge its * +// * use in resulting scientific publications, and indicate your * +// * acceptance of all terms of the Geant4 Software license. * +// ******************************************************************** +// +// +// $Id: G4DEAbsorption.hh 69576 2016-12-19 13:48:13Z gcosmo $ +// +//////////////////////////////////////////////////////////////////////// +// Optical DriftElectron Absorption Implementation +//////////////////////////////////////////////////////////////////////// +// +// File: G4DEAbsorption.hh +// Description: Discrete Process -- Bulk absorption of Drift Electron +// Version: 1.0 +// Created: 2016-12-19 +// Author: Adrien Matta +// Updated: +// +// +// +// +// mail: matta@lpccaen.in2p3.fr +// +//////////////////////////////////////////////////////////////////////// + +#include "G4ios.hh" +#include "G4DEProcessSubType.hh" + +#include "G4DEAbsorption.hh" + +///////////////////////// +// Class Implementation +///////////////////////// + +////////////// +// Operators +////////////// + +// G4DEAbsorption::operator=(const G4DEAbsorption &right) +// { +// } + +///////////////// +// Constructors +///////////////// + +G4DEAbsorption::G4DEAbsorption(const G4String& processName, G4ProcessType type) + : G4VDiscreteProcess(processName, type) +{ + if (verboseLevel>0) { + G4cout << GetProcessName() << " is created " << G4endl; + } + + SetProcessSubType(fDEAbsorption); +} + +// G4DEAbsorption::G4DEAbsorption(const G4OpAbsorpton &right) +// { +// } + +//////////////// +// Destructors +//////////////// + +G4DEAbsorption::~G4DEAbsorption(){} + +//////////// +// Methods +//////////// + +// PostStepDoIt +// ------------- +// + G4VParticleChange* +G4DEAbsorption::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep) +{ + aParticleChange.Initialize(aTrack); + + const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle(); + G4double theDriftElectronMomentum = aParticle->GetTotalMomentum(); + + aParticleChange.ProposeLocalEnergyDeposit(theDriftElectronMomentum); + aParticleChange.ProposeTrackStatus(fStopAndKill); + if (verboseLevel>0) { + G4cout << "\n** Drift electron absorbed! **" << G4endl; + } + return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep); +} + + +// GetMeanFreePath +// --------------- +// +G4double G4DEAbsorption::GetMeanFreePath(const G4Track& aTrack, + G4double , + G4ForceCondition* ) +{ + const G4Material* aMaterial = aTrack.GetMaterial(); + + G4MaterialPropertiesTable* aMaterialPropertyTable; + + // If none provided, the absorbtion is maximum + G4double AttenuationLength = 0; + + aMaterialPropertyTable = aMaterial->GetMaterialPropertiesTable(); + + if ( aMaterialPropertyTable ){ + if( aMaterialPropertyTable->ConstPropertyExists("DE_ABSLENGTH")) + AttenuationLength = aMaterialPropertyTable->GetConstProperty("DE_ABSLENGTH"); + } + return AttenuationLength; +} diff --git a/NPSimulation/Process/G4DEAbsorption.hh b/NPSimulation/Process/G4DEAbsorption.hh new file mode 100644 index 000000000..4e62552d3 --- /dev/null +++ b/NPSimulation/Process/G4DEAbsorption.hh @@ -0,0 +1,127 @@ +// +// ******************************************************************** +// * License and Disclaimer * +// * * +// * The Geant4 software is copyright of the Copyright Holders of * +// * the Geant4 Collaboration. It is provided under the terms and * +// * conditions of the Geant4 Software License, included in the file * +// * LICENSE and available at http://cern.ch/geant4/license . These * +// * include a list of copyright holders. * +// * * +// * Neither the authors of this software system, nor their employing * +// * institutes,nor the agencies providing financial support for this * +// * work make any representation or warranty, express or implied, * +// * regarding this software system or assume any liability for its * +// * use. Please see the license in the file LICENSE and URL above * +// * for the full disclaimer and the limitation of liability. * +// * * +// * This code implementation is the result of the scientific and * +// * technical work of the GEANT4 collaboration. * +// * By using, copying, modifying or distributing the software (or * +// * any work based on the software) you agree to acknowledge its * +// * use in resulting scientific publications, and indicate your * +// * acceptance of all terms of the Geant4 Software license. * +// ******************************************************************** +// +// +// $Id: G4DEAbsorption.hh 69576 2013-05-08 13:48:13Z gcosmo $ +// +//////////////////////////////////////////////////////////////////////// +// Optical Photon Absorption Class Definition +//////////////////////////////////////////////////////////////////////// +// +// File: G4DEAbsorption.hh +// Description: Discrete Process -- Bulk absorption of Drift Electron +// Version: 1.0 +// Created: 2016-12-19 +// Author: Adrien Matta +// Updated: +// +// +// +// +// mail: matta@lpccaen.in2p3.fr +// +//////////////////////////////////////////////////////////////////////// + +#ifndef G4DEAbsorption_h +#define G4DEAbsorption_h 1 + +///////////// +// Includes +///////////// + +#include "globals.hh" +#include "templates.hh" +#include "Randomize.hh" +#include "G4Step.hh" +#include "G4VDiscreteProcess.hh" +#include "G4DynamicParticle.hh" +#include "G4Material.hh" +#include "G4DriftElectron.hh" + +// Class Description: +// Discrete Process -- Bulk absorption of Drift Electron. +// Class inherits publicly from G4VDiscreteProcess +// Class Description - End: + +///////////////////// +// Class Definition +///////////////////// + +class G4DEAbsorption : public G4VDiscreteProcess +{ + +public: + + //////////////////////////////// + // Constructors and Destructor + //////////////////////////////// + + G4DEAbsorption(const G4String& processName = "DEAbsorption", + G4ProcessType type = fTransportation); + ~G4DEAbsorption(); + +private: + + G4DEAbsorption(const G4DEAbsorption &right); + + ////////////// + // Operators + ////////////// + + G4DEAbsorption& operator=(const G4DEAbsorption &right); + +public: + + //////////// + // Methods + //////////// + + G4bool IsApplicable(const G4ParticleDefinition& aParticleType); + // Returns true -> 'is applicable' only for an drift electron. + + G4double GetMeanFreePath(const G4Track& aTrack, + G4double , + G4ForceCondition* ); + // Returns the absorption length for bulk absorption of drift + // electron in media with a specified attenuation length. + + G4VParticleChange* PostStepDoIt(const G4Track& aTrack, + const G4Step& aStep); + // This is the method implementing bulk absorption of drift + // electron. + +}; + +//////////////////// +// Inline methods +//////////////////// + +inline +G4bool G4DEAbsorption::IsApplicable(const G4ParticleDefinition& aParticleType) +{ + return ( &aParticleType == G4DriftElectron::DriftElectron() ); +} + +#endif /* G4DEAbsorption_h */ diff --git a/NPSimulation/Process/G4DEAmplification.cc b/NPSimulation/Process/G4DEAmplification.cc new file mode 100644 index 000000000..42f9a55d9 --- /dev/null +++ b/NPSimulation/Process/G4DEAmplification.cc @@ -0,0 +1,200 @@ +// +// ******************************************************************** +// * License and Disclaimer * +// * * +// * The Geant4 software is copyright of the Copyright Holders of * +// * the Geant4 Collaboration. It is provided under the terms and * +// * conditions of the Geant4 Software License, included in the file * +// * LICENSE and available at http://cern.ch/geant4/license . These * +// * include a list of copyright holders. * +// * * +// * Neither the authors of this software system, nor their employing * +// * institutes,nor the agencies providing financial support for this * +// * work make any representation or warranty, express or implied, * +// * regarding this software system or assume any liability for its * +// * use. Please see the license in the file LICENSE and URL above * +// * for the full disclaimer and the limitation of liability. * +// * * +// * This code implementation is the result of the scientific and * +// * technical work of the GEANT4 collaboration. * +// * By using, copying, modifying or distributing the software (or * +// * any work based on the software) you agree to acknowledge its * +// * use in resulting scientific publications, and indicate your * +// * acceptance of all terms of the Geant4 Software license. * +// ******************************************************************** +// +// +// $Id: G4DEAmplification.hh 69576 2016-12-19 13:48:13Z gcosmo $ +// +//////////////////////////////////////////////////////////////////////// +// Optical DriftElectron Amplification Implementation +//////////////////////////////////////////////////////////////////////// +// +// File: G4DEAmplification.hh +// Description: Discrete Process -- Bulk amplification of Drift Electron +// Version: 1.0 +// Created: 2016-12-19 +// Author: Adrien Matta +// Updated: +// +// +// +// +// mail: matta@lpccaen.in2p3.fr +// +//////////////////////////////////////////////////////////////////////// + +#include "G4ios.hh" +#include "G4DEProcessSubType.hh" + +#include "G4DEAmplification.hh" +#include "G4PhysicalConstants.hh" +#include <iostream> +using namespace std; +///////////////////////// +// Class Implementation +///////////////////////// + +////////////// +// Operators +////////////// + +// G4DEAmplification::operator=(const G4DEAmplification &right) +// { +// } + +///////////////// +// Constructors +///////////////// + +G4DEAmplification::G4DEAmplification(const G4String& processName, G4ProcessType type) + : G4VRestDiscreteProcess(processName, type) +{ + if (verboseLevel>0) { + G4cout << GetProcessName() << " is created " << G4endl; + } + + SetProcessSubType(fDEAmplification); +} + +// G4DEAmplification::G4DEAmplification(const G4OpAbsorpton &right) +// { +// } + +//////////////// +// Destructors +//////////////// + +G4DEAmplification::~G4DEAmplification(){} + +//////////// +// Methods +//////////// + +// PostStepDoIt +// ------------- +// + G4VParticleChange* +G4DEAmplification::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep) +{ + // Get the primary track + aParticleChange.Initialize(aTrack); + + // Check that the parent process is not amplification + if(aTrack.GetCreatorProcess()->GetProcessName()=="DEAmplification" ) + return G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep); + + const G4Material* aMaterial = aTrack.GetMaterial(); + + G4StepPoint* pPreStepPoint = aStep.GetPreStepPoint(); + + + + G4ThreeVector x0 = pPreStepPoint->GetPosition(); + G4ThreeVector p0 = aStep.GetDeltaPosition().unit(); + G4double t0 = pPreStepPoint->GetGlobalTime(); + + // Get the material table + G4MaterialPropertiesTable* aMaterialPropertiesTable = + aMaterial->GetMaterialPropertiesTable(); + if (!aMaterialPropertiesTable){ // Does the table exist + return G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep); + } + + if(!aMaterialPropertiesTable->ConstPropertyExists("DE_PAIRENERGY") || + !aMaterialPropertiesTable->ConstPropertyExists("DE_AMPLIFICATION") || + !aMaterialPropertiesTable->ConstPropertyExists("DE_YIELD") || + !aMaterialPropertiesTable->ConstPropertyExists("DE_DRIFTSPEED") || + !aMaterialPropertiesTable->ConstPropertyExists("DE_TRANSVERSALSPREAD") || + !aMaterialPropertiesTable->ConstPropertyExists("DE_LONGITUDINALSPREAD")) + return G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep); + + G4double pair_energy = 0; + pair_energy = + aMaterialPropertiesTable->GetConstProperty("DE_PAIRENERGY"); + + G4double amplification=0; + amplification= + aMaterialPropertiesTable->GetConstProperty("DE_AMPLIFICATION"); + G4double Yield = 0; + Yield= + aMaterialPropertiesTable->GetConstProperty("DE_YIELD"); + + G4int number_electron = amplification*Yield; + //if no electron leave + if(number_electron<1){ + aParticleChange.SetNumberOfSecondaries(0); + return G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep); + } + + aParticleChange.SetNumberOfSecondaries(number_electron); + + // Create the secondary tracks + for(G4int i = 0 ; i < number_electron ; i++){ + // Electron follow the field direction + G4ThreeVector field(0,1,0); + G4double v_drift = aMaterialPropertiesTable->GetConstProperty("DE_DRIFTSPEED"); + G4double v_long = G4RandGauss::shoot(0,aMaterialPropertiesTable->GetConstProperty("DE_LONGITUDINALSPREAD")/v_drift); + G4double v_trans = G4RandGauss::shoot(0,aMaterialPropertiesTable->GetConstProperty("DE_TRANSVERSALSPREAD")/v_drift); +cout << v_drift << " " << v_long << " " << v_trans << endl; + G4ThreeVector v = v_drift*field+v_long*G4ThreeVector(0,0,1)+v_trans*G4ThreeVector(1,0,0); + + // Random Position along the step with matching time + G4double rand = G4UniformRand(); + G4ThreeVector pos = x0 + rand * aStep.GetDeltaPosition(); + G4double time = t0+ rand* aStep.GetDeltaTime(); + + G4DynamicParticle* particle = new G4DynamicParticle(G4DriftElectron::DriftElectron(),v.unit(), electron_mass_c2*v.mag()/c_squared); + G4Track* aSecondaryTrack = new G4Track(particle,time,pos); + + aSecondaryTrack->SetTouchableHandle( + aStep.GetPreStepPoint()->GetTouchableHandle()); + + aSecondaryTrack->SetParentID(aTrack.GetTrackID()); + aSecondaryTrack->SetTouchableHandle(aStep.GetPreStepPoint()->GetTouchableHandle()); + aParticleChange.AddSecondary(aSecondaryTrack); + } + + return G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep); +} + + +// GetMeanFreePath +// --------------- +// +G4double G4DEAmplification::GetMeanFreePath(const G4Track& , + G4double , + G4ForceCondition* condition) +{ + *condition = StronglyForced; + return DBL_MAX; +} +// GetMeanLifeTime +// --------------- +// +G4double G4DEAmplification::GetMeanLifeTime(const G4Track& , + G4ForceCondition* condition) +{ + *condition = StronglyForced; + return DBL_MAX; +} diff --git a/NPSimulation/Process/G4DEAmplification.hh b/NPSimulation/Process/G4DEAmplification.hh new file mode 100644 index 000000000..4d31bac16 --- /dev/null +++ b/NPSimulation/Process/G4DEAmplification.hh @@ -0,0 +1,130 @@ +// +// ******************************************************************** +// * License and Disclaimer * +// * * +// * The Geant4 software is copyright of the Copyright Holders of * +// * the Geant4 Collaboration. It is provided under the terms and * +// * conditions of the Geant4 Software License, included in the file * +// * LICENSE and available at http://cern.ch/geant4/license . These * +// * include a list of copyright holders. * +// * * +// * Neither the authors of this software system, nor their employing * +// * institutes,nor the agencies providing financial support for this * +// * work make any representation or warranty, express or implied, * +// * regarding this software system or assume any liability for its * +// * use. Please see the license in the file LICENSE and URL above * +// * for the full disclaimer and the limitation of liability. * +// * * +// * This code implementation is the result of the scientific and * +// * technical work of the GEANT4 collaboration. * +// * By using, copying, modifying or distributing the software (or * +// * any work based on the software) you agree to acknowledge its * +// * use in resulting scientific publications, and indicate your * +// * acceptance of all terms of the Geant4 Software license. * +// ******************************************************************** +// +// +// $Id: G4DEAmplification.hh 69576 2013-05-08 13:48:13Z gcosmo $ +// +//////////////////////////////////////////////////////////////////////// +// Optical Photon Amplification Class Definition +//////////////////////////////////////////////////////////////////////// +// +// File: G4DEAmplification.hh +// Description: Discrete Process -- Bulk absorption of Drift Electron +// Version: 1.0 +// Created: 2016-12-19 +// Author: Adrien Matta +// Updated: +// +// +// +// +// mail: matta@lpccaen.in2p3.fr +// +//////////////////////////////////////////////////////////////////////// + +#ifndef G4DEAmplification_h +#define G4DEAmplification_h 1 + +///////////// +// Includes +///////////// + +#include "globals.hh" +#include "templates.hh" +#include "Randomize.hh" +#include "G4Step.hh" +#include "G4VRestDiscreteProcess.hh" +#include "G4DynamicParticle.hh" +#include "G4Material.hh" +#include "G4DriftElectron.hh" + +// Class Description: +// Discrete Process -- Bulk amplification of Drift Electron. +// Class inherits publicly from G4VDiscreteProcess +// Class Description - End: + +///////////////////// +// Class Definition +///////////////////// + +class G4DEAmplification : public G4VRestDiscreteProcess +{ + +public: + + //////////////////////////////// + // Constructors and Destructor + //////////////////////////////// + + G4DEAmplification(const G4String& processName = "DEAmplification", + G4ProcessType type = fTransportation); + ~G4DEAmplification(); + +private: + + G4DEAmplification(const G4DEAmplification &right); + + ////////////// + // Operators + ////////////// + + G4DEAmplification& operator=(const G4DEAmplification &right); + +public: + + //////////// + // Methods + //////////// + + G4bool IsApplicable(const G4ParticleDefinition& aParticleType); + // Returns true -> 'is applicable' only for an drift electron. + + G4double GetMeanFreePath(const G4Track& aTrack, + G4double , + G4ForceCondition* ); + // Returns true and strongly enforced + + G4double GetMeanLifeTime(const G4Track& aTrack, + G4ForceCondition* ); + // Returns true and strongly enforced + + G4VParticleChange* PostStepDoIt(const G4Track& aTrack, + const G4Step& aStep); + // This is the method implementing bulk amplification of drift + // electron. + +}; + +//////////////////// +// Inline methods +//////////////////// + +inline +G4bool G4DEAmplification::IsApplicable(const G4ParticleDefinition& aParticleType) +{ + return ( &aParticleType == G4DriftElectron::DriftElectron() ); +} + +#endif /* G4DEAmplification_h */ diff --git a/NPSimulation/Process/G4DEProcessSubType.hh b/NPSimulation/Process/G4DEProcessSubType.hh new file mode 100644 index 000000000..584896708 --- /dev/null +++ b/NPSimulation/Process/G4DEProcessSubType.hh @@ -0,0 +1,52 @@ +// +// ******************************************************************** +// * License and Disclaimer * +// * * +// * The Geant4 software is copyright of the Copyright Holders of * +// * the Geant4 Collaboration. It is provided under the terms and * +// * conditions of the Geant4 Software License, included in the file * +// * LICENSE and available at http://cern.ch/geant4/license . These * +// * include a list of copyright holders. * +// * * +// * Neither the authors of this software system, nor their employing * +// * institutes,nor the agencies providing financial support for this * +// * work make any representation or warranty, express or implied, * +// * regarding this software system or assume any liability for its * +// * use. Please see the license in the file LICENSE and URL above * +// * for the full disclaimer and the limitation of liability. * +// * * +// * This code implementation is the result of the scientific and * +// * technical work of the GEANT4 collaboration. * +// * By using, copying, modifying or distributing the software (or * +// * any work based on the software) you agree to acknowledge its * +// * use in resulting scientific publications, and indicate your * +// * acceptance of all terms of the Geant4 Software license. * +// ******************************************************************** +// +// +// $Id: G4OpProcessSubType.hh 69576 2013-05-08 13:48:13Z gcosmo $ +// +//--------------------------------------------------------------- +// +// G4OpProcessSubType.hh +// +// Class Description: +// This is an enumerator to define sub-type of optical processes +// +// Creation date: 23.09.2008 +// Modifications: +// +//--------------------------------------------------------------- + +#ifndef G4DEProcessSubType_h +#define G4DEProcessSubType_h 1 + +enum G4DEProcessSubType +{ + fDEAbsorption = 31, + fDETransport= 32, + fDEAmplification = 33 + +}; + +#endif diff --git a/NPSimulation/Process/G4DETransport.cc b/NPSimulation/Process/G4DETransport.cc new file mode 100644 index 000000000..87578f5d8 --- /dev/null +++ b/NPSimulation/Process/G4DETransport.cc @@ -0,0 +1,165 @@ +// +// ******************************************************************** +// * License and Disclaimer * +// * * +// * The Geant4 software is copyright of the Copyright Holders of * +// * the Geant4 Collaboration. It is provided under the terms and * +// * conditions of the Geant4 Software License, included in the file * +// * LICENSE and available at http://cern.ch/geant4/license . These * +// * include a list of copyright holders. * +// * * +// * Neither the authors of this software system, nor their employing * +// * institutes,nor the agencies providing financial support for this * +// * work make any representation or warranty, express or implied, * +// * regarding this software system or assume any liability for its * +// * use. Please see the license in the file LICENSE and URL above * +// * for the full disclaimer and the limitation of liability. * +// * * +// * This code implementation is the result of the scientific and * +// * technical work of the GEANT4 collaboration. * +// * By using, copying, modifying or distributing the software (or * +// * any work based on the software) you agree to acknowledge its * +// * use in resulting scientific publications, and indicate your * +// * acceptance of all terms of the Geant4 Software license. * +// ******************************************************************** +// +// +// $Id: G4DETransport.hh 69576 2016-12-19 13:48:13Z gcosmo $ +// +//////////////////////////////////////////////////////////////////////// +// Optical DriftElectron Transport Implementation +//////////////////////////////////////////////////////////////////////// +// +// File: G4DETransport.cc +// Description: Discrete Process -- Bulk transport of Drift Electron +// Version: 1.0 +// Created: 2016-12-19 +// Author: Adrien Matta +// Updated: +// +// +// +// +// mail: matta@lpccaen.in2p3.fr +// +//////////////////////////////////////////////////////////////////////// + +#include "G4ios.hh" +#include "G4DEProcessSubType.hh" +#include "G4PhysicalConstants.hh" +#include "G4SystemOfUnits.hh" +#include "G4DETransport.hh" + +///////////////////////// +// Class Implementation +///////////////////////// + +////////////// +// Operators +////////////// + +// G4DETransport::operator=(const G4DETransport &right) +// { +// } + +///////////////// +// Constructors +///////////////// + +G4DETransport::G4DETransport(const G4String& processName, G4ProcessType type) + : G4VDiscreteProcess(processName, type) +{ + if (verboseLevel>0) { + G4cout << GetProcessName() << " is created " << G4endl; + } + + SetProcessSubType(fDETransport); +} + +// G4DETransport::G4DETransport(const G4OpAbsorpton &right) +// { +// } + +//////////////// +// Destructors +//////////////// + +G4DETransport::~G4DETransport(){} + +//////////// +// Methods +//////////// + +// PostStepDoIt +// ------------- +// +G4VParticleChange* +G4DETransport::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep) +{ + aParticleChange.Initialize(aTrack); + G4StepPoint* pPreStepPoint = aStep.GetPreStepPoint(); + + G4ThreeVector x0 = pPreStepPoint->GetPosition(); + G4ThreeVector p0 = aStep.GetDeltaPosition().unit(); + G4double t0 = pPreStepPoint->GetGlobalTime(); + + G4double Energy = pPreStepPoint->GetKineticEnergy(); + + // Get the material table + const G4Material* aMaterial = aTrack.GetMaterial(); + + G4MaterialPropertiesTable* aMaterialPropertiesTable = + aMaterial->GetMaterialPropertiesTable(); + if (!aMaterialPropertiesTable) + return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep); + + + if(!aMaterialPropertiesTable->ConstPropertyExists("DE_PAIRENERGY") || + !aMaterialPropertiesTable->ConstPropertyExists("DE_YIELD") || + !aMaterialPropertiesTable->ConstPropertyExists("DE_TRANSVERSALSPREAD") || + !aMaterialPropertiesTable->ConstPropertyExists("DE_LONGITUDINALSPREAD") || + !aMaterialPropertiesTable->ConstPropertyExists("DE_DRIFTSPEED") ) + return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep); + + + // Electron follow the field direction + G4ThreeVector driftDir(0,1,0); + G4double v_drift = aMaterialPropertiesTable->GetConstProperty("DE_DRIFTSPEED"); + // The time scale is imposed by the drift speed + G4double step = aStep.GetDeltaTime(); + if(!step) + return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep); + + // Should be equal to delta pos y + G4double step_length = step*v_drift; + G4double sigmaTrans = sqrt(2*step_length*aMaterialPropertiesTable->GetConstProperty("DE_TRANSVERSALSPREAD")/v_drift); + G4double sigmaLong = sqrt(2*step_length*aMaterialPropertiesTable->GetConstProperty("DE_LONGITUDINALSPREAD")/v_drift); + + G4double d_long = G4RandGauss::shoot(0,sigmaLong); + G4double d_trans = G4RandGauss::shoot(0,sigmaTrans); + + G4ThreeVector trans(G4RandGauss::shoot(0,d_trans),0,0); + trans.rotateY(twopi*G4UniformRand()); + G4ThreeVector d_Pos = (step_length+d_long)*driftDir+trans; + // Random Position along the step with matching time + G4ThreeVector pos = x0 + d_Pos; + G4double time = t0 + step; + + aParticleChange.ProposeMomentumDirection(d_Pos.unit()); + aParticleChange.ProposeEnergy(Energy); + aParticleChange.ProposePosition(pos); + aParticleChange.ProposeGlobalTime(time); + + return &aParticleChange; +} + + +// GetMeanFreePath +// --------------- +// +G4double G4DETransport::GetMeanFreePath(const G4Track& , + G4double , + G4ForceCondition* ) +{ + return 0.5*mm; +} diff --git a/NPSimulation/Process/G4DETransport.hh b/NPSimulation/Process/G4DETransport.hh new file mode 100644 index 000000000..6505959f1 --- /dev/null +++ b/NPSimulation/Process/G4DETransport.hh @@ -0,0 +1,127 @@ +// +// ******************************************************************** +// * License and Disclaimer * +// * * +// * The Geant4 software is copyright of the Copyright Holders of * +// * the Geant4 Collaboration. It is provided under the terms and * +// * conditions of the Geant4 Software License, included in the file * +// * LICENSE and available at http://cern.ch/geant4/license . These * +// * include a list of copyright holders. * +// * * +// * Neither the authors of this software system, nor their employing * +// * institutes,nor the agencies providing financial support for this * +// * work make any representation or warranty, express or implied, * +// * regarding this software system or assume any liability for its * +// * use. Please see the license in the file LICENSE and URL above * +// * for the full disclaimer and the limitation of liability. * +// * * +// * This code implementation is the result of the scientific and * +// * technical work of the GEANT4 collaboration. * +// * By using, copying, modifying or distributing the software (or * +// * any work based on the software) you agree to acknowledge its * +// * use in resulting scientific publications, and indicate your * +// * acceptance of all terms of the Geant4 Software license. * +// ******************************************************************** +// +// +// $Id: G4DETransport.hh 69576 2013-05-08 13:48:13Z gcosmo $ +// +//////////////////////////////////////////////////////////////////////// +// Optical Photon Transport Class Definition +//////////////////////////////////////////////////////////////////////// +// +// File: G4DETransport.hh +// Description: Discrete Process -- Bulk transport of Drift Electron +// Version: 1.0 +// Created: 2016-12-19 +// Author: Adrien Matta +// Updated: +// +// +// +// +// mail: matta@lpccaen.in2p3.fr +// +//////////////////////////////////////////////////////////////////////// + +#ifndef G4DETransport_h +#define G4DETransport_h 1 + +///////////// +// Includes +///////////// + +#include "globals.hh" +#include "templates.hh" +#include "Randomize.hh" +#include "G4Step.hh" +#include "G4VDiscreteProcess.hh" +#include "G4DynamicParticle.hh" +#include "G4Material.hh" +#include "G4DriftElectron.hh" + +// Class Description: +// Discrete Process -- Bulk transport of Drift Electron. +// Class inherits publicly from G4VDiscreteProcess +// Class Description - End: + +///////////////////// +// Class Definition +///////////////////// + +class G4DETransport : public G4VDiscreteProcess +{ + +public: + + //////////////////////////////// + // Constructors and Destructor + //////////////////////////////// + + G4DETransport(const G4String& processName = "DETransport", + G4ProcessType type = fTransportation); + ~G4DETransport(); + +private: + + G4DETransport(const G4DETransport &right); + + ////////////// + // Operators + ////////////// + + G4DETransport& operator=(const G4DETransport &right); + +public: + + //////////// + // Methods + //////////// + + G4bool IsApplicable(const G4ParticleDefinition& aParticleType); + // Returns true -> 'is applicable' only for an drift electron. + + G4double GetMeanFreePath(const G4Track& aTrack, + G4double , + G4ForceCondition* ); + // Returns the transport length for bulk transport of drift + // electron in media with a specified attenuation length. + + G4VParticleChange* PostStepDoIt(const G4Track& aTrack, + const G4Step& aStep); + // This is the method implementing bulk transport of drift + // electron. + +}; + +//////////////////// +// Inline methods +//////////////////// + +inline +G4bool G4DETransport::IsApplicable(const G4ParticleDefinition& aParticleType) +{ + return ( &aParticleType == G4DriftElectron::DriftElectron() ); +} + +#endif /* G4DETransport_h */ diff --git a/NPSimulation/Process/G4DriftElectron.cc b/NPSimulation/Process/G4DriftElectron.cc index 486cfe8c4..04a09993a 100644 --- a/NPSimulation/Process/G4DriftElectron.cc +++ b/NPSimulation/Process/G4DriftElectron.cc @@ -24,7 +24,7 @@ // ******************************************************************** // // -// $Id: G4DriftElectron.cc 67971 2013-03-13 10:13:24Z gcosmo $ +// $Id: G4DriftElectron.cc 67971 2016-12-19 10:00:30Z gcosmo $ // // // ---------------------------------------------------------------------- @@ -41,7 +41,7 @@ #include "G4ParticleTable.hh" // ###################################################################### -// ### OPTICAL PHOTON ### +// ### Drift Electron ### // ###################################################################### G4DriftElectron* G4DriftElectron::theInstance = 0; @@ -55,23 +55,23 @@ G4DriftElectron* G4DriftElectron::Definition() G4ParticleDefinition* anInstance = pTable->FindParticle(name); if (anInstance ==0) { - // create particle - // - // Arguments for constructor are as follows - // name mass width charge - // 2*spin parity C-conjugation - // 2*Isospin 2*Isospin3 G-parity - // type lepton number baryon number PDG encoding - // stable lifetime decay table - // shortlived subType anti_encoding + // create particle + // + // Arguments for constructor are as follows + // name mass width charge + // 2*spin parity C-conjugation + // 2*Isospin 2*Isospin3 G-parity + // type lepton number baryon number PDG encoding + // stable lifetime decay table + // shortlived subType anti_encoding anInstance = new G4ParticleDefinition( - name, electron_mass_c2, 0.0*MeV, -1.*eplus, - 1, 0, 0, - 0, 0, 0, - "lepton", 1, 0, 11, - true, -1.0, NULL, - false, "e" - ); + name, electron_mass_c2, 0.0*MeV, -1.*eplus, + 1, 0, 0, + 0, 0, 0, + "lepton", 1, 0, 11, + true, -1.0, NULL, + false, "e" + ); } theInstance = reinterpret_cast<G4DriftElectron*>(anInstance); return theInstance; diff --git a/NPSimulation/Process/G4DriftElectronPhysics.cc b/NPSimulation/Process/G4DriftElectronPhysics.cc new file mode 100644 index 000000000..1c9f6e954 --- /dev/null +++ b/NPSimulation/Process/G4DriftElectronPhysics.cc @@ -0,0 +1,253 @@ +// +// ******************************************************************** +// * License and Disclaimer * +// * * +// * The Geant4 software is copyright of the Copyright Holders of * +// * the Geant4 Collaboration. It is provided under the terms and * +// * conditions of the Geant4 Software License, included in the file * +// * LICENSE and available at http://cern.ch/geant4/license . These * +// * include a list of copyright holders. * +// * * +// * Neither the authors of this software system, nor their employing * +// * institutes,nor the agencies providing financial support for this * +// * work make any representation or warranty, express or implied, * +// * regarding this software system or assume any liability for its * +// * use. Please see the license in the file LICENSE and URL above * +// * for the full disclaimer and the limitation of liability. * +// * * +// * This code implementation is the result of the scientific and * +// * technical work of the GEANT4 collaboration. * +// * By using, copying, modifying or distributing the software (or * +// * any work based on the software) you agree to acknowledge its * +// * use in resulting scientific publications, and indicate your * +// * acceptance of all terms of the Geant4 Software license. * +// ******************************************************************** +// +// +//--------------------------------------------------------------------------- +// +// ClassName: G4DriftElectronPhysics +// +// Author: Adrie Matta 2016-12-19 +// +// Modified: +// +// +//---------------------------------------------------------------------------- +// + +#include "G4DriftElectronPhysics.hh" + +#include "G4DEAbsorption.hh" +#include "G4IonizationWithDE.hh" +#include "G4DEAmplification.hh" +#include "G4DETransport.hh" +#include "G4LossTableManager.hh" +#include "G4EmSaturation.hh" +#include<iostream> +// factory +#include "G4PhysicsConstructorFactory.hh" +// +G4_DECLARE_PHYSCONSTR_FACTORY(G4DriftElectronPhysics); + + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +G4DriftElectronPhysics::G4DriftElectronPhysics(G4int verbose, const G4String& name) + : G4VPhysicsConstructor(name), + + fMaxNumDriftElectron(100), + fMaxBetaChange(10.0), + fYieldFactor(1.) +{ + verboseLevel = verbose; + + for ( G4int i=0; i<kDENoProcess; i++ ) { + fProcessUse.push_back(true); + fProcessTrackSecondariesFirst.push_back(true); + } +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +G4DriftElectronPhysics::~G4DriftElectronPhysics() +{ +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +void G4DriftElectronPhysics::PrintStatistics() const +{ +// Print all processes activation and their parameters + + for ( G4int i=0; i<kDENoProcess; i++ ) { + G4cout << " " << G4DriftElectronProcessName(i) << " process: "; + if ( ! fProcessUse[i] ) { + G4cout << "not used" << G4endl; + } + else { + G4cout << "used" << G4endl; + if ( i == kDEIonization ) { + if ( fProcessTrackSecondariesFirst[kDEIonization] ) { + G4cout << " Track secondaries first: activated" << G4endl; + } + else { + G4cout << " Track secondaries first: inactivated" << G4endl; + } + } + if ( i == kDEAmplification ) { + if ( fProcessTrackSecondariesFirst[kDEAmplification] ) { + G4cout << " Track secondaries first: activated" << G4endl; + } + else { + G4cout << " Track secondaries first: inactivated" << G4endl; + } + } + if ( i == kDETransport) { + if ( fProcessTrackSecondariesFirst[kDETransport] ) { + G4cout << " Track secondaries first: activated" << G4endl; + } + else { + G4cout << " Track secondaries first: inactivated" << G4endl; + } + } + if ( i == kDEAbsorption ) { + if ( fProcessTrackSecondariesFirst[kDEAbsorption] ) { + G4cout << " Track secondaries first: activated" << G4endl; + } + else { + G4cout << " Track secondaries first: inactivated" << G4endl; + } + } + } + } +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +#include "G4DriftElectron.hh" + +void G4DriftElectronPhysics::ConstructParticle() +{ +/// Instantiate particles. + + // Drift Electron + G4DriftElectron::DriftElectronDefinition(); +} + + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +#include "G4Threading.hh" +#include "G4ParticleDefinition.hh" +#include "G4ProcessManager.hh" +#include "G4AutoDelete.hh" +#include<iostream> +using namespace std; +void G4DriftElectronPhysics::ConstructProcess() +{ +// Construct drift electron processes. + + if(verboseLevel>0) + G4cout <<"G4DriftElectronPhysics:: Add DriftElectron Physics Processes"<< G4endl; + +// A vector of optical processes + std::vector<G4VProcess*> DEProcesses; + + for ( G4int i=0; i<kDENoProcess; i++ ) DEProcesses.push_back(NULL); + + // Add DriftElectron Processes + G4IonizationWithDE* DEIonizationProcess = new G4IonizationWithDE(); + DEProcesses[kDEIonization] = DEIonizationProcess; + + G4DEAbsorption* DEAbsorptionProcess = new G4DEAbsorption(); + DEProcesses[kDEAbsorption] = DEAbsorptionProcess; + + G4DETransport * DETransportProcess = new G4DETransport(); + DEProcesses[kDETransport] = DETransportProcess; + + G4DEAmplification* DEAmplificationProcess = new G4DEAmplification(); + DEProcesses[kDEAmplification] = DEAmplificationProcess; + + G4ProcessManager * pManager = 0; + pManager = G4DriftElectron::DriftElectron()->GetProcessManager(); + + if (!pManager) { + std::ostringstream o; + o << "Drift electron without a Process Manager"; + G4Exception("G4DriftElectronPhysics::ConstructProcess()","", + FatalException,o.str().c_str()); + return; + } + + for ( G4int i=kDEIonization; i<=kDEAbsorption; i++ ) { + if ( fProcessUse[i] ) { + pManager->AddDiscreteProcess(DEProcesses[i]); + } + } + + G4IonizationWithDE* IonizationProcess = new G4IonizationWithDE(); + IonizationProcess->SetTrackSecondariesFirst(fProcessTrackSecondariesFirst[kDEIonization]); + G4EmSaturation* emSaturation = G4LossTableManager::Instance()->EmSaturation(); + IonizationProcess->AddSaturation(emSaturation); + DEProcesses[kDEIonization] = IonizationProcess; + + auto myParticleIterator=GetParticleIterator(); + myParticleIterator->reset(); + + while( (*myParticleIterator)() ){ + + G4ParticleDefinition* particle = myParticleIterator->value(); + G4String particleName = particle->GetParticleName(); + + pManager = particle->GetProcessManager(); + if (!pManager) { + std::ostringstream o; + o << "Particle " << particleName << "without a Process Manager"; + G4Exception("G4DriftElectronPhysics::ConstructProcess()","", + FatalException,o.str().c_str()); + return; // else coverity complains for pManager use below + } + + if( IonizationProcess->IsApplicable(*particle) && + fProcessUse[kDEIonization] ){ + pManager->AddProcess(IonizationProcess); + pManager->SetProcessOrderingToLast(IonizationProcess,idxAtRest); + pManager->SetProcessOrderingToLast(IonizationProcess,idxPostStep); + } + + } + + // Add verbose + for ( G4int i=0; i<kDENoProcess; i++ ) { + if ( fProcessUse[i] ) DEProcesses[i]->SetVerboseLevel(verboseLevel); + } + + if (verboseLevel > 1) PrintStatistics(); + if (verboseLevel > 0) + G4cout << "### " << namePhysics << " physics constructed." << G4endl; +} + +void G4DriftElectronPhysics::SetMaxNumDriftElectronPerStep(G4int maxNumDriftElectron) +{ +/// Limit step to the specified maximum number of Cherenkov photons + + fMaxNumDriftElectron = maxNumDriftElectron; +} + +void G4DriftElectronPhysics::SetMaxBetaChangePerStep(G4double maxBetaChange) +{ +/// Limit step to the specified maximum change of beta of the parent particle + + fMaxBetaChange = maxBetaChange; +} + +void G4DriftElectronPhysics::SetTrackSecondariesFirst(G4DriftElectronProcessIndex index, + G4bool trackSecondariesFirst) +{ + if ( index >= kDENoProcess ) return; + if ( fProcessTrackSecondariesFirst[index] == trackSecondariesFirst ) return; + fProcessTrackSecondariesFirst[index] = trackSecondariesFirst; +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... diff --git a/NPSimulation/Process/G4DriftElectronPhysics.hh b/NPSimulation/Process/G4DriftElectronPhysics.hh new file mode 100644 index 000000000..d15b3175c --- /dev/null +++ b/NPSimulation/Process/G4DriftElectronPhysics.hh @@ -0,0 +1,116 @@ +// +// ******************************************************************** +// * License and Disclaimer * +// * * +// * The Geant4 software is copyright of the Copyright Holders of * +// * the Geant4 Collaboration. It is provided under the terms and * +// * conditions of the Geant4 Software License, included in the file * +// * LICENSE and available at http://cern.ch/geant4/license . These * +// * include a list of copyright holders. * +// * * +// * Neither the authors of this software system, nor their employing * +// * institutes,nor the agencies providing financial support for this * +// * work make any representation or warranty, express or implied, * +// * regarding this software system or assume any liability for its * +// * use. Please see the license in the file LICENSE and URL above * +// * for the full disclaimer and the limitation of liability. * +// * * +// * This code implementation is the result of the scientific and * +// * technical work of the GEANT4 collaboration. * +// * By using, copying, modifying or distributing the software (or * +// * any work based on the software) you agree to acknowledge its * +// * use in resulting scientific publications, and indicate your * +// * acceptance of all terms of the Geant4 Software license. * +// ******************************************************************** +// +// +// +//--------------------------------------------------------------------------- +// +// ClassName: G4DriftElectronPhysics +// +// Author: Adrien Matta 2016-12-19 +// +// Modified: +// +// +//--------------------------------------------------------------------------- +// +// This class provides construction of default optical physics +// + +#ifndef G4DriftElectronPhysics_h +#define G4DriftElectronPhysics_h 1 + +#include "globals.hh" + +#include "G4VPhysicsConstructor.hh" + +#include "G4DriftElectronProcessIndex.hh" + +#include <vector> + +class G4VProcess; +class G4EmSaturation; + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +class G4DriftElectronPhysics : public G4VPhysicsConstructor +{ + public: + + G4DriftElectronPhysics(G4int verbose = 0, const G4String& name = "DriftElectron"); + virtual ~G4DriftElectronPhysics(); + + protected: + + // construct particle and physics + virtual void ConstructParticle(); + virtual void ConstructProcess(); + + private: + + /// Not implemented + G4DriftElectronPhysics(const G4DriftElectronPhysics& right); + /// Not implemented + G4DriftElectronPhysics& operator=(const G4DriftElectronPhysics& right); + + public: + + // configure G4DriftElectronPhysics builder + void Configure(G4DriftElectronProcessIndex, G4bool ); + + void SetMaxNumDriftElectronPerStep(G4int ); + void SetMaxBetaChangePerStep(G4double ); + void SetElectronYieldFactor(G4double ); + void SetTrackSecondariesFirst(G4DriftElectronProcessIndex, G4bool ); + + private: + + // methods + void PrintStatistics() const; + + // The vector of process configuration + std::vector<G4bool> fProcessUse; + + // The vector of track secondaries options; + // the option to track secondaries before finishing their parent track + std::vector<G4bool> fProcessTrackSecondariesFirst; + + /// max number of Cerenkov photons per step + G4int fMaxNumDriftElectron; + + /// max change of beta per step + G4double fMaxBetaChange; + + /// scintillation yield factor + G4double fYieldFactor; + + /// option to allow for G4ElectronTrackInformation + /// to be attached to a drift electron's track + G4bool fElectronTrackInfo; +}; + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +#endif // G4DriftElectronPhysics_h diff --git a/NPSimulation/Process/G4DriftElectronProcessIndex.hh b/NPSimulation/Process/G4DriftElectronProcessIndex.hh new file mode 100644 index 000000000..9593a740b --- /dev/null +++ b/NPSimulation/Process/G4DriftElectronProcessIndex.hh @@ -0,0 +1,76 @@ +// +// ******************************************************************** +// * License and Disclaimer * +// * * +// * The Geant4 software is copyright of the Copyright Holders of * +// * the Geant4 Collaboration. It is provided under the terms and * +// * conditions of the Geant4 Software License, included in the file * +// * LICENSE and available at http://cern.ch/geant4/license . These * +// * include a list of copyright holders. * +// * * +// * Neither the authors of this software system, nor their employing * +// * institutes,nor the agencies providing financial support for this * +// * work make any representation or warranty, express or implied, * +// * regarding this software system or assume any liability for its * +// * use. Please see the license in the file LICENSE and URL above * +// * for the full disclaimer and the limitation of liability. * +// * * +// * This code implementation is the result of the scientific and * +// * technical work of the GEANT4 collaboration. * +// * By using, copying, modifying or distributing the software (or * +// * any work based on the software) you agree to acknowledge its * +// * use in resulting scientific publications, and indicate your * +// * acceptance of all terms of the Geant4 Software license. * +// ******************************************************************** +// +// +// +//--------------------------------------------------------------------------- +// +// ClassName: G4DriftElectronProcessIndex +// +// Author: Adrien Matta 2016-12-19 +// +// Modified: +// +//---------------------------------------------------------------------------- +// +// Enumeration for processes defined in G4DriftElectronPhysics +// +// The enumeration constants are used as the indices of the instatianted +// processes in the vector of processes; needed to configure G4DriftElectronPhysics +// (PhysicsList) according to selected process choices. +// + +#ifndef G4DriftElectronProcessIndex_h +#define G4DriftElectronProcessIndex_h 1 + +#include "globals.hh" + +enum G4DriftElectronProcessIndex { + kDEIonization, ///< ionisation process index + kDEAmplification, /// < Amplification of electron + kDETransport, // < Transport of Electron in medium + kDEAbsorption, ///< Absorption process index + kDENoProcess ///< no selected process +}; + +/// Return the name for a given optical process index +G4String G4DriftElectronProcessName(G4int ); + +//////////////////// +// Inline methods +//////////////////// + +inline +G4String G4DriftElectronProcessName(G4int processNumber) +{ + switch ( processNumber ) { + case kDEIonization: return "Ionization"; + case kDEAmplification: return "Amplification"; + case kDEAbsorption: return "Absorption"; + default: return "NoProcess"; + } +} + +#endif // G4DriftElectronProcessIndex_h diff --git a/NPSimulation/Process/G4IonizationWithDE.cc b/NPSimulation/Process/G4IonizationWithDE.cc new file mode 100644 index 000000000..2856283b4 --- /dev/null +++ b/NPSimulation/Process/G4IonizationWithDE.cc @@ -0,0 +1,248 @@ +// +// ******************************************************************** +// * License and Disclaimer * +// * * +// * The Geant4 software is copyright of the Copyright Holders of * +// * the Geant4 Collaboration. It is provided under the terms and * +// * conditions of the Geant4 Software License, included in the file * +// * LICENSE and available at http://cern.ch/geant4/license . These * +// * include a list of copyright holders. * +// * * +// * Neither the authors of this software system, nor their employing * +// * institutes,nor the agencies providing financial support for this * +// * work make any representation or warranty, express or implied, * +// * regarding this software system or assume any liability for its * +// * use. Please see the license in the file LICENSE and URL above * +// * for the full disclaimer and the limitation of liability. * +// * * +// * This code implementation is the result of the scientific and * +// * technical work of the GEANT4 collaboration. * +// * By using, copying, modifying or distributing the software (or * +// * any work based on the software) you agree to acknowledge its * +// * use in resulting scientific publications, and indicate your * +// * acceptance of all terms of the Geant4 Software license. * +// ******************************************************************** +// +// +// $Id: G4IonizationWithDE.hh 98002 2016-06-30 13:03:36Z gcosmo $ +// +// +//////////////////////////////////////////////////////////////////////// +// IonizationWithDE Light Class Definition +//////////////////////////////////////////////////////////////////////// +// +// File: G4IonizationWithDE.hh +// Description: Discrete Process - Generation of Drift Electron +// Version: 1.0 +// Created: 2016-12-19 +// Author: Adrien Matta +// Updated: +// +// +// mail: matta@lpccaen.in2p3.fr +// +//////////////////////////////////////////////////////////////////////// + +#include "G4ios.hh" +#include "globals.hh" +#include "G4PhysicalConstants.hh" +#include "G4SystemOfUnits.hh" +#include "G4ParticleTypes.hh" +#include "G4EmProcessSubType.hh" + +#include "G4IonizationWithDE.hh" +#include <iostream> +using namespace std; +///////////////////////// +// Class Implementation +///////////////////////// + +////////////////////// +// static data members +////////////////////// + + +//G4bool G4IonizationWithDE::fTrackSecondariesFirst = false; +//G4bool G4IonizationWithDE::fIonizationWithDEByParticleType = false; +//G4bool G4IonizationWithDE::fStackingFlag = true; +//G4EmSaturation* G4IonizationWithDE::fEmSaturation = NULL; + +////////////// +// Operators +////////////// + +// G4IonizationWithDE::operator=(const G4IonizationWithDE &right) +// { +// } + +///////////////// +// Constructors +///////////////// + +G4IonizationWithDE::G4IonizationWithDE(const G4String& processName, + G4ProcessType type) + : G4VRestDiscreteProcess(processName, type) , + fTrackSecondariesFirst(false), + fStackingFlag(true), + fNumDriftElectron(1e6), + fEmSaturation(nullptr) +{ + SetProcessSubType(fIonisation); // left as it is + + if (verboseLevel>0) { + G4cout << GetProcessName() << " is created " << G4endl; + } +} + +//////////////// +// Destructors +//////////////// + +G4IonizationWithDE::~G4IonizationWithDE() +{ +} + +//////////// +// Methods +//////////// + +void G4IonizationWithDE::BuildPhysicsTable(const G4ParticleDefinition&) +{ +} + + +// AtRestDoIt +// ---------- +// + G4VParticleChange* +G4IonizationWithDE::AtRestDoIt(const G4Track& aTrack, const G4Step& aStep) + + // This routine simply calls the equivalent PostStepDoIt since all the + // necessary information resides in aStep.GetTotalEnergyDeposit() + +{ + return G4IonizationWithDE::PostStepDoIt(aTrack, aStep); +} + +// PostStepDoIt +// ------------- +// + G4VParticleChange* +G4IonizationWithDE::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep) + + // This routine is called for each tracking step of a charged particle + // in a drift chamber gas. an electron is created for each pair energy deposited in + // the material. + +{ + // Get the primary track + aParticleChange.Initialize(aTrack); + + const G4Material* aMaterial = aTrack.GetMaterial(); + + G4StepPoint* pPreStepPoint = aStep.GetPreStepPoint(); + + G4ThreeVector x0 = pPreStepPoint->GetPosition(); + G4ThreeVector p0 = aStep.GetDeltaPosition().unit(); + G4double t0 = pPreStepPoint->GetGlobalTime(); + + G4double TotalEnergyDeposit = aStep.GetTotalEnergyDeposit(); + + // Get the material table + G4MaterialPropertiesTable* aMaterialPropertiesTable = + aMaterial->GetMaterialPropertiesTable(); + if (!aMaterialPropertiesTable) + return G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep); + + + if(!aMaterialPropertiesTable->ConstPropertyExists("DE_PAIRENERGY") || + !aMaterialPropertiesTable->ConstPropertyExists("DE_YIELD") || + !aMaterialPropertiesTable->ConstPropertyExists("DE_TRANSVERSALSPREAD") || + !aMaterialPropertiesTable->ConstPropertyExists("DE_LONGITUDINALSPREAD") || + !aMaterialPropertiesTable->ConstPropertyExists("DE_DRIFTSPEED") ) + return G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep); + + + G4double pair_energy=0; + pair_energy= + aMaterialPropertiesTable->GetConstProperty("DE_PAIRENERGY"); + G4double IonizationWithDEYield = 0; + IonizationWithDEYield= + aMaterialPropertiesTable->GetConstProperty("DE_YIELD"); + + G4int number_electron = IonizationWithDEYield*TotalEnergyDeposit/pair_energy; + //if no electron leave + if(number_electron<1){ + aParticleChange.SetNumberOfSecondaries(0); + return G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep); + } + + aParticleChange.SetNumberOfSecondaries(number_electron); + + // Create the secondary tracks + for(G4int i = 0 ; i < number_electron ; i++){ + // Random direction at creation + G4double cost = 1-2*G4UniformRand(); + G4double theta = acos(cost); + G4double phi = twopi*G4UniformRand(); + G4ThreeVector p; + p.setRThetaPhi(1,theta,phi); + + // Random Position along the step with matching time + G4double rand = G4UniformRand(); + G4ThreeVector pos = x0 + rand * aStep.GetDeltaPosition(); + G4double time = t0+ rand* aStep.GetDeltaTime(); + + G4DynamicParticle* particle = new G4DynamicParticle(G4DriftElectron::DriftElectron(),p, pair_energy); + G4Track* aSecondaryTrack = new G4Track(particle,time,pos); + + aSecondaryTrack->SetTouchableHandle( + aStep.GetPreStepPoint()->GetTouchableHandle()); + + aSecondaryTrack->SetParentID(aTrack.GetTrackID()); + aSecondaryTrack->SetTouchableHandle(aStep.GetPreStepPoint()->GetTouchableHandle()); + aParticleChange.AddSecondary(aSecondaryTrack); + } + + if (fTrackSecondariesFirst) { + if (aTrack.GetTrackStatus() == fAlive ) + aParticleChange.ProposeTrackStatus(fSuspend); + } + return G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep); +} + +// BuildThePhysicsTable for the scintillation process +// -------------------------------------------------- +// + +void G4IonizationWithDE::BuildThePhysicsTable() +{ +} + +// GetMeanFreePath +// --------------- +// + +G4double G4IonizationWithDE::GetMeanFreePath(const G4Track&, + G4double , + G4ForceCondition* condition) +{ + *condition = StronglyForced; + + return DBL_MAX; + +} + +// GetMeanLifeTime +// --------------- +// + +G4double G4IonizationWithDE::GetMeanLifeTime(const G4Track&, + G4ForceCondition* condition) +{ + *condition = StronglyForced; + + return DBL_MAX; + +} + diff --git a/NPSimulation/Process/G4IonizationWithDE.hh b/NPSimulation/Process/G4IonizationWithDE.hh new file mode 100644 index 000000000..8d2c2862a --- /dev/null +++ b/NPSimulation/Process/G4IonizationWithDE.hh @@ -0,0 +1,254 @@ +// +// ******************************************************************** +// * License and Disclaimer * +// * * +// * The Geant4 software is copyright of the Copyright Holders of * +// * the Geant4 Collaboration. It is provided under the terms and * +// * conditions of the Geant4 Software License, included in the file * +// * LICENSE and available at http://cern.ch/geant4/license . These * +// * include a list of copyright holders. * +// * * +// * Neither the authors of this software system, nor their employing * +// * institutes,nor the agencies providing financial support for this * +// * work make any representation or warranty, express or implied, * +// * regarding this software system or assume any liability for its * +// * use. Please see the license in the file LICENSE and URL above * +// * for the full disclaimer and the limitation of liability. * +// * * +// * This code implementation is the result of the scientific and * +// * technical work of the GEANT4 collaboration. * +// * By using, copying, modifying or distributing the software (or * +// * any work based on the software) you agree to acknowledge its * +// * use in resulting scientific publications, and indicate your * +// * acceptance of all terms of the Geant4 Software license. * +// ******************************************************************** +// +// +// $Id: G4IonizationWithDE.hh 98002 2016-06-30 13:03:36Z gcosmo $ +// +// +//////////////////////////////////////////////////////////////////////// +// IonizationWithDE Light Class Definition +//////////////////////////////////////////////////////////////////////// +// +// File: G4IonizationWithDE.hh +// Description: Discrete Process - Generation of Drift Electron +// Version: 1.0 +// Created: 2016-12-19 +// Author: Adrien Matta +// Updated: +// +// +// mail: matta@lpccaen.in2p3.fr +// +//////////////////////////////////////////////////////////////////////// + +#ifndef G4IonizationWithDE_h +#define G4IonizationWithDE_h 1 + +///////////// +// Includes +///////////// + +#include "globals.hh" +#include "templates.hh" +#include "Randomize.hh" +#include "G4Poisson.hh" +#include "G4ThreeVector.hh" +#include "G4ParticleMomentum.hh" +#include "G4Step.hh" +#include "G4VRestDiscreteProcess.hh" +#include "G4DriftElectron.hh" +#include "G4DynamicParticle.hh" +#include "G4Material.hh" +#include "G4PhysicsTable.hh" +#include "G4MaterialPropertiesTable.hh" +#include "G4PhysicsOrderedFreeVector.hh" + +#include "G4EmSaturation.hh" +#include<iostream> +// Class Description: +// RestDiscrete Process - Generation of IonizationWithDE Photons. +// Class inherits publicly from G4VRestDiscreteProcess. +// Class Description - End: + +///////////////////// +// Class Definition +///////////////////// + +class G4IonizationWithDE : public G4VRestDiscreteProcess +{ + +public: + + //////////////////////////////// + // Constructors and Destructor + //////////////////////////////// + + explicit G4IonizationWithDE(const G4String& processName = "IonizationWithDE", + G4ProcessType type = fElectromagnetic); + ~G4IonizationWithDE(); + +private: + + G4IonizationWithDE(const G4IonizationWithDE &right) = delete; + + ////////////// + // Operators + ////////////// + + G4IonizationWithDE& operator=(const G4IonizationWithDE &right) = delete; + +public: + + //////////// + // Methods + //////////// + + // G4IonizationWithDE Process has both PostStepDoIt (for energy + // deposition of particles in flight) and AtRestDoIt (for energy + // given to the medium by particles at rest) + + G4bool IsApplicable( + const G4ParticleDefinition& aParticleType) override; + // Returns true -> 'is applicable', for any particle type except + // for an 'driftelectron' and for short-lived particles + + void BuildPhysicsTable( + const G4ParticleDefinition& aParticleType) override; + // Build table at the right time + + G4double GetMeanFreePath(const G4Track& aTrack, + G4double , + G4ForceCondition* ) override; + // Returns infinity; i. e. the process does not limit the step, + // but sets the 'StronglyForced' condition for the DoIt to be + // invoked at every step. + + G4double GetMeanLifeTime(const G4Track& aTrack, + G4ForceCondition* ) override; + // Returns infinity; i. e. the process does not limit the time, + // but sets the 'StronglyForced' condition for the DoIt to be + // invoked at every step. + + G4VParticleChange* PostStepDoIt(const G4Track& aTrack, + const G4Step& aStep) override; + G4VParticleChange* AtRestDoIt (const G4Track& aTrack, + const G4Step& aStep) override; + + // These are the methods implementing the drift electron process. + + void SetTrackSecondariesFirst(const G4bool state); + // If set, the primary particle tracking is interrupted and any + // produced drift electron are tracked next. When all + // have been tracked, the tracking of the primary resumes. + + G4bool GetTrackSecondariesFirst() const; + // Returns the boolean flag for tracking secondaries first. + + G4PhysicsTable* GetIntegralTable() const; + // Returns the address of the ionisation integral table. + + void AddSaturation(G4EmSaturation* sat); + // Adds Birks Saturation to the process. + + void RemoveSaturation(); + // Removes the Birks Saturation from the process. + + G4EmSaturation* GetSaturation() const; + // Returns the Birks Saturation. + + void SetStackDriftElectron(const G4bool ); + // Call by the user to set the flag for stacking the drift electron + + G4bool GetStackDriftElectron() const; + // Return the boolean for whether or not the drift electron are stacked + + G4int GetNumDriftElectron() const; + // Returns the current number of drift electron (after PostStepDoIt) + + void DumpPhysicsTable() const; + // Prints the drift electron integral tables. + +protected: + + void BuildThePhysicsTable(); + // It builds either the drift electronintegral table; + + /////////////////////// + // Class Data Members + /////////////////////// + +private: + + G4bool fTrackSecondariesFirst; + G4bool fStackingFlag; + G4int fNumDriftElectron; + + G4double single_exp(G4double t, G4double tau2); + G4double bi_exp(G4double t, G4double tau1, G4double tau2); + G4EmSaturation* fEmSaturation; + +}; + +//////////////////// +// Inline methods +//////////////////// + +inline +G4bool G4IonizationWithDE::IsApplicable(const G4ParticleDefinition& aParticleType) +{ + if (aParticleType.GetParticleName() == "driftelectron") return false; + return true; +} + +inline +void G4IonizationWithDE::SetTrackSecondariesFirst(const G4bool state) +{ + fTrackSecondariesFirst = state; +} + +inline +G4bool G4IonizationWithDE::GetTrackSecondariesFirst() const +{ + return fTrackSecondariesFirst; +} + +inline +void G4IonizationWithDE::AddSaturation(G4EmSaturation* sat) +{ + fEmSaturation = sat; +} + +inline +void G4IonizationWithDE::RemoveSaturation() +{ + fEmSaturation = nullptr; +} + +inline +G4EmSaturation* G4IonizationWithDE::GetSaturation() const +{ + return fEmSaturation; +} + +inline +void G4IonizationWithDE::SetStackDriftElectron(const G4bool stackingFlag) +{ + fStackingFlag = stackingFlag; +} + +inline +G4bool G4IonizationWithDE::GetStackDriftElectron() const +{ + return fStackingFlag; +} + +inline +G4int G4IonizationWithDE::GetNumDriftElectron() const +{ + return fNumDriftElectron; +} + + +#endif /* G4IonizationWithDE_h */ diff --git a/NPSimulation/Process/PhysicsList.cc b/NPSimulation/Process/PhysicsList.cc index 38329420e..6cc59d9a5 100644 --- a/NPSimulation/Process/PhysicsList.cc +++ b/NPSimulation/Process/PhysicsList.cc @@ -115,7 +115,12 @@ PhysicsList::PhysicsList() : G4VUserPhysicsList(){ opticalPhysicsList->SetTrackSecondariesFirst(kScintillation,true); opticalPhysicsList->SetTrackSecondariesFirst(kCerenkov,true); } - + + if(m_DriftElectronPhysics){ + + driftElectronPhysicsList = new G4DriftElectronPhysics(0); + driftElectronPhysicsList->SetMaxNumDriftElectronPerStep(1e6); + } // Decay physics // Add Radioactive decay @@ -142,6 +147,7 @@ void PhysicsList::ReadConfiguration(std::string filename){ m_NPIonInelasticPhysics = 0; m_StoppingPhysics = 0; m_OpticalPhysics = 0; + m_DriftElectronPhysics = 0; m_HadronPhysicsQGSP_BIC_HP = 0; m_HadronPhysicsINCLXX = 0; m_Decay = 0; @@ -175,6 +181,8 @@ void PhysicsList::ReadConfiguration(std::string filename){ m_StoppingPhysics= value; else if (name == "OpticalPhysics") m_OpticalPhysics= value; + else if (name == "DriftElectronPhysics") + m_DriftElectronPhysics= value; else if (name == "HadronPhysicsQGSP_BIC_HP") m_HadronPhysicsQGSP_BIC_HP= value; else if (name == "HadronPhysicsINCLXX") @@ -204,7 +212,11 @@ void PhysicsList::ConstructParticle(){ ((G4VPhysicsConstructor*) opticalPhysicsList)->ConstructParticle(); } - + + if(driftElectronPhysicsList){ + ((G4VPhysicsConstructor*) driftElectronPhysicsList)->ConstructParticle(); + } + if(decay_List){ decay_List -> ConstructParticle(); radioactiveDecay_List->ConstructParticle(); @@ -244,7 +256,6 @@ void PhysicsList::ConstructParticle(){ ionConstructor.ConstructParticle() ; } } - ///////////////////////////////////////////////////////////////////////////// void PhysicsList::ConstructProcess(){ // Transportation @@ -255,6 +266,10 @@ void PhysicsList::ConstructProcess(){ if(opticalPhysicsList){ ((G4VPhysicsConstructor*) opticalPhysicsList)->ConstructProcess(); } + if(driftElectronPhysicsList){ + ((G4VPhysicsConstructor*) driftElectronPhysicsList)->ConstructProcess(); + } + // Hadronic physics std::map<std::string,G4VPhysicsConstructor*>::iterator it; for(it = m_PhysList.begin(); it!= m_PhysList.end(); it++){ diff --git a/NPSimulation/Process/PhysicsList.hh b/NPSimulation/Process/PhysicsList.hh index f291e314c..c171bc627 100644 --- a/NPSimulation/Process/PhysicsList.hh +++ b/NPSimulation/Process/PhysicsList.hh @@ -52,6 +52,9 @@ //Optical #include "G4OpticalPhysics.hh" +// Drift Electron +#include "G4DriftElectronPhysics.hh" + //Hadronique #include "G4IonElasticPhysics.hh" @@ -89,6 +92,7 @@ class PhysicsList: public G4VUserPhysicsList{ private: // Physics List G4OpticalPhysics* opticalPhysicsList; + G4DriftElectronPhysics* driftElectronPhysicsList; G4VPhysicsConstructor* emPhysicsList; G4VPhysicsConstructor* decay_List; G4VPhysicsConstructor* radioactiveDecay_List; @@ -102,6 +106,7 @@ class PhysicsList: public G4VUserPhysicsList{ double m_HadronInelasticPhysics; double m_StoppingPhysics; double m_OpticalPhysics; + double m_DriftElectronPhysics; double m_HadronPhysicsQGSP_BIC_HP; double m_HadronPhysicsINCLXX; double m_Decay; diff --git a/Projects/BeLise/34Si.beam b/Projects/BeLise/34Si.beam new file mode 100644 index 000000000..6a4ae9bac --- /dev/null +++ b/Projects/BeLise/34Si.beam @@ -0,0 +1,19 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%% Reaction file for 11Li(d,3He)10He reaction %%%%%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +Beam + Particle= 34Si + Energy= 6800 MeV + SigmaEnergy= 0 MeV + SigmaThetaX= 0. deg + SigmaPhiY= 0. deg + SigmaX= 0 mm + SigmaY= 0 mm + MeanThetaX= 0 deg + MeanPhiY= 0 deg + MeanX= 0 mm + MeanY= 0 mm + %EnergyProfilePath= + %XThetaXProfilePath= + %YPhiYProfilePath= +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% diff --git a/Projects/BeLise/chio.detector b/Projects/BeLise/chio.detector new file mode 100644 index 000000000..c7cfe4ce9 --- /dev/null +++ b/Projects/BeLise/chio.detector @@ -0,0 +1,15 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +Target + THICKNESS= 10 micrometer + RADIUS= 20 mm + MATERIAL= CD2 + ANGLE= 0 deg + X= 0 mm + Y= 0 mm + Z= 0 mm +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +Chio + POS= 0 0 350 mm + Shape= Square +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + -- GitLab