| Geant4 Cross Reference |
1 //******************************************** 1 //***********************************************************************************************************
2 // GPSPointLoop.C 2 // GPSPointLoop.C
3 // Root command file 3 // Root command file
4 // Type: root GPSPointLoop.C 4 // Type: root GPSPointLoop.C
5 // 5 //
6 // It generates a macro file to run the simula 6 // It generates a macro file to run the simulation.
7 // 7 //
8 // More information is available in UserGuide 8 // More information is available in UserGuide
9 // Created by Z.LI LP2i Bordeaux 2022 9 // Created by Z.LI LP2i Bordeaux 2022
10 //******************************************** 10 //***********************************************************************************************************
11 11
12 // include <stdio.h> 12 // include <stdio.h>
13 // include <string.h> 13 // include <string.h>
14 // include <stdint.h> 14 // include <stdint.h>
15 // include <vector> 15 // include <vector>
16 // include <math.h> 16 // include <math.h>
17 // using namespace std; 17 // using namespace std;
18 18
19 void GPSPointLoop() 19 void GPSPointLoop()
20 { 20 {
21 gSystem->CopyFile("pixe3d_initial.mac", "pix 21 gSystem->CopyFile("pixe3d_initial.mac", "pixe3d.mac", true);
22 FILE* pfile = fopen("pixe3d.mac", "a+"); 22 FILE* pfile = fopen("pixe3d.mac", "a+");
23 23
24 // gSystem->CopyFile("pixe3d_initial.mac" 24 // gSystem->CopyFile("pixe3d_initial.mac", "pixe3d_stim.mac", true);
25 // FILE* pfile = fopen("pixe3d_stim.mac", 25 // FILE* pfile = fopen("pixe3d_stim.mac", "a+");
26 26
27 //****************************************** 27 //***********************************************
28 //***(begin)** Define scan parameters******* 28 //***(begin)** Define scan parameters************
29 //****************************************** 29 //***********************************************
30 //****************************************** 30 //***********************************************
31 int NumberOfProjections = 10; // Define the 31 int NumberOfProjections = 10; // Define the number of Projections from zero to TotalAngleSpan
32 // (last valu 32 // (last value "TotalAngleSpan" is excluded)
33 int NumberOfSlices = 1; // Define the numbe 33 int NumberOfSlices = 1; // Define the number of Slices
34 int NumberOfPixels = 20; // Define the numb 34 int NumberOfPixels = 20; // Define the number of Pixels for square YZ Scan
35 double TotalAngleSpan = 180; // scan angula 35 double TotalAngleSpan = 180; // scan angular range in degrees
36 double ScanSize = 40 * 1.8; // unit um, sca 36 double ScanSize = 40 * 1.8; // unit um, scan size for cube of 40 um
37 // double ScanSize = 42.48*1.8; // 37 // double ScanSize = 42.48*1.8; // unit um, scan size for C.elegans
38 // double ScanSize = 500; // unit u 38 // double ScanSize = 500; // unit um, scan size for GDP
39 double ScanHeight = ScanSize; // Height of 39 double ScanHeight = ScanSize; // Height of the scan, it depends on the need
40 // double ScanHeight = 201.127; //Height o 40 // double ScanHeight = 201.127; //Height of the scan for STIM-T simulatio of C. elegans, for 128
41 // slices 41 // slices
42 int NbParticles = 1000000; 42 int NbParticles = 1000000;
43 double energy = 1.5; // MeV 43 double energy = 1.5; // MeV
44 char typeParticle[10] = "proton"; 44 char typeParticle[10] = "proton";
45 45
46 double PixelWidth = 1. * ScanSize / NumberOf 46 double PixelWidth = 1. * ScanSize / NumberOfPixels; // Width of each pixel
47 double SliceHeight = 1. * ScanHeight / Numbe 47 double SliceHeight = 1. * ScanHeight / NumberOfSlices; // Height of each
48 48 // slice
49 double AngleStep = 49 double AngleStep =
50 1. * TotalAngleSpan / NumberOfProjections; 50 1. * TotalAngleSpan / NumberOfProjections; // angular increment (in degrees)
51 51 // between two consecutive projections
52 // 52 //
53 // The beam position is at the center of eac 53 // The beam position is at the center of each pixel
54 // Starting position of the beam = StartScan 54 // Starting position of the beam = StartScan + 0.5 x PixelWidth
55 // The scan starts from the bottom left of t 55 // The scan starts from the bottom left of the square
56 // 56 //
57 double StartScanXY = -0.5 * ScanSize; 57 double StartScanXY = -0.5 * ScanSize;
58 double StartScanZ = -0.5 * ScanHeight; 58 double StartScanZ = -0.5 * ScanHeight;
59 // double StartScanZ = 0; 59 // double StartScanZ = 0;
60 60
61 bool isInterrupted = false; 61 bool isInterrupted = false;
62 int P_interrupt = 0; // the start of projec 62 int P_interrupt = 0; // the start of projection index to resume a simulation
63 63
64 //****************************************** 64 //***********************************************
65 //***(end)** Define scan parameters********* 65 //***(end)** Define scan parameters**************
66 //****************************************** 66 //***********************************************
67 67
68 //************************************ 68 //************************************
69 //***(begin)** SCAN IMPLEMENTATION *** 69 //***(begin)** SCAN IMPLEMENTATION ***
70 //************************************ 70 //************************************
71 fprintf(pfile, "/tomography/run/scanParamete 71 fprintf(pfile, "/tomography/run/scanParameters %d %d %d\n", NumberOfProjections, NumberOfSlices,
72 NumberOfPixels); 72 NumberOfPixels);
73 fprintf(pfile, "#\n"); 73 fprintf(pfile, "#\n");
74 74
75 if (isInterrupted) { 75 if (isInterrupted) {
76 fprintf(pfile, "/tomography/run/resumeSimu 76 fprintf(pfile, "/tomography/run/resumeSimulation true\n");
77 fprintf(pfile, "/tomography/run/resumeProj 77 fprintf(pfile, "/tomography/run/resumeProjectionIndex %d\n", P_interrupt);
78 fprintf(pfile, "#\n"); 78 fprintf(pfile, "#\n");
79 } 79 }
80 fprintf(pfile, "/run/initialize\n"); 80 fprintf(pfile, "/run/initialize\n");
81 fprintf(pfile, "#\n"); 81 fprintf(pfile, "#\n");
82 // fprintf(pfile, "/run/printProgress 500000\ 82 // fprintf(pfile, "/run/printProgress 500000\n");
83 // fprintf(pfile, "#\n"); 83 // fprintf(pfile, "#\n");
84 fprintf(pfile, "# Source definition : energy 84 fprintf(pfile, "# Source definition : energy, type\n");
85 fprintf(pfile, "#\n"); 85 fprintf(pfile, "#\n");
86 fprintf(pfile, "/gps/energy %.2f MeV\n", ene 86 fprintf(pfile, "/gps/energy %.2f MeV\n", energy);
87 fprintf(pfile, "/gps/particle %s\n", typePar 87 fprintf(pfile, "/gps/particle %s\n", typeParticle);
88 fprintf(pfile, "#\n"); 88 fprintf(pfile, "#\n");
89 fprintf(pfile, "# SOURCE POSITION AND DIRECT 89 fprintf(pfile, "# SOURCE POSITION AND DIRECTION\n");
90 fprintf(pfile, "#\n"); 90 fprintf(pfile, "#\n");
91 for (int projectionIndex = 0; projectionInde 91 for (int projectionIndex = 0; projectionIndex < NumberOfProjections;
92 ++projectionIndex) // projections 92 ++projectionIndex) // projections
93 { 93 {
94 if (isInterrupted) { 94 if (isInterrupted) {
95 if (projectionIndex < P_interrupt) conti 95 if (projectionIndex < P_interrupt) continue;
96 } 96 }
97 97
98 for (int sliceIndex = 0; sliceIndex < Numb 98 for (int sliceIndex = 0; sliceIndex < NumberOfSlices; ++sliceIndex) // slices
99 { 99 {
100 // if(sliceIndex<15) continue; 100 // if(sliceIndex<15) continue;
101 for (int pixelIndex = 0; pixelIndex < Nu 101 for (int pixelIndex = 0; pixelIndex < NumberOfPixels; ++pixelIndex) // pixels
102 { 102 {
103 double px = cos(projectionIndex * Angl 103 double px = cos(projectionIndex * AngleStep * TMath::DegToRad()); // beam direction
104 double py = sin(projectionIndex * Angl 104 double py = sin(projectionIndex * AngleStep * TMath::DegToRad());
105 double pz = 0.0; 105 double pz = 0.0;
106 double x = 106 double x =
107 StartScanXY * px - (StartScanXY + (p 107 StartScanXY * px - (StartScanXY + (pixelIndex + 0.5) * PixelWidth) * py; // beam position
108 double y = StartScanXY * py + (StartSc 108 double y = StartScanXY * py + (StartScanXY + (pixelIndex + 0.5) * PixelWidth) * px;
109 double z = StartScanZ + (sliceIndex + 109 double z = StartScanZ + (sliceIndex + 0.5) * SliceHeight;
110 // z = 18.07; //if z is fixed 110 // z = 18.07; //if z is fixed
111 // z = z + 1.953125; 111 // z = z + 1.953125;
112 // z = z + 3.90625; 112 // z = z + 3.90625;
113 fprintf(pfile, "/gps/direction %f %f % 113 fprintf(pfile, "/gps/direction %f %f %f\n", px, py, pz);
114 // fprintf(pfile, "/gps/pos/centre %.6 114 // fprintf(pfile, "/gps/pos/centre %.6f %.6f %.6f um\n",x, y, z );
115 fprintf(pfile, "/gps/pos/centre %f %f 115 fprintf(pfile, "/gps/pos/centre %f %f %f um\n", x, y, z);
116 fprintf(pfile, "/run/beamOn %d\n", NbP 116 fprintf(pfile, "/run/beamOn %d\n", NbParticles);
117 fprintf(pfile, "#\n"); 117 fprintf(pfile, "#\n");
118 } 118 }
119 } 119 }
120 } 120 }
121 fclose(pfile); 121 fclose(pfile);
122 122
123 //************************************ 123 //************************************
124 //***(end)** SCAN IMPLEMENTATION *** 124 //***(end)** SCAN IMPLEMENTATION ***
125 //************************************ 125 //************************************
126 } 126 }
127 127