Gabriele Cosmo
306 lines
11 KiB
C++
306 lines
11 KiB
C++
//***********************************************************************************************************
|
|
// BinToStd_proton_position.C
|
|
// Root command file
|
|
// Type: root BinToStd_proton_position.C
|
|
//
|
|
// Read the X-ray output file that is generated by Geant4 tomography
|
|
// simulation. It reads gamma information, either at creation, or at exit, and rewrite the events
|
|
// in a binary file StimEvent_std.DAT
|
|
//
|
|
// More information is available in UserGuide
|
|
// Created by Z.LI LP2i Bordeaux 2022
|
|
//***********************************************************************************************************
|
|
|
|
#include <math.h>
|
|
#include <stdint.h>
|
|
#include <stdio.h>
|
|
#include <string.h>
|
|
|
|
#include <vector>
|
|
// using namespace std;
|
|
|
|
bool IsEqual(double a, double b, double eps, double releps)
|
|
{
|
|
if (a == b) {
|
|
return true;
|
|
}
|
|
|
|
if (fabs(a - b) <= releps * fabs(b)) {
|
|
return true;
|
|
}
|
|
|
|
if (fabs(a - b) < eps) {
|
|
return true;
|
|
}
|
|
|
|
return false;
|
|
}
|
|
double eps = 1e-20; // absolut difference
|
|
double releps = 1e-10; // relative difference
|
|
|
|
// Define a structure to read and write each event in the required binary format
|
|
struct StimEvent
|
|
{
|
|
uint16_t energy_keV; // different from Pixe Event, it is in keV
|
|
uint16_t pixelIndex;
|
|
uint16_t sliceIndex;
|
|
uint8_t projectionIndex;
|
|
};
|
|
struct ParticleInfo
|
|
{
|
|
float energy_keV;
|
|
float mx;
|
|
float my;
|
|
float mz;
|
|
float x;
|
|
float y;
|
|
float z;
|
|
};
|
|
struct RunInfo
|
|
{
|
|
// uint_16t
|
|
uint8_t projectionIndex; // 1 byte
|
|
uint16_t sliceIndex; //
|
|
uint16_t pixelIndex;
|
|
uint32_t nbParticle; // 4 bytes int
|
|
};
|
|
struct Point
|
|
{
|
|
double m_x;
|
|
double m_y;
|
|
double m_z;
|
|
};
|
|
|
|
bool IsDetected(Point poi1, Point poi2, double theta)
|
|
{
|
|
double a = (poi1.m_x * poi2.m_x + poi1.m_y * poi2.m_y + poi1.m_z * poi2.m_z)
|
|
/ sqrt(poi1.m_x * poi1.m_x + poi1.m_y * poi1.m_y + poi1.m_z * poi1.m_z)
|
|
/ sqrt(poi2.m_x * poi2.m_x + poi2.m_y * poi2.m_y + poi2.m_z * poi2.m_z);
|
|
if (a > 1.0) a = 1;
|
|
if (a < -1.0) a = -1;
|
|
double r = acos(a);
|
|
if (r > theta)
|
|
return false;
|
|
else {
|
|
// printf(" acos: %f, radius: %f\n", r, theta);
|
|
return true;
|
|
}
|
|
}
|
|
|
|
bool IsDetected_position(Point poi1, Point poi2, double r)
|
|
{
|
|
double a = sqrt((poi1.m_x - poi2.m_x) * (poi1.m_x - poi2.m_x)
|
|
+ (poi1.m_y - poi2.m_y) * (poi1.m_y - poi2.m_y)
|
|
+ (poi1.m_z - poi2.m_z) * (poi1.m_z - poi2.m_z));
|
|
|
|
// if(a <= r) return true;
|
|
if (a > r)
|
|
return false;
|
|
|
|
else {
|
|
// printf(" distance of two points: %f, radius: %f\n", a, r);
|
|
return true;
|
|
}
|
|
}
|
|
|
|
void BinToStd_proton_position()
|
|
{
|
|
// printf("%f %f %f\n", acos(1), acos(-1), acos(0));
|
|
// return;
|
|
|
|
//***********************************************************************
|
|
//**************************Detection parameters (begin)*****************
|
|
//***********************************************************************
|
|
|
|
const int nbProjection = 1;
|
|
const int nbSlice = 1;
|
|
const int nbPixel = 1;
|
|
double totalAngleSpan = 180.; // in degree
|
|
|
|
double angleOfDetector = 0.; // angle of detector relative to the incident
|
|
|
|
double distanceObjectDetector = 22000.; // um
|
|
|
|
// double theta = atan(radiusOfDetector/distanceObjectDetector); //half apex
|
|
// angle of the right circular cone in radian
|
|
double theta = 10.2 * TMath::DegToRad(); // in radian
|
|
double radiusOfDetector = distanceObjectDetector * tan(theta);
|
|
bool usePosition = true;
|
|
|
|
//***********************************************************************
|
|
//**************************Detection parameters (end)*******************
|
|
//***********************************************************************
|
|
|
|
FILE* input = fopen("../build/ProtonAtExit.dat", "rb");
|
|
FILE* out = fopen("../build/StimEvent_std", "wb");
|
|
|
|
if (input == NULL) {
|
|
printf("error for opening the input file\n");
|
|
return;
|
|
}
|
|
|
|
RunInfo runInfo;
|
|
StimEvent stimEvent;
|
|
Point centerOfDetector;
|
|
Point protonMomentum;
|
|
Point protonPosition;
|
|
Point intersectionPoint;
|
|
long long count = 0;
|
|
int runID = -1; // index of simulations, namely runID, starting from 0
|
|
|
|
// while(!feof(input)) //if not the end, read
|
|
while (fread(&runInfo, sizeof(RunInfo), 1, input)) {
|
|
runID++;
|
|
int nbParticle = runInfo.nbParticle;
|
|
|
|
//(begin)****************************************************************
|
|
// the following codes are used only when in the simulation
|
|
// the index of projection, slice and pixel is not
|
|
// correctly configured
|
|
runInfo.projectionIndex = runID / (nbSlice * nbPixel);
|
|
int remain = runID % (nbSlice * nbPixel);
|
|
runInfo.sliceIndex = remain / nbPixel;
|
|
runInfo.pixelIndex = remain % nbPixel;
|
|
|
|
//(end)*******************************************************************
|
|
|
|
//***********************************************************************
|
|
//**************************Print information (begin)********************
|
|
//***********************************************************************
|
|
|
|
printf(
|
|
"---------RunID=%d:\nProjectionIndex=%d, SliceIndex=%d, PixelIndex=%d,"
|
|
"nbParticle = %d\n",
|
|
runID, runInfo.projectionIndex, runInfo.sliceIndex, runInfo.pixelIndex, nbParticle);
|
|
|
|
//***********************************************************************
|
|
//**************************Print information (end)**********************
|
|
//***********************************************************************
|
|
|
|
if (!nbParticle) continue;
|
|
std::vector<ParticleInfo> protonAtExit(nbParticle);
|
|
fread(&protonAtExit[0], sizeof(ParticleInfo), nbParticle, input);
|
|
|
|
// angleOfDetector+totalAngleSpan/nbProjection*runInfo.projectionIndex means
|
|
// the angle between source direction and detector, which should be constant
|
|
// when source is rotating
|
|
double ra = TMath::DegToRad()
|
|
* (angleOfDetector + totalAngleSpan / nbProjection * runInfo.projectionIndex);
|
|
centerOfDetector.m_x = distanceObjectDetector * cos(ra);
|
|
centerOfDetector.m_y = distanceObjectDetector * sin(ra);
|
|
centerOfDetector.m_z = 0;
|
|
|
|
for (int i = 0; i < nbParticle; ++i) {
|
|
// proton selection: energy should be lower than 4095 keV
|
|
if (protonAtExit[i].energy_keV >= 4095) continue;
|
|
|
|
protonMomentum.m_x = protonAtExit[i].mx;
|
|
protonMomentum.m_y = protonAtExit[i].my;
|
|
protonMomentum.m_z = protonAtExit[i].mz;
|
|
|
|
if (!usePosition) {
|
|
if (!IsDetected(centerOfDetector, protonMomentum, theta)) continue;
|
|
}
|
|
else {
|
|
double c =
|
|
distanceObjectDetector * (protonMomentum.m_x * cos(ra) + protonMomentum.m_y * sin(ra));
|
|
if (IsEqual(0, c, eps, releps)) continue; // parallel
|
|
|
|
protonPosition.m_x = protonAtExit[i].x;
|
|
protonPosition.m_y = protonAtExit[i].y;
|
|
protonPosition.m_z = protonAtExit[i].z;
|
|
|
|
double t = (distanceObjectDetector * distanceObjectDetector
|
|
- protonPosition.m_x * distanceObjectDetector * cos(ra)
|
|
- protonPosition.m_y * distanceObjectDetector * sin(ra))
|
|
/ c;
|
|
|
|
intersectionPoint.m_x = protonPosition.m_x + protonMomentum.m_x * t;
|
|
intersectionPoint.m_y = protonPosition.m_y + protonMomentum.m_y * t;
|
|
intersectionPoint.m_z = protonPosition.m_z + protonMomentum.m_z * t;
|
|
|
|
if (!IsDetected_position(centerOfDetector, intersectionPoint, radiusOfDetector)) continue;
|
|
|
|
// printf(" t = %f, intersection point: (%f, %f, %f) centor of detector: (%f, %f, %f)
|
|
// 111=%f, 222=%f \n", t, intersectionPoint.m_x,intersectionPoint.m_y,intersectionPoint.m_z,
|
|
// centerOfDetector.m_x,centerOfDetector.m_y,centerOfDetector.m_z,
|
|
// (distanceObjectDetector*distanceObjectDetector-protonPosition.m_x*distanceObjectDetector*cos(ra)
|
|
// -protonPosition.m_y*distanceObjectDetector*sin(ra)), c);
|
|
|
|
// printf(" distanceObjectDetector = %f, protonPosition.m_x=%f,
|
|
// distanceObjectDetector*cos(ra)=%f, protonPosition.m_y=%f,
|
|
// distanceObjectDetector*sin(ra)=%f\n", distanceObjectDetector, protonPosition.m_x,
|
|
// distanceObjectDetector*cos(ra),
|
|
// protonPosition.m_y,
|
|
// distanceObjectDetector*sin(ra));
|
|
|
|
double tt = (intersectionPoint.m_x - protonPosition.m_x) * protonMomentum.m_x
|
|
+ (intersectionPoint.m_y - protonPosition.m_y) * protonMomentum.m_y
|
|
+ (intersectionPoint.m_z - protonPosition.m_z) * protonMomentum.m_z;
|
|
if (tt < 0) continue;
|
|
}
|
|
|
|
stimEvent.energy_10eV = floor(100 * protonAtExit[i].energy_keV + 0.5);
|
|
stimEvent.projectionIndex = runInfo.projectionIndex;
|
|
stimEvent.sliceIndex = runInfo.sliceIndex;
|
|
stimEvent.pixelIndex = runInfo.pixelIndex;
|
|
fwrite(&stimEvent, 7, 1, out);
|
|
count++;
|
|
|
|
//***********************************************************************
|
|
//**************************Print information
|
|
//(begin)********************
|
|
//***********************************************************************
|
|
|
|
if (!usePosition) {
|
|
printf(
|
|
"---------id = %d, RunID=%d ProjectionIndex=%d, SliceIndex=%d, PixelIndex=%d, momentum: "
|
|
"(%f, %f, %f), energy: %f keV\n",
|
|
i, runID, runInfo.projectionIndex, runInfo.sliceIndex, runInfo.pixelIndex,
|
|
protonAtExit[i].mx, protonAtExit[i].my, protonAtExit[i].mz, protonAtExit[i].energy_keV);
|
|
}
|
|
else {
|
|
// printf("---------id = %d, RunID=%d ProjectionIndex=%d, SliceIndex=%d, PixelIndex=%d,
|
|
// momentum: (%f, %f, %f), energy: %f keV, position: (%f, %f, %f)\n", i, runID,
|
|
// runInfo.projectionIndex, runInfo.sliceIndex, runInfo.pixelIndex, protonAtExit[i].mx,
|
|
// protonAtExit[i].my, protonAtExit[i].mz, protonAtExit[i].energy_keV, protonAtExit[i].x,
|
|
// protonAtExit[i].y, protonAtExit[i].z);
|
|
|
|
printf(
|
|
"---------id = %d, RunID=%d ProjectionIndex=%d, SliceIndex=%d, PixelIndex=%d, momentum: "
|
|
"(%f, %f, %f), energy: %f keV\n",
|
|
i, runID, runInfo.projectionIndex, runInfo.sliceIndex, runInfo.pixelIndex,
|
|
protonAtExit[i].mx, protonAtExit[i].my, protonAtExit[i].mz, protonAtExit[i].energy_keV);
|
|
}
|
|
|
|
//***********************************************************************
|
|
//**************************Print information (end)**********************
|
|
//***********************************************************************
|
|
}
|
|
}
|
|
printf(
|
|
"\n---------------Number of StimEvent in total: "
|
|
"%lld------------------------\n",
|
|
count);
|
|
fclose(input);
|
|
fclose(out);
|
|
|
|
// FILE* input2;
|
|
// input2 = fopen("StimEvent_std.DAT","rb");
|
|
// StimEvent p;
|
|
// double eventId = -1;
|
|
// while(fread(&p, 7, 1, input2))
|
|
// {
|
|
|
|
// if(p.projectionIndex == 8 &&p.sliceIndex ==64 && p.pixelIndex==64)
|
|
// {
|
|
// eventId++;
|
|
// printf("StimEvent_%.0f ProjectionIndex=%d, SliceIndex=%d, PixelIndex=%d, Energy_keV=%d keV\n",
|
|
// eventId, p.projectionIndex, p.sliceIndex, p.pixelIndex, p.energy_keV);
|
|
|
|
// }
|
|
|
|
// }
|
|
// fclose(input2);
|
|
}
|