Вычисление методом Монте-Карло. ВЫЧИСЛЕНИЕ МЕТОДОМ МОНТЕ-КАРЛО КОЛИЧЕСТВА ПЕРВИЧНЫХ ДЕФЕКТОВ В К. вычисление методом монтекарло количества первичных дефектов в кристаллическом кремнии, создаваемых протонами и нейтронами различных энергий

• Название: вычисление методом монтекарло количества первичных дефектов в кристаллическом кремнии, создаваемых протонами и нейтронами различных энергий
• Анкор: Вычисление методом Монте-Карло
• Дата: 13.02.2022
• Размер: 0.62 Mb.
• Имя файла: ВЫЧИСЛЕНИЕ МЕТОДОМ МОНТЕ-КАРЛО КОЛИЧЕСТВА ПЕРВИЧНЫХ ДЕФЕКТОВ В К.doc
• Тип: Пояснительная записка; #360814
• страница: 5 из 6

ExN02SteppingVerbose.hh:

 

class ExN02SteppingVerbose;

 

#ifndef ExN02SteppingVerbose_h

#define ExN02SteppingVerbose_h 1

 

#include "G4SteppingVerbose.hh"

 

class ExN02SteppingVerbose : public G4SteppingVerbose

{

 public:

  ExN02SteppingVerbose();

 ExN02SteppingVerbose();

  void StepInfo();

  void TrackingStarted();

};

#endif

 

ExN02TrackingAction.cc:

 

#include "ExN02TrackingAction.hh"

#include "G4TrackingManager.hh"

#include "G4Track.hh"

#include "G4UnitsTable.hh"

 

ExN02TrackingAction::ExN02TrackingAction()

{

        //Инициализация массивов для спектров

    N = 0;

    E_total=0;

    for(int j=0; j
    O16_spectr[j]=0;

    Ne20_spectr[j]=0; Ne22_spectr[j]=0;

    Na22_spectr[j]=0; Na23_spectr[j]=0;

    Mg23_spectr[j]=0; Mg24_spectr[j]=0; Mg25_spectr[j]=0; Mg26_spectr[j]=0; Mg27_spectr[j]=0;

    Al26_spectr[j]=0; Al27_spectr[j]=0; Al28_spectr[j]=0; Al29_spectr[j]=0;

    Si27_spectr[j]=0; Si28_spectr[j]=0; Si29_spectr[j]=0; Si30_spectr[j]=0; Si31_spectr[j]=0;

    P29_spectr[j]=0;

    Elastic_spectr[j]=0;

    Inelastic_spectr[j]=0;

    }

}

 

void ExN02TrackingAction::PreUserTrackingAction(const G4Track* aTrack)

{

    double energy;

    G4Track* Track = (G4Track*)aTrack;

 

    name = aTrack->GetDefinition()->GetParticleName();

    energy = aTrack->GetKineticEnergy();

    if ((name=="neutron" && energy <= 1.*keV) ||

     (name!="neutron" && name!="proton")) {

        Track->SetKineticEnergy(0);

    }

 

    if (name == "O16[0.0]") {

        PutIntoSpectr(aTrack, O16_spectr, energy);

    } else

    if (name == "Ne20[0.0]") {

        PutIntoSpectr(aTrack, Ne20_spectr, energy);

    } else

    if (name == "Ne22[0.0]") {

        PutIntoSpectr(aTrack, Ne22_spectr, energy);

    } else

    if (name == "Na22[0.0]") {

        PutIntoSpectr(aTrack, Na22_spectr, energy);

    } else

    if (name == "Na23[0.0]") {

        PutIntoSpectr(aTrack, Na23_spectr, energy);

    } else

    if (name == "Mg23[0.0]") {

        PutIntoSpectr(aTrack, Mg23_spectr, energy);

    } else

    if (name == "Mg24[0.0]") {

        PutIntoSpectr(aTrack, Mg24_spectr, energy);

    } else

    if (name == "Mg25[0.0]") {

        PutIntoSpectr(aTrack, Mg25_spectr, energy);

    } else

    if (name == "Mg26[0.0]") {

        PutIntoSpectr(aTrack, Mg26_spectr, energy);

    } else

    if (name == "Mg27[0.0]") {

        PutIntoSpectr(aTrack, Mg27_spectr, energy);

    } else

    if (name == "Al26[0.0]") {

        PutIntoSpectr(aTrack, Al26_spectr, energy);

    } else

    if (name == "Al27[0.0]") {

        PutIntoSpectr(aTrack, Al27_spectr, energy);

    } else

    if (name == "Al28[0.0]") {

        PutIntoSpectr(aTrack, Al28_spectr, energy);

    } else

    if (name == "Al29[0.0]") {

        PutIntoSpectr(aTrack, Al29_spectr, energy);

    } else

    if (name == "Si27[0.0]") {

        PutIntoSpectr(aTrack, Si27_spectr, energy);

    } else

    if (name == "Si28[0.0]") {

        PutIntoSpectr(aTrack, Si28_spectr, energy);

    } else

    if (name == "Si29[0.0]") {

        PutIntoSpectr(aTrack, Si29_spectr, energy);

    } else

    if (name == "Si30[0.0]") {

        PutIntoSpectr(aTrack, Si30_spectr, energy);

    } else

    if (name == "Si31[0.0]") {

        PutIntoSpectr(aTrack, Si31_spectr, energy);

    } else

    if (name == "P29[0.0]") {

        PutIntoSpectr(aTrack, P29_spectr, energy);

    }   

}

 

