Geant4 LXR

Cross-Referencing   Geant4
Geant4/examples/advanced/stim_pixe_tomography/GPSPointLoop.C

   //***********************************************************************************************************
   // GPSPointLoop.C
   // Root command file
   // Type: root GPSPointLoop.C
   //
   // It generates a macro file to run the simulation.
   //
   // More information is available in UserGuide
   // Created by Z.LI LP2i Bordeaux 2022
  //***********************************************************************************************************
  
  // include <stdio.h>
  // include <string.h>
  // include <stdint.h>
  // include <vector>
  // include <math.h>
  // using namespace std;
  
  void GPSPointLoop()
  {
    gSystem->CopyFile("pixe3d_initial.mac", "pixe3d.mac", true);
    FILE* pfile = fopen("pixe3d.mac", "a+");
  
    //    gSystem->CopyFile("pixe3d_initial.mac", "pixe3d_stim.mac", true);
    //    FILE* pfile = fopen("pixe3d_stim.mac", "a+");
  
    //***********************************************
    //***(begin)** Define scan parameters************
    //***********************************************
    //***********************************************
    int NumberOfProjections = 10;  // Define the number of Projections from zero to TotalAngleSpan
                                   // (last value "TotalAngleSpan" is excluded)
    int NumberOfSlices = 1;  // Define the number of Slices
    int NumberOfPixels = 20;  // Define the number of Pixels for square YZ Scan
    double TotalAngleSpan = 180;  // scan angular range in degrees
    double ScanSize = 40 * 1.8;  // unit um, scan size for cube of 40 um
    //  double ScanSize   = 42.48*1.8;       // unit um, scan size for C.elegans
    //  double ScanSize   = 500;       // unit um, scan size for GDP
    double ScanHeight = ScanSize;  // Height of the scan, it depends on the need
    // double ScanHeight = 201.127;   //Height of the scan for STIM-T simulatio of C. elegans, for 128
    // slices
    int NbParticles = 1000000;
    double energy = 1.5;  // MeV
    char typeParticle[10] = "proton";
  
    double PixelWidth = 1. * ScanSize / NumberOfPixels;  // Width of each pixel
    double SliceHeight = 1. * ScanHeight / NumberOfSlices;  // Height of each
                                                            // slice
    double AngleStep =
      1. * TotalAngleSpan / NumberOfProjections;  // angular increment (in degrees)
                                                  // between two consecutive projections
    //
    // The beam position is at the center of each pixel
    // Starting position of the beam = StartScan + 0.5 x PixelWidth
    // The scan starts from the bottom left of the square
    //
    double StartScanXY = -0.5 * ScanSize;
    double StartScanZ = -0.5 * ScanHeight;
    // double StartScanZ = 0;
  
    bool isInterrupted = false;
    int P_interrupt = 0;  // the start of projection index to resume a simulation
  
    //***********************************************
    //***(end)** Define scan parameters**************
    //***********************************************
  
    //************************************
    //***(begin)** SCAN IMPLEMENTATION ***
    //************************************
    fprintf(pfile, "/tomography/run/scanParameters %d %d %d\n", NumberOfProjections, NumberOfSlices,
            NumberOfPixels);
    fprintf(pfile, "#\n");
  
    if (isInterrupted) {
      fprintf(pfile, "/tomography/run/resumeSimulation true\n");
      fprintf(pfile, "/tomography/run/resumeProjectionIndex %d\n", P_interrupt);
      fprintf(pfile, "#\n");
    }
    fprintf(pfile, "/run/initialize\n");
    fprintf(pfile, "#\n");
  //  fprintf(pfile, "/run/printProgress 500000\n");
  //  fprintf(pfile, "#\n");
    fprintf(pfile, "# Source definition : energy, type\n");
    fprintf(pfile, "#\n");
    fprintf(pfile, "/gps/energy %.2f MeV\n", energy);
    fprintf(pfile, "/gps/particle %s\n", typeParticle);
    fprintf(pfile, "#\n");
    fprintf(pfile, "# SOURCE POSITION AND DIRECTION\n");
    fprintf(pfile, "#\n");
    for (int projectionIndex = 0; projectionIndex < NumberOfProjections;
         ++projectionIndex)  // projections
    {
      if (isInterrupted) {
        if (projectionIndex < P_interrupt) continue;
      }
  
      for (int sliceIndex = 0; sliceIndex < NumberOfSlices; ++sliceIndex)  // slices
      {
       // if(sliceIndex<15) continue;
       for (int pixelIndex = 0; pixelIndex < NumberOfPixels; ++pixelIndex)  // pixels
       {
         double px = cos(projectionIndex * AngleStep * TMath::DegToRad());  // beam direction
         double py = sin(projectionIndex * AngleStep * TMath::DegToRad());
         double pz = 0.0;
         double x =
           StartScanXY * px - (StartScanXY + (pixelIndex + 0.5) * PixelWidth) * py;  // beam position
         double y = StartScanXY * py + (StartScanXY + (pixelIndex + 0.5) * PixelWidth) * px;
         double z = StartScanZ + (sliceIndex + 0.5) * SliceHeight;
         // z = 18.07;  //if z is fixed
         //  z = z + 1.953125;
         //  z = z + 3.90625;
         fprintf(pfile, "/gps/direction %f %f %f\n", px, py, pz);
         // fprintf(pfile, "/gps/pos/centre %.6f %.6f %.6f um\n",x, y, z );
         fprintf(pfile, "/gps/pos/centre %f %f %f um\n", x, y, z);
         fprintf(pfile, "/run/beamOn %d\n", NbParticles);
         fprintf(pfile, "#\n");
       }
     }
   }
   fclose(pfile);
 
   //************************************
   //***(end)** SCAN IMPLEMENTATION ***
   //************************************
 }