Gabriele Cosmo
297 lines
11 KiB
C++
297 lines
11 KiB
C++
//***********************************************************************************************************
|
|
// BinToStd_gamma_position.C
|
|
// Root command file
|
|
// Type: root BinToStd_gamma_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 PixeEvent_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 PixeEvent
|
|
{
|
|
uint16_t energy_10eV;
|
|
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_gamma_position()
|
|
{
|
|
//***********************************************************************
|
|
//**************************Detection parameters (begin)*****************
|
|
//***********************************************************************
|
|
|
|
const int nbProjection = 1;
|
|
const int nbSlice = 1;
|
|
const int nbPixel = 1;
|
|
double totalAngleSpan = 180.; // in degree
|
|
|
|
double angleOfDetector = 135.; // 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 = 14.726 * TMath::DegToRad(); // in radian
|
|
double radiusOfDetector = distanceObjectDetector * tan(theta);
|
|
bool usePosition = true;
|
|
|
|
//***********************************************************************
|
|
//**************************Detection parameters (end)*******************
|
|
//***********************************************************************
|
|
|
|
FILE* input = fopen("../build/GammaAtExit.dat", "rb");
|
|
FILE* out = fopen("../build/PixeEvent_std_AtExit.DAT", "wb");
|
|
|
|
if (input == NULL) {
|
|
printf("error for opening the input file\n");
|
|
return;
|
|
}
|
|
|
|
RunInfo runInfo;
|
|
PixeEvent pixeEvent;
|
|
Point centerOfDetector;
|
|
Point gammaMomentum;
|
|
Point gammaPosition;
|
|
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> gammaAtExit(nbParticle);
|
|
fread(&gammaAtExit[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) {
|
|
// gamma selection: energy should be lower than 4095*10eV = 49.45 keV
|
|
if (gammaAtExit[i].energy_keV >= 40.95 || gammaAtExit[i].energy_keV <= 0.9) continue;
|
|
|
|
gammaMomentum.m_x = gammaAtExit[i].mx;
|
|
gammaMomentum.m_y = gammaAtExit[i].my;
|
|
gammaMomentum.m_z = gammaAtExit[i].mz;
|
|
|
|
if (!usePosition) {
|
|
if (!IsDetected(centerOfDetector, gammaMomentum, theta)) continue;
|
|
}
|
|
else {
|
|
double c =
|
|
distanceObjectDetector * (gammaMomentum.m_x * cos(ra) + gammaMomentum.m_y * sin(ra));
|
|
if (IsEqual(0, c, eps, releps)) continue; // parallel
|
|
|
|
gammaPosition.m_x = gammaAtExit[i].x;
|
|
gammaPosition.m_y = gammaAtExit[i].y;
|
|
gammaPosition.m_z = gammaAtExit[i].z;
|
|
|
|
double t = (distanceObjectDetector * distanceObjectDetector
|
|
- gammaPosition.m_x * distanceObjectDetector * cos(ra)
|
|
- gammaPosition.m_y * distanceObjectDetector * sin(ra))
|
|
/ c;
|
|
|
|
intersectionPoint.m_x = gammaPosition.m_x + gammaMomentum.m_x * t;
|
|
intersectionPoint.m_y = gammaPosition.m_y + gammaMomentum.m_y * t;
|
|
intersectionPoint.m_z = gammaPosition.m_z + gammaMomentum.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-gammaPosition.m_x*distanceObjectDetector*cos(ra)
|
|
// -gammaPosition.m_y*distanceObjectDetector*sin(ra)), c);
|
|
|
|
// printf(" distanceObjectDetector = %f, gammaPosition.m_x=%f,
|
|
// distanceObjectDetector*cos(ra)=%f, gammaPosition.m_y=%f,
|
|
// distanceObjectDetector*sin(ra)=%f\n", distanceObjectDetector, gammaPosition.m_x,
|
|
// distanceObjectDetector*cos(ra),
|
|
// gammaPosition.m_y,
|
|
// distanceObjectDetector*sin(ra));
|
|
|
|
double tt = (intersectionPoint.m_x - gammaPosition.m_x) * gammaMomentum.m_x
|
|
+ (intersectionPoint.m_y - gammaPosition.m_y) * gammaMomentum.m_y
|
|
+ (intersectionPoint.m_z - gammaPosition.m_z) * gammaMomentum.m_z;
|
|
if (tt < 0) continue;
|
|
}
|
|
|
|
pixeEvent.energy_10eV = floor(100 * gammaAtExit[i].energy_keV + 0.5);
|
|
pixeEvent.projectionIndex = runInfo.projectionIndex;
|
|
pixeEvent.sliceIndex = runInfo.sliceIndex;
|
|
pixeEvent.pixelIndex = runInfo.pixelIndex;
|
|
fwrite(&pixeEvent, 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,
|
|
gammaAtExit[i].mx, gammaAtExit[i].my, gammaAtExit[i].mz, gammaAtExit[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, gammaAtExit[i].mx,
|
|
// gammaAtExit[i].my, gammaAtExit[i].mz, gammaAtExit[i].energy_keV, gammaAtExit[i].x,
|
|
// gammaAtExit[i].y, gammaAtExit[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,
|
|
gammaAtExit[i].mx, gammaAtExit[i].my, gammaAtExit[i].mz, gammaAtExit[i].energy_keV);
|
|
}
|
|
|
|
//***********************************************************************
|
|
//**************************Print information (end)**********************
|
|
//***********************************************************************
|
|
}
|
|
}
|
|
printf(
|
|
"\n---------------Number of PixeEvent in total: "
|
|
"%lld------------------------\n",
|
|
count);
|
|
fclose(input);
|
|
fclose(out);
|
|
|
|
// Recheck the output file in case
|
|
// FILE* input2;
|
|
// input2 = fopen("PixeEvent_std_AtExit.DAT","rb");
|
|
// PixeEvent p;
|
|
// while(fread(&p, 7, 1, input2))
|
|
// {
|
|
// printf("__ProjectionIndex=%d, SliceIndex=%d, PixelIndex=%d,
|
|
// Energy_10eV=%d\n", p.projectionIndex, p.sliceIndex, p.pixelIndex,
|
|
// p.energy_10eV);
|
|
|
|
// }
|
|
// fclose(input2);
|
|
}
|