advanced/stim_pixe_tomography/GPSPointLoop.C

0001 //***********************************************************************************************************
0002 // GPSPointLoop.C
0003 // Root command file
0004 // Type: root GPSPointLoop.C
0005 //
0006 // It generates a macro file to run the simulation.
0007 //
0008 // More information is available in UserGuide
0009 // Created by Z.LI LP2i Bordeaux 2022
0010 //***********************************************************************************************************
0011
0012 // include <stdio.h>
0013 // include <string.h>
0014 // include <stdint.h>
0015 // include <vector>
0016 // include <math.h>
0017 // using namespace std;
0018
0019 void GPSPointLoop()
0020 {
0021   gSystem->CopyFile("pixe3d_initial.mac", "pixe3d.mac", true);
0022   FILE* pfile = fopen("pixe3d.mac", "a+");
0023
0024   //    gSystem->CopyFile("pixe3d_initial.mac", "pixe3d_stim.mac", true);
0025   //    FILE* pfile = fopen("pixe3d_stim.mac", "a+");
0026
0027   //***********************************************
0028   //***(begin)** Define scan parameters************
0029   //***********************************************
0030   //***********************************************
0031   int NumberOfProjections = 10;  // Define the number of Projections from zero to TotalAngleSpan
0032                                  // (last value "TotalAngleSpan" is excluded)
0033   int NumberOfSlices = 1;  // Define the number of Slices
0034   int NumberOfPixels = 20;  // Define the number of Pixels for square YZ Scan
0035   double TotalAngleSpan = 180;  // scan angular range in degrees
0036   double ScanSize = 40 * 1.8;  // unit um, scan size for cube of 40 um
0037   //  double ScanSize   = 42.48*1.8;       // unit um, scan size for C.elegans
0038   //  double ScanSize   = 500;       // unit um, scan size for GDP
0039   double ScanHeight = ScanSize;  // Height of the scan, it depends on the need
0040   // double ScanHeight = 201.127;   //Height of the scan for STIM-T simulatio of C. elegans, for 128
0041   // slices
0042   int NbParticles = 1000000;
0043   double energy = 1.5;  // MeV
0044   char typeParticle[10] = "proton";
0045
0046   double PixelWidth = 1. * ScanSize / NumberOfPixels;  // Width of each pixel
0047   double SliceHeight = 1. * ScanHeight / NumberOfSlices;  // Height of each
0048                                                           // slice
0049   double AngleStep =
0050     1. * TotalAngleSpan / NumberOfProjections;  // angular increment (in degrees)
0051                                                 // between two consecutive projections
0052   //
0053   // The beam position is at the center of each pixel
0054   // Starting position of the beam = StartScan + 0.5 x PixelWidth
0055   // The scan starts from the bottom left of the square
0056   //
0057   double StartScanXY = -0.5 * ScanSize;
0058   double StartScanZ = -0.5 * ScanHeight;
0059   // double StartScanZ = 0;
0060
0061   bool isInterrupted = false;
0062   int P_interrupt = 0;  // the start of projection index to resume a simulation
0063
0064   //***********************************************
0065   //***(end)** Define scan parameters**************
0066   //***********************************************
0067
0068   //************************************
0069   //***(begin)** SCAN IMPLEMENTATION ***
0070   //************************************
0071   fprintf(pfile, "/tomography/run/scanParameters %d %d %d\n", NumberOfProjections, NumberOfSlices,
0072           NumberOfPixels);
0073   fprintf(pfile, "#\n");
0074
0075   if (isInterrupted) {
0076     fprintf(pfile, "/tomography/run/resumeSimulation true\n");
0077     fprintf(pfile, "/tomography/run/resumeProjectionIndex %d\n", P_interrupt);
0078     fprintf(pfile, "#\n");
0079   }
0080   fprintf(pfile, "/run/initialize\n");
0081   fprintf(pfile, "#\n");
0082 //  fprintf(pfile, "/run/printProgress 500000\n");
0083 //  fprintf(pfile, "#\n");
0084   fprintf(pfile, "# Source definition : energy, type\n");
0085   fprintf(pfile, "#\n");
0086   fprintf(pfile, "/gps/energy %.2f MeV\n", energy);
0087   fprintf(pfile, "/gps/particle %s\n", typeParticle);
0088   fprintf(pfile, "#\n");
0089   fprintf(pfile, "# SOURCE POSITION AND DIRECTION\n");
0090   fprintf(pfile, "#\n");
0091   for (int projectionIndex = 0; projectionIndex < NumberOfProjections;
0092        ++projectionIndex)  // projections
0093   {
0094     if (isInterrupted) {
0095       if (projectionIndex < P_interrupt) continue;
0096     }
0097
0098     for (int sliceIndex = 0; sliceIndex < NumberOfSlices; ++sliceIndex)  // slices
0099     {
0100       // if(sliceIndex<15) continue;
0101       for (int pixelIndex = 0; pixelIndex < NumberOfPixels; ++pixelIndex)  // pixels
0102       {
0103         double px = cos(projectionIndex * AngleStep * TMath::DegToRad());  // beam direction
0104         double py = sin(projectionIndex * AngleStep * TMath::DegToRad());
0105         double pz = 0.0;
0106         double x =
0107           StartScanXY * px - (StartScanXY + (pixelIndex + 0.5) * PixelWidth) * py;  // beam position
0108         double y = StartScanXY * py + (StartScanXY + (pixelIndex + 0.5) * PixelWidth) * px;
0109         double z = StartScanZ + (sliceIndex + 0.5) * SliceHeight;
0110         // z = 18.07;  //if z is fixed
0111         //  z = z + 1.953125;
0112         //  z = z + 3.90625;
0113         fprintf(pfile, "/gps/direction %f %f %f\n", px, py, pz);
0114         // fprintf(pfile, "/gps/pos/centre %.6f %.6f %.6f um\n",x, y, z );
0115         fprintf(pfile, "/gps/pos/centre %f %f %f um\n", x, y, z);
0116         fprintf(pfile, "/run/beamOn %d\n", NbParticles);
0117         fprintf(pfile, "#\n");
0118       }
0119     }
0120   }
0121   fclose(pfile);
0122
0123   //************************************
0124   //***(end)** SCAN IMPLEMENTATION ***
0125   //************************************
0126 }