Gabriele Cosmo
277 lines
10 KiB
C++
277 lines
10 KiB
C++
//***********************************************************************************************************
|
|
// Concatenate_BinToStd_ProtonAtExit.C
|
|
// Root command file
|
|
// Type: root Concatenate_BinToStd_ProtonAtExit.C
|
|
//
|
|
// It is used in case of interruption
|
|
// Read 2 output files ProtonAtExit_1.dat and ProtonAtExit_2.dat that are generated by Geant4
|
|
// tomography simulation It reads protons at exit information, 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;
|
|
|
|
// 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;
|
|
};
|
|
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
|
|
return true;
|
|
}
|
|
|
|
void Recheck()
|
|
{
|
|
// Recheck the output file in case
|
|
FILE* input3 = fopen("../build/StimEvent_std_Detector0_Aperture10.2.DAT", "rb");
|
|
StimEvent p;
|
|
double eventId = -1;
|
|
while (fread(&p, 7, 1, input3)) {
|
|
if (p.projectionIndex == 8 && p.sliceIndex == 64 && p.pixelIndex == 10) {
|
|
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(input3);
|
|
}
|
|
void Concatenate_BinToStd_ProtonAtExit()
|
|
{
|
|
// Recheck();
|
|
// return;
|
|
|
|
//***********************************************************************
|
|
//**************************Detection parameters (begin)*****************
|
|
//***********************************************************************
|
|
|
|
const int nbProjection = 10;
|
|
const int nbSlice = 128;
|
|
const int nbPixel = 20;
|
|
double totalAngleSpan = 180.; // in degree
|
|
|
|
// angle of detector relative to the incident direction of the primary protons at first projection
|
|
// for proton, it is fixed to 0 degree, namely opposite to the source
|
|
double angleOfDetector = 0.;
|
|
double distanceObjectDetector = 22.; // 22 mm
|
|
double radiusOfDetector = 5.; // 5 mm
|
|
// double theta = atan(radiusOfDetector/distanceObjectDetector); //half apex angle of the right
|
|
// circular cone in radian
|
|
double theta = 10.2 * TMath::DegToRad(); // in radian
|
|
|
|
int P_interrupt = 2; // Projection of interruption
|
|
|
|
//***********************************************************************
|
|
//**************************Detection parameters (end)*******************
|
|
//***********************************************************************
|
|
|
|
// assuming there is one interruption
|
|
FILE* input1 = fopen("../build/ProtonAtExit_1.dat", "rb");
|
|
FILE* input2 = fopen("../build/ProtonAtExit_2.dat", "rb");
|
|
FILE* out = fopen("../build/StimEvent_std.DAT", "wb");
|
|
|
|
if (input1 == NULL) {
|
|
printf("error for opening the input ProtonAtExit_1.dat file\n");
|
|
return;
|
|
}
|
|
if (input2 == NULL) {
|
|
printf("error for opening the input ProtonAtExit_2.dat file\n");
|
|
return;
|
|
}
|
|
|
|
RunInfo runInfo;
|
|
StimEvent stimEvent;
|
|
Point centerOfDetector;
|
|
Point protonMomentum;
|
|
|
|
long long count1 = 0;
|
|
long long count2 = 0;
|
|
int runID = -1; // index of simulations, namely runID, starting from 0
|
|
|
|
// ************************************************************(begin)
|
|
// **********************READ FIRST FILE***********************
|
|
// ************************************************************
|
|
while (fread(&runInfo, sizeof(RunInfo), 1, input1)) {
|
|
runID++;
|
|
runInfo.projectionIndex = runID / (nbSlice * nbPixel);
|
|
int remain = runID % (nbSlice * nbPixel);
|
|
runInfo.sliceIndex = remain / nbPixel;
|
|
runInfo.pixelIndex = remain % nbPixel;
|
|
if (runInfo.projectionIndex == P_interrupt) {
|
|
runID--;
|
|
break;
|
|
}
|
|
|
|
int nbParticle = runInfo.nbParticle;
|
|
|
|
//***********************************************************************
|
|
//**************************Print information (begin)********************
|
|
//***********************************************************************
|
|
|
|
printf("-1--runId %d, ProjectionIndex=%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, input1);
|
|
|
|
// if(runInfo.sliceIndex!=1) continue;
|
|
// if(runInfo.sliceIndex!=31&&runInfo.sliceIndex!=32) continue;
|
|
// if(runInfo.sliceIndex!=31) continue;
|
|
|
|
// 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; // proton selection
|
|
|
|
protonMomentum.m_x = protonAtExit[i].mx;
|
|
protonMomentum.m_y = protonAtExit[i].my;
|
|
protonMomentum.m_z = protonAtExit[i].mz;
|
|
|
|
if (!IsDetected(centerOfDetector, protonMomentum, theta))
|
|
continue;
|
|
else {
|
|
stimEvent.energy_keV = floor(protonAtExit[i].energy_keV + 0.5);
|
|
stimEvent.projectionIndex = runInfo.projectionIndex;
|
|
stimEvent.sliceIndex = runInfo.sliceIndex;
|
|
stimEvent.pixelIndex = runInfo.pixelIndex;
|
|
fwrite(&stimEvent, 7, 1, out);
|
|
count1++;
|
|
}
|
|
}
|
|
}
|
|
printf("---------------Number of StimEvent in the first file: %lld------------------------\n",
|
|
count1);
|
|
fclose(input1);
|
|
|
|
// ************************************************************
|
|
// **********************READ FIRST FILE (end)*****************
|
|
// ************************************************************
|
|
|
|
// ************************************************************
|
|
// **********************READ SECOND FILE (begin)**************
|
|
// ************************************************************
|
|
while (fread(&runInfo, sizeof(RunInfo), 1, input2)) {
|
|
runID++;
|
|
runInfo.projectionIndex = runID / (nbSlice * nbPixel);
|
|
int remain = runID % (nbSlice * nbPixel);
|
|
runInfo.sliceIndex = remain / nbPixel;
|
|
runInfo.pixelIndex = remain % nbPixel;
|
|
|
|
int nbParticle = runInfo.nbParticle;
|
|
|
|
//***********************************************************************
|
|
//**************************Print information (begin)********************
|
|
//***********************************************************************
|
|
|
|
printf("-2--runId %d, ProjectionIndex=%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, input2);
|
|
|
|
// if(runInfo.sliceIndex!=1) continue;
|
|
// if(runInfo.sliceIndex!=31) continue;
|
|
// if(runInfo.sliceIndex!=31&&runInfo.sliceIndex!=32) continue;
|
|
|
|
// 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; // proton selection
|
|
|
|
protonMomentum.m_x = protonAtExit[i].mx;
|
|
protonMomentum.m_y = protonAtExit[i].my;
|
|
protonMomentum.m_z = protonAtExit[i].mz;
|
|
|
|
if (!IsDetected(centerOfDetector, protonMomentum, theta))
|
|
continue;
|
|
else {
|
|
stimEvent.energy_keV = floor(protonAtExit[i].energy_keV + 0.5);
|
|
stimEvent.projectionIndex = runInfo.projectionIndex;
|
|
stimEvent.sliceIndex = runInfo.sliceIndex;
|
|
stimEvent.pixelIndex = runInfo.pixelIndex;
|
|
fwrite(&stimEvent, 7, 1, out);
|
|
count2++;
|
|
}
|
|
}
|
|
}
|
|
printf("---------------Number of StimEvent in in the second file: %lld------------------------\n",
|
|
count2);
|
|
|
|
// ************************************************************
|
|
// **********************READ SECOND FILE (end)****************
|
|
// ************************************************************
|
|
|
|
printf("---------------Number of StimEvent in total: %lld------------------------\n",
|
|
count1 + count2);
|
|
fclose(input2);
|
|
fclose(out);
|
|
}
|