void ExN02TrackingAction::PutIntoSpectr(const G4Track* aTrack, t_spectr spectr, double energy)

{

    FILE *f;

    int bin;

    N++;

   

    if (aTrack->GetCreatorProcess())

        proc = aTrack->GetCreatorProcess()->GetProcessName();  //Получили имя родительского процесса

    else

        proc = "starting";

   

    E_total = E_total + energy; //Полные потери энергии

    energy = energy * 1000 / BINSIZE;

    bin = (int)energy;  //Определили бин в спектре

 

    //заполняем спектры

    if (bin >=0 && bin <= MAX-1)

    if (proc=="LElastic") {

        Elastic_spectr[bin] = Elastic_spectr[bin] + 1;

    }

    else

    if (proc=="ProtonInelastic" || proc=="NeutronInelastic") {

        spectr[bin] = spectr[bin] + 1;

        Inelastic_spectr[bin] = Inelastic_spectr[bin] + 1;

    }

 

    if (N%10000 == 0) { //Если накопилось 10000 новых событий, то скидываем спектры в файл

        if ((f=fopen("spectr.txt", "w"))==NULL) {

            G4cout << "error Output!\n";

            exit(1);

        }

        G4cout << N << " putted.\n";

        fprintf(f, "%i\nE=%f MeV\n", N, E_total);

        fprintf(f,"Elastic(Si28,Si29,Si30):\n");

        for(int i=0; i
        fprintf(f,"Inelastic(All):\n");

        for(i=0; i
        fprintf(f,"Inelastic(O16):\n");

        for(i=0; i
        fprintf(f,"Inelastic(Ne20):\n");

        for(i=0; i
        fprintf(f,"Inelastic(Ne22):\n");

        for(i=0; i
        fprintf(f,"Inelastic(Na22):\n");

        for(i=0; i
        fprintf(f,"Inelastic(Na23):\n");

        for(i=0; i
        fprintf(f,"Inelastic(Mg23):\n");

        for(i=0; i
        fprintf(f,"Inelastic(Mg24):\n");

        for(i=0; i
        fprintf(f,"Inelastic(Mg25):\n");

        for(i=0; i
        fprintf(f,"Inelastic(Mg26):\n");

        for(i=0; i
        fprintf(f,"Inelastic(Mg27):\n");

        for(i=0; i
        fprintf(f,"Inelastic(Al26):\n");

        for(i=0; i
        fprintf(f,"Inelastic(Al27):\n");

        for(i=0; i
        fprintf(f,"Inelastic(Al28):\n");

        for(i=0; i
        fprintf(f,"Inelastic(Al29):\n");

        for(i=0; i
        fprintf(f,"Inelastic(Si27):\n");

        for(i=0; i
        fprintf(f,"Inelastic(Si28):\n");

        for(i=0; i
        fprintf(f,"Inelastic(Si29):\n");

        for(i=0; i
        fprintf(f,"Inelastic(Si30):\n");

        for(i=0; i
        fprintf(f,"Inelastic(Si31):\n");

        for(i=0; i
        fprintf(f,"Inelastic(P29):\n");

        for(i=0; i
        fclose(f);

    }

   

    if (N>=EVMAX) { //Если набрали нужное количество событий, то выходим из программы

        G4cout << "Work done ! :)\n";

        exit(0);

    }

}

 

ExN02TrackingAction.hh:

 

#ifndef ExN02TrackingAction_h

#define ExN02TrackingAction_h 1

#include "G4UserTrackingAction.hh"

#include "G4String.hh"

#define MAX 300                            // Число каналов в спектре

#define EVMAX 1000000                               // общее число событий

#define BINSIZE 10                         // размер канала в кевах

 

typedef int t_spectr[MAX];

 

class ExN02TrackingAction : public G4UserTrackingAction {

 

  public:

    ExN02TrackingAction();

    virtual ExN02TrackingAction(){};

    virtual void PreUserTrackingAction(const G4Track* aTrack);

  private:

    void PutIntoSpectr(const G4Track* aTrack, t_spectr spectr, double energy);

    int N;

    double E_total;

    G4String name;

    G4String proc;

    t_spectr O16_spectr;

    t_spectr Ne20_spectr, Ne22_spectr;

    t_spectr Na22_spectr, Na23_spectr;

    t_spectr Mg23_spectr, Mg24_spectr, Mg25_spectr, Mg26_spectr, Mg27_spectr;

    t_spectr Al26_spectr, Al27_spectr, Al28_spectr, Al29_spectr;

    t_spectr Si27_spectr, Si28_spectr, Si29_spectr, Si30_spectr, Si31_spectr;

    t_spectr P29_spectr;

    t_spectr Elastic_spectr;

    t_spectr Inelastic_spectr;

};

#endif

 

ExN02VisManager.cc:

#ifdef G4VIS_USE

 

#include "ExN02VisManager.hh"

 

#include "G4ASCIITree.hh"

#include "G4DAWNFILE.hh"

#include "G4GAGTree.hh"

#include "G4HepRepFile.hh"

#include "G4RayTracer.hh"

#include "G4VRML1File.hh"

#include "G4VRML2File.hh"

 

// Внешние графические пакеты добавляются если определены соответствующие переменные окружения

#ifdef G4VIS_USE_DAWN

#include "G4FukuiRenderer.hh"

#endif

 

#ifdef G4VIS_USE_OPACS

#include "G4Wo.hh"

#include "G4Xo.hh"

#endif

 

#ifdef G4VIS_USE_OPENGLX

#include "G4OpenGLImmediateX.hh"

#include "G4OpenGLStoredX.hh"

#endif

 

#ifdef G4VIS_USE_OPENGLWIN32

#include "G4OpenGLImmediateWin32.hh"

#include "G4OpenGLStoredWin32.hh"

#endif

 

#ifdef G4VIS_USE_OPENGLXM

#include "G4OpenGLImmediateXm.hh"

#include "G4OpenGLStoredXm.hh"

#endif

 

#ifdef G4VIS_USE_OIX

#include "G4OpenInventorX.hh"

#endif

 

#ifdef G4VIS_USE_OIWIN32

#include "G4OpenInventorWin32.hh"

#endif

 

#ifdef G4VIS_USE_VRML

#include "G4VRML1.hh"

#include "G4VRML2.hh"

#endif

 

 

ExN02VisManager::ExN02VisManager () {}

 

void ExN02VisManager::RegisterGraphicsSystems () {

 

  RegisterGraphicsSystem (new G4ASCIITree);

  RegisterGraphicsSystem (new G4DAWNFILE);

  RegisterGraphicsSystem (new G4GAGTree);

  RegisterGraphicsSystem (new G4HepRepFile);

  RegisterGraphicsSystem (new G4RayTracer);

  RegisterGraphicsSystem (new G4VRML1File);

  RegisterGraphicsSystem (new G4VRML2File);

 

 

#ifdef G4VIS_USE_DAWN

  RegisterGraphicsSystem (new G4FukuiRenderer);

#endif

 

#ifdef G4VIS_USE_OPACS

  RegisterGraphicsSystem (new G4Wo);

  RegisterGraphicsSystem (new G4Xo);

#endif

 

#ifdef G4VIS_USE_OPENGLX

  RegisterGraphicsSystem (new G4OpenGLImmediateX);

  RegisterGraphicsSystem (new G4OpenGLStoredX);

#endif

 

#ifdef G4VIS_USE_OPENGLWIN32

  RegisterGraphicsSystem (new G4OpenGLImmediateWin32);

  RegisterGraphicsSystem (new G4OpenGLStoredWin32);

#endif

 

#ifdef G4VIS_USE_OPENGLXM

  RegisterGraphicsSystem (new G4OpenGLImmediateXm);

  RegisterGraphicsSystem (new G4OpenGLStoredXm);

#endif

 

#ifdef G4VIS_USE_OIX

  RegisterGraphicsSystem (new G4OpenInventorX);

#endif

 

#ifdef G4VIS_USE_OIWIN32

  RegisterGraphicsSystem (new G4OpenInventorWin32);

#endif

 

#ifdef G4VIS_USE_VRML

  RegisterGraphicsSystem (new G4VRML1);

  RegisterGraphicsSystem (new G4VRML2);

#endif

 

  if (fVerbose > 0) {

    G4cout <<

      "\nYou have successfully chosen to use the following graphics systems."

         << G4endl;

    PrintAvailableGraphicsSystems ();

  }

}

#endif

 

ExN02VisManager.hh:

#ifndef ExN02VisManager_h

#define ExN02VisManager_h 1

#ifdef G4VIS_USE

#include "G4VisManager.hh"

 

class ExN02VisManager: public G4VisManager {

public:

  ExN02VisManager ();

private:

  void RegisterGraphicsSystems ();

};

#endif

#endif

 

PhysicsList.cc:

Использует дополнительную библиотеку для описания взаимодействий протонов и нетронов средних энергий с веществом: LHEP_PRECO_HP (http://cmsdoc.cern.ch/hpw/GHAD/HomePage/lists/LHEP_PRECO_HP/index.html)

 

#include "PhysicsList.hh"

#include "globals.hh"

#include "G4ParticleDefinition.hh"

#include "G4ParticleWithCuts.hh"

#include "G4ProcessManager.hh"

#include "G4ProcessVector.hh"

#include "G4ParticleTypes.hh"

#include "G4ParticleTable.hh"

 

#include "G4Material.hh"

#include "G4MaterialTable.hh"

#include "G4ios.hh"

#include "g4std/iomanip"  

 

#include "GeneralPhysics.hh"

#include "EMPhysics.hh"

#include "MuonPhysics.hh"

#include "HadronPhysicsLHEP_PRECO_HP.hh"

#include "IonPhysics.hh"

 

PhysicsList::PhysicsList():  G4VModularPhysicsList()

{

  defaultCutValue = 1.0*mm; //граница обрезания по умолчанию

  SetVerboseLevel(1);

  RegisterPhysics( new GeneralPhysics("general") ); //Основные физические процессы

  RegisterPhysics( new EMPhysics("standard EM")); //Электро-магнитные взаимодействия

  RegisterPhysics(  new MuonPhysics("muon")); //Мюонные взаимодействия

  RegisterPhysics(  new HadronPhysicsLHEP_PRECO_HP("hadron")); //Физические процессы для нейтронов и протонов средних энергий

  RegisterPhysics( new IonPhysics("ion")); //Взаимодействия для ионов

}

PhysicsList::PhysicsList()

{

}

 

PhysicsList.hh:

#ifndef PhysicsList_h

#define PhysicsList_h 1

#include "G4VModularPhysicsList.hh"

#include "globals.hh"

#define kNuCut  5*m

class PhysicsList: public G4VModularPhysicsList

{

public:

  PhysicsList();

  virtual PhysicsList();

public:

  virtual void SetCuts();

};

#endif

 

ExN02SteppingAction.cc (new_cos):

#include "ExN02SteppingAction.hh"

#include "G4SteppingManager.hh"

 

 

ExN02SteppingAction::ExN02SteppingAction()

{

}

 

void ExN02SteppingAction::UserSteppingAction(const G4Step* aStep)

{

  double cost;

  int q, bin, i;

  FILE *fs;

  G4String proc, name;

 

    name = aStep->GetTrack()->GetDefinition()->GetParticleName();

   

    if (Act==0) {

        if (name==PARTICLE_INCIDENT) {

            proc = aStep->GetPostStepPoint()->GetProcessDefinedStep()->GetProcessName();

            if (proc == PROCESS_NAME) {

                        X1 = aStep->GetTrack()->GetPosition().x();

                        Y1 = aStep->GetTrack()->GetPosition().y();

                        Z1 = aStep->GetTrack()->GetPosition().z();

                        Act = 1;

            }

        }

    } else

 

    if (Act==1)

    {

        if (name==PARTICLE_PRODUCT1

#ifdef PARTICLE_PRODUCT2

            || name==PARTICLE_PRODUCT2

#endif

#ifdef PARTICLE_PRODUCT3

            || name==PARTICLE_PRODUCT3

#endif

            ) {

            X2 = aStep->GetTrack()->GetPosition().x();

            Y2 = aStep->GetTrack()->GetPosition().y();

            Z2 = aStep->GetTrack()->GetPosition().z();

            Act = 2;

            cost = (X2 - X1)/sqrt(sqr(X2-X1) + sqr(Y2-Y1) + sqr(Z2-Z1));

            if (cost >= 1.0000001 || cost <= -1.0000001) {

                        fprintf(stdout, "Super Mega Hyper Error ! cos=%f\n", cost);

                        fflush(stdout);

                        exit(1);

            }

            bin = cost*100 + 100;

            if (bin>=0 && bin
            SPECTR_cost[bin] = SPECTR_cost[bin] + 1;

            N++;

 

            if (N%100==0) {

                        G4cout << N << " putted.\n";

                        if ((fs=fopen("spectr_cos.txt", "w"))==NULL) {

                            printf("I/O Error !\n");

                            exit(1);

                        };

                        fprintf(fs, "Cos_t(%i): \n", N);

                        for(i=0; i
                          fprintf(fs, "%i\n", SPECTR_cost[i]);

                        fclose(fs);

                        if (N>=MAXEVT) {

                            printf("Work done !!!\n");

                            exit(0);

                        }

            }

        }

    }

}