|
Вычисление методом Монте-Карло. ВЫЧИСЛЕНИЕ МЕТОДОМ МОНТЕ-КАРЛО КОЛИЧЕСТВА ПЕРВИЧНЫХ ДЕФЕКТОВ В К. вычисление методом монтекарло количества первичных дефектов в кристаллическом кремнии, создаваемых протонами и нейтронами различных энергий
new_phys.cc:
#include "ExN02DetectorConstruction.hh"
#include "PhysicsList.hh"
#include "ExN02PrimaryGeneratorAction.hh"
#include "ExN02EventAction.hh"
#include "ExN02TrackingAction.hh"
#include "ExN02SteppingVerbose.hh"
#include "G4RunManager.hh"
#include "G4UImanager.hh"
#include "G4UIterminal.hh"
#include "G4UItcsh.hh"
#ifdef G4VIS_USE //Если установлена переменная окружения G4VIS_USE, то используем визуализацию
#include "ExN02VisManager.hh"
#endif
int main(int argc,char** argv) { //argc – число параметров командной строки, argv – массив строк с параметрами
srand((unsigned int)time((time_t *)NULL)); //Инициализация генератора случайных чисел
HepRandom::setTheSeed(rand()); // Инициализация генератора HepRandom для Geant4
G4VSteppingVerbose::SetInstance(new ExN02SteppingVerbose); //Указание способа формирования отчета
G4RunManager * runManager = new G4RunManager; //Инициализация менеджера загрузки
ExN02DetectorConstruction* ExN02detector = new ExN02DetectorConstruction; //Инициализация описания детектора
runManager->SetUserInitialization(ExN02detector);
runManager->SetUserInitialization(new PhysicsList);
#ifdef G4VIS_USE
G4VisManager* visManager = new ExN02VisManager; //Инициализация менеджера визуализации
visManager->Initialize();
#endif
// Инициализация пользовательских классов
runManager->SetUserAction(new ExN02PrimaryGeneratorAction(ExN02detector));
runManager->SetUserAction(new ExN02EventAction);
runManager->SetUserAction(new ExN02TrackingAction);
//Инициализация ядра Geant4
runManager->Initialize();
//получить указатель на менеджер пользовательского интерфейса
G4UImanager * UI = G4UImanager::GetUIpointer();
if(argc==1) //если в качестве параметров командной строки ничего не указано
// Определяем терминал пользователя для интерактивного режима
{
G4UIsession * session = 0;
session = new G4UIterminal();
UI->ApplyCommand("/control/execute vis.mac"); //Инициализируем визуализацию из командного
session->SessionStart();
delete session;
}
else
// иначе командный режим
{
G4String command = "/control/execute ";
G4String fileName = argv[1]; //получаем имя файла из командной строки
UI->ApplyCommand(command+fileName); //и выполняем его
}
#ifdef G4VIS_USE //Если используем визуализацию, то удаляем из памяти менеджер визуализации
delete visManager;
#endif
delete runManager;
return 0;
}
ExN02DetectorConstruction.cc:
#include "ExN02DetectorConstruction.hh"
#include "G4Material.hh"
#include "G4Box.hh"
#include "G4LogicalVolume.hh"
#include "G4PVPlacement.hh"
#include "G4PVParameterised.hh"
#include "G4SDManager.hh"
#include "G4UserLimits.hh"
#include "G4VisAttributes.hh"
#include "G4Colour.hh"
#include "G4ios.hh"
ExN02DetectorConstruction::ExN02DetectorConstruction() //Конструктор класса описания детектора
:solidWorld(0), logicWorld(0), physiWorld(0),
solidTarget(0), logicTarget(0), physiTarget(0),
solidTracker(0),logicTracker(0),physiTracker(0),
TargetMater(0), fWorldLength(0.), fTargetLength(0.), fTrackerLength(0.)
{
}
ExN02DetectorConstruction::ExN02DetectorConstruction()
{
}
G4VPhysicalVolume* ExN02DetectorConstruction::Construct()
{
//Определение материала детектора
G4double a, iz, z, density;
G4String name, symbol;
G4double temperature, pressure;
G4int nel;
density = 2.33*g/cm3; //плотность
a = 28.09*g/mole; //массовое число
// Природный кремний
G4Material * Si28nat = new G4Material(name="Silicon", z=14., a, density);
G4Box* solidWorld= new G4Box("world", 25*cm, 25*cm, 25*cm); //геометрический объем
logicWorld= new G4LogicalVolume( solidWorld, Si28nat, "World", 0, 0, 0); //логический объем
//физический объем
physiWorld = new G4PVPlacement(0, // нет вращения
G4ThreeVector(), // в точку (0,0,0)
"World", // имя объема
logicWorld, // логический объем
0, // не имеет родительских объемов
false, // логические операции не определены
0); // нет поля
//--------- Атрибуты визуализации -------------------------------
G4VisAttributes* BoxVisAtt= new G4VisAttributes(G4Colour(1.0,1.0,1.0));
logicWorld->SetVisAttributes(BoxVisAtt);
return physiWorld;
}
ExN02DetectorConstruction.hh:
#ifndef ExN02DetectorConstruction_h //проверка на то, что этот файл включен другими файлами проекта не больше чем 1 раз
#define ExN02DetectorConstruction_h 1
#include "globals.hh"
#include "G4VUserDetectorConstruction.hh"
class G4Box;
class G4LogicalVolume;
class G4VPhysicalVolume;
class G4Material;
class ExN02DetectorConstruction : public G4VUserDetectorConstruction
{
public:
ExN02DetectorConstruction();
ExN02DetectorConstruction();
public:
G4VPhysicalVolume* Construct();
const
G4VPhysicalVolume* GetTracker() {return physiTracker;};
G4double GetTrackerFullLength() {return fTrackerLength;};
G4double GetTargetFullLength() {return fTargetLength;};
G4double GetWorldFullLength() {return fWorldLength;};
private:
G4Box* solidWorld; // указатель на геометрический объем
G4LogicalVolume* logicWorld; // указатель на логический объем
G4VPhysicalVolume* physiWorld; // указатель на физический объем
};
#endif
ExN02EventAction.cc:
#include "ExN02EventAction.hh"
#include "G4Event.hh"
#include "G4EventManager.hh"
#include "G4TrajectoryContainer.hh"
#include "G4Trajectory.hh"
#include "G4VVisManager.hh"
#include "G4ios.hh"
ExN02EventAction::ExN02EventAction() //Конструктор
{}
ExN02EventAction::ExN02EventAction()//Деструктор
{}
void ExN02EventAction::BeginOfEventAction(const G4Event*) //Действия в начале события
{}
void ExN02EventAction::EndOfEventAction(const G4Event* evt) //Действия в конце события
{
G4int event_id = evt->GetEventID();
G4TrajectoryContainer* trajectoryContainer = evt->GetTrajectoryContainer();
G4int n_trajectories = 0;
if (trajectoryContainer) n_trajectories = trajectoryContainer->entries();
// вывод номера события с периодом 100
if (event_id%1000 == 0) {
G4cout << ">>> Event " << evt->GetEventID() << G4endl;
}
// нарисовать траектории с помощью менеджера визуализации
if (G4VVisManager::GetConcreteInstance())
{
for (G4int i=0; i { G4Trajectory* trj = (G4Trajectory*)
((*(evt->GetTrajectoryContainer()))[i]);
trj->DrawTrajectory(50);
}
}
}
ExN02EventAction.hh:
#ifndef ExN02EventAction_h
#define ExN02EventAction_h 1
#include "G4UserEventAction.hh"
class G4Event;
class ExN02EventAction : public G4UserEventAction
{
public:
ExN02EventAction();
ExN02EventAction();
public:
void BeginOfEventAction(const G4Event*);
void EndOfEventAction(const G4Event*);
};
#endif
ExN02PrimaryGeneratorAction.cc:
#include "ExN02PrimaryGeneratorAction.hh"
#include "ExN02DetectorConstruction.hh"
#include "G4Event.hh"
#include "G4ParticleGun.hh"
#include "G4ParticleTable.hh"
#include "G4ParticleDefinition.hh"
#include "globals.hh"
ExN02PrimaryGeneratorAction::ExN02PrimaryGeneratorAction(
ExN02DetectorConstruction* myDC)
:myDetector(myDC)
{
G4int n_particle = 1;
particleGun = new G4ParticleGun(n_particle);
G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();
G4ParticleDefinition* particle = particleTable->FindParticle("neutron"); //Задаем тип начальной частицы
particleGun->SetParticleDefinition(particle);
particleGun->SetParticleMomentumDirection(G4ThreeVector(1.,0.,0.)); //Задаем направление вылета начальной частицы
particleGun->SetParticleEnergy(14.0*MeV); //Задаем энергию начальной частицы
}
ExN02PrimaryGeneratorAction::ExN02PrimaryGeneratorAction()
{
delete particleGun;
}
void ExN02PrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent)
{
particleGun->SetParticlePosition(G4ThreeVector(-25*cm, 0.*cm, 0.*cm)); //Задаем позицию начальной частицы
particleGun->GeneratePrimaryVertex(anEvent);
}
ExN02PrimaryGeneratorAction.hh:
#ifndef ExN02PrimaryGeneratorAction_h
#define ExN02PrimaryGeneratorAction_h 1
#include "G4VUserPrimaryGeneratorAction.hh"
class ExN02DetectorConstruction;
class G4ParticleGun;
class G4Event;
class ExN02PrimaryGeneratorAction : public G4VUserPrimaryGeneratorAction
{
public:
ExN02PrimaryGeneratorAction(ExN02DetectorConstruction*);
ExN02PrimaryGeneratorAction();
public:
void GeneratePrimaries(G4Event*);
private:
G4ParticleGun* particleGun;
ExN02DetectorConstruction* myDetector;
};
#endif
ExN02SteppingVerbose.cc:
#include "ExN02SteppingVerbose.hh"
#include "G4SteppingManager.hh"
#include "G4UnitsTable.hh"
ExN02SteppingVerbose::ExN02SteppingVerbose()
{}
ExN02SteppingVerbose::ExN02SteppingVerbose()
{}
void ExN02SteppingVerbose::StepInfo() //Вывод информации о частице
{
CopyState();
G4int prec = G4cout.precision(3);
if( verboseLevel >= 1 ){
if( verboseLevel >= 4 ) VerboseTrack();
if( verboseLevel >= 3 ){ G4cout << G4endl;
G4cout << G4std::setw( 5) << "#Step#" << " "
<< G4std::setw( 6) << "X" << " "
<< G4std::setw( 6) << "Y" << " "
<< G4std::setw( 6) << "Z" << " "
<< G4std::setw( 9) << "KineE" << " "
<< G4std::setw( 9) << "dEStep" << " "
<< G4std::setw(10) << "StepLeng"
<< G4std::setw(10) << "TrakLeng"
<< G4std::setw(10) << "Volume" << " "
<< G4std::setw(10) << "Process" << G4endl;
}
G4cout << G4std::setw(5) << fTrack->GetCurrentStepNumber() << " "
<< G4std::setw(6) << G4BestUnit(fTrack->GetPosition().x(),"Length")
<< G4std::setw(6) << G4BestUnit(fTrack->GetPosition().y(),"Length")
<< G4std::setw(6) << G4BestUnit(fTrack->GetPosition().z(),"Length")
<< G4std::setw(6) << G4BestUnit(fTrack->GetKineticEnergy(),"Energy")
<< G4std::setw(6) << G4BestUnit(fStep->GetTotalEnergyDeposit(),"Energy")
<< G4std::setw(6) << G4BestUnit(fStep->GetStepLength(),"Length")
<< G4std::setw(6) << G4BestUnit(fTrack->GetTrackLength(),"Length");
if( fTrack->GetNextVolume() != 0 ) {
G4cout << G4std::setw(10) << fTrack->GetVolume()->GetName();
} else {
G4cout << G4std::setw(10) << "OutOfWorld";
}
if(fStep->GetPostStepPoint()->GetProcessDefinedStep() != NULL){
G4cout << " "
<< G4std::setw(10) << fStep->GetPostStepPoint()->GetProcessDefinedStep()
->GetProcessName();
} else {
G4cout << " UserLimit";
}
G4cout << G4endl;
if( verboseLevel == 2 ){
G4int tN2ndariesTot = fN2ndariesAtRestDoIt +
fN2ndariesAlongStepDoIt +
fN2ndariesPostStepDoIt;
if(tN2ndariesTot>0){
G4cout << " :----- List of 2ndaries - "
<< "#SpawnInStep=" << G4std::setw(3) << tN2ndariesTot
<< "(Rest=" << G4std::setw(2) << fN2ndariesAtRestDoIt
<< ",Along=" << G4std::setw(2) << fN2ndariesAlongStepDoIt
<< ",Post=" << G4std::setw(2) << fN2ndariesPostStepDoIt
<< "), "
<< "#SpawnTotal=" << G4std::setw(3) << (*fSecondary).size()
<< " ---------------"
<< G4endl;
for(size_t lp1=(*fSecondary).size()-tN2ndariesTot;
lp1<(*fSecondary).size(); lp1++){
G4cout << " : "
<< G4std::setw(6)
<< G4BestUnit((*fSecondary)[lp1]->GetPosition().x(),"Length")
<< G4std::setw(6)
<< G4BestUnit((*fSecondary)[lp1]->GetPosition().y(),"Length")
<< G4std::setw(6)
<< G4BestUnit((*fSecondary)[lp1]->GetPosition().z(),"Length")
<< G4std::setw(6)
<< G4BestUnit((*fSecondary)[lp1]->GetKineticEnergy(),"Energy")
<< G4std::setw(10)
<< (*fSecondary)[lp1]->GetDefinition()->GetParticleName();
G4cout << G4endl;
}
G4cout << " :-----------------------------"
<< "----------------------------------"
<< "-- EndOf2ndaries Info ---------------"
<< G4endl;
}
}
}
G4cout.precision(prec);
}
void ExN02SteppingVerbose::TrackingStarted()
{
CopyState();
G4int prec = G4cout.precision(3);
if( verboseLevel > 0 ){
G4cout << G4std::setw( 5) << "Step#" << " "
<< G4std::setw( 6) << "X" << " "
<< G4std::setw( 6) << "Y" << " "
<< G4std::setw( 6) << "Z" << " "
<< G4std::setw( 9) << "KineE" << " "
<< G4std::setw( 9) << "dEStep" << " "
<< G4std::setw(10) << "StepLeng"
<< G4std::setw(10) << "TrakLeng"
<< G4std::setw(10) << "Volume" << " "
<< G4std::setw(10) << "Process" << G4endl;
G4cout << G4std::setw(5) << fTrack->GetCurrentStepNumber() << " "
<< G4std::setw(6) << G4BestUnit(fTrack->GetPosition().x(),"Length")
<< G4std::setw(6) << G4BestUnit(fTrack->GetPosition().y(),"Length")
<< G4std::setw(6) << G4BestUnit(fTrack->GetPosition().z(),"Length")
<< G4std::setw(6) << G4BestUnit(fTrack->GetKineticEnergy(),"Energy")
<< G4std::setw(6) << G4BestUnit(fStep->GetTotalEnergyDeposit(),"Energy")
<< G4std::setw(6) << G4BestUnit(fStep->GetStepLength(),"Length")
<< G4std::setw(6) << G4BestUnit(fTrack->GetTrackLength(),"Length");
if(fTrack->GetNextVolume()){
G4cout << G4std::setw(10) << fTrack->GetVolume()->GetName();
} else {
G4cout << G4std::setw(10) << "OutOfWorld";
}
G4cout << " initStep" << G4endl;
}
G4cout.precision(prec);
}
|
|
|