ProtonMC/ProtonMC.cc at master · particleCT/ProtonMC · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
#include "PhysicsList.hh"
#include "PrimaryGeneratorAction.hh"
#include "OrganicMaterial.hh"
#include "G4NistManager.hh"
#include "DetectorConstruction.hh"
#include "SteppingAction.hh"
#include "G4EmCalculator.hh"
#include "G4IonTable.hh"
#include "G4RunManager.hh"
#include <TROOT.h>
#include "G4ExceptionHandler.hh"
#include "G4Material.hh"
#include "G4NistManager.hh"
#include "G4EmCalculator.hh"
#include "G4Proton.hh"
#include "G4UnitsTable.hh"
#include "Analysis.hh"
#include "G4SystemOfUnits.hh"
#include "G4PhysicalConstants.hh"
#include "G4UImanager.hh"

#ifdef VIS
#include "G4VisExecutive.hh"

#endif
#include "G4UIExecutive.hh"
#include <iostream>
#include <fstream>
using namespace std;
void calcRSP(PrimaryGeneratorAction *);
void calcStoppingPower(PrimaryGeneratorAction *);
int main(int argc,char** argv) {
  gROOT->ProcessLine("#include <vector>");
  time_t seed;
  time(&seed);
  CLHEP::RanecuEngine *theRanGenerator = new CLHEP::RanecuEngine;
  theRanGenerator->setSeed(seed);
  CLHEP::HepRandom::setTheEngine(theRanGenerator);
  if(argc<=6) {
    cout<<"Please input the following arguments : NParticules Energy Model Angle Thick Thread ANumber"<<endl;
    return 0;
  }
  G4int nProtons  = atoi(argv[1]);
  G4double Energy = atof(argv[2]); // MeV
  G4String Model  = argv[3];
  G4double angle  = atof(argv[4]); // Degrees
  G4double thick  = atof(argv[5]); // cm
  G4int   thread  = atoi(argv[6]);
  G4int  ANumber  = atoi(argv[7]); //Atomic Number

  G4String paraWorldName = "";
  G4RunManager* runManager = new G4RunManager;
  //G4ExceptionHandler* handler = new G4ExceptionHandler();
  runManager->SetUserInitialization(new PhysicsList(paraWorldName));
  DetectorConstruction* myDC = new DetectorConstruction(Model,angle,thick);
  PrimaryGeneratorAction *theGenerator =  new PrimaryGeneratorAction(Energy,ANumber);
  Analysis* theAnalysis = new Analysis(thread,angle,Model, nProtons);

  runManager->SetUserAction(theGenerator);
  runManager->SetUserAction( new SteppingAction() );
  runManager->SetUserInitialization( myDC );


  runManager->SetVerboseLevel(0);
  runManager->Initialize();

  //G4UImanager * UImanager = G4UImanager::GetUIpointer();
  //UImanager->ApplyCommand("/process/inactivate msc all");
  //UImanager->ApplyCommand("/process/inactivate had all");
  //UImanager->ApplyCommand("/process/eLoss/fluct false");
  //UImanager->ApplyCommand("/process/eLoss/CSDARange true");
  //UImanager->ApplyCommand("/run/verbose 1");
  //UImanager->ApplyCommand("/event/verbose 1");
  //UImanager->ApplyCommand("/tracking/verbose 2");
  //G4UIExecutive * ui = new G4UIExecutive(argc,argv);

  #ifdef VIS

  G4VisManager* visManager = new G4VisExecutive;
  visManager->SetVerboseLevel(0);
  visManager->Initialize();
  G4UImanager * UImanager = G4UImanager::GetUIpointer();
  G4cout << " UI session starts ..." << G4endl;
  G4UIExecutive* ui = new G4UIExecutive(argc, argv);
  UImanager->ApplyCommand("/control/execute vis.mac");
  ui->SessionStart();
  #endif

  runManager->BeamOn( nProtons );
  theAnalysis->Save();
  //calcRSP(theGenerator);
  calcStoppingPower(theGenerator);
  //delete visManager;
  return 0;
  delete runManager;
  }

void calcRSP(PrimaryGeneratorAction* theGenerator){
  G4ParticleDefinition* particle = theGenerator->particle;
  G4EmCalculator* emCal = new G4EmCalculator;
  ofstream myfile;
  myfile.open ("RSP.txt");
  myfile<<"Density RelativeElectronDensity IValue RSP Name"<<endl;
  G4MaterialTable *theMaterialTable = G4Material::GetMaterialTable();
  G4double waterelectrondensity = 0;
  for(size_t i =0;i<theMaterialTable->size();i++){

    G4int I = 0;
    G4double tot =0;
    G4Material* water = theMaterialTable->at(0);
    waterelectrondensity = water->GetElectronDensity()/(g/cm3);
    for(int j=1;j<50000;j++){
      G4double dedx_w = emCal->ComputeElectronicDEDX( double(j)/10*MeV,particle,water);
      G4double dedx_b = emCal->ComputeElectronicDEDX( double(j)/10*MeV,particle,theMaterialTable->at(i));
      tot +=dedx_b/dedx_w;
      I+=1;
    }
    G4double RSP = tot/I;
    G4double electrondensity = theMaterialTable->at(i)->GetElectronDensity()/(g/cm3);
    myfile<<theMaterialTable->at(i)->GetDensity()/(g/cm3)<<" "<<electrondensity/waterelectrondensity<<" "<<theMaterialTable->at(i)->GetIonisation()->GetMeanExcitationEnergy()/eV<<" "<<RSP<<" "<<theMaterialTable->at(i)->GetName()<<endl;
    }
  myfile.close();

}

void calcStoppingPower(PrimaryGeneratorAction* theGenerator){
  G4ParticleDefinition* particle = theGenerator->particle;
  G4EmCalculator* emCal = new G4EmCalculator;
  ofstream myfile;
  myfile.open ("Water_Geant4.dat");
  G4MaterialTable *theMaterialTable = G4Material::GetMaterialTable();
  G4Material* water = theMaterialTable->at(0);
  for(int j=1;j<50000;j++){
    G4double dedx_w = emCal->ComputeElectronicDEDX( double(j)/10*MeV,particle,water);
    myfile<<double(j)/10*MeV<<" "<<dedx_w*MeV/mm<<" "<<endl;
  }
  myfile.close();

}