Skip to content
Snippets Groups Projects
Commit f60c174d authored by mingf2's avatar mingf2
Browse files

postprocess script

parent ae698014
No related branches found
No related tags found
No related merge requests found
Expand all Hide whitespace changes
Inline Side-by-side
Showing
with 878 additions and 0 deletions
Postprocessing/.gitignore 0 → 100644
+ 1
− 0
View file @ f60c174d
polimi
\ No newline at end of file
Postprocessing/CMakeLists.txt 0 → 100644
+ 13
− 0
View file @ f60c174d
cmake_minimum_required(VERSION 3.16.3)
project(Postprocessing CXX)
set(CMAKE_CXX_STANDARD 11)
set(CMAKE_POSITION_INDEPENDENT_CODE ON)
include_directories(${CMAKE_CURRENT_SOURCE_DIR}/Header)
enable_testing()
add_subdirectory(googletest)
add_subdirectory(Sources)
add_subdirectory(Test)
\ No newline at end of file
Postprocessing/Header/reader.h 0 → 100644
+ 115
− 0
View file @ f60c174d
#pragma once
#include <sstream>
#include <string>
#include <vector>
#include <cmath>
class Vector3D
{
public:
double X=0;
double Y=0;
double Z=0;
Vector3D():Vector3D(0, 0, 0) {}
Vector3D(const double x, const double y, const double z): X(x), Y(y), Z(z) {}
Vector3D operator+(const Vector3D& r) const {
return Vector3D(this->X + r.X, this->Y + r.Y, this->Z + r.Z);
}
Vector3D operator-(const Vector3D& r) const {
return Vector3D(this->X - r.X, this->Y - r.Y, this->Z - r.Z);
}
double operator*(const Vector3D& r) const {
return this->X * r.X + this->Y * r.Y + this->Z * r.Z;
}
Vector3D operator*(const double& d) const {
return Vector3D(this->X * d, this->Y * d, this->Z * d);
}
Vector3D operator/(const double& d) const {
return Vector3D(this->X / d, this->Y / d, this->Z / d);
}
};
inline Vector3D operator*(const double& d, const Vector3D& r){
return r * d;
}
inline double getCosAngle(const Vector3D& l, const Vector3D& r) {
return l*r / std::sqrt((l*l) * (r*r));
}
class Cone
{
public:
// interaction site 1 is the apex
Vector3D apex;
// site 2 - site 1 is the axis
Vector3D axis;
// half angle = abs (scattering angle)
// for neutron, 0 < half angle < pi /2
// for gamma, 0 < half angle < pi
double cosHalfAngle;
// initial energy
double E0;
// energy deposited in compoton scattering
double Edpst;
Cone(): Cone(Vector3D(0,0,0), Vector3D(0,0,0), 0, 0, 0) {}
Cone(const Vector3D point, const Vector3D line, const double cosTheta, const double E0_, const double Edpst_):
apex(point), axis(line), cosHalfAngle(cosTheta), E0(E0_), Edpst(Edpst_) {}
Cone(const std::string& record) {
apex = Vector3D(std::stod(record.substr(0, 8)),
std::stod(record.substr(8, 8)),
std::stod(record.substr(16, 8)));
axis = Vector3D(std::stod(record.substr(24, 8)),
std::stod(record.substr(32, 8)),
std::stod(record.substr(40, 8)));
cosHalfAngle = std::stod(record.substr(48, 13));
E0 = std::stod(record.substr(61, 13));
Edpst = std::stod(record.substr(74, 13));
}
};
class Pulse
{
private:
/* data */
public:
ulong histNo=0; // history number
double height=0; // pulse amplitude, MeV
ushort reaction=0; // reaction type, 1 or 3
int cellNo=0; // cell number
double ergIn=0; // energy of incoming photon, MeV
Pulse() {}
Pulse(std::string record) {
histNo = std::stoi(record.substr(56, 8));
height = std::stod(record.substr(77, 9));
reaction = std::stoi(record.substr(88, 4));
cellNo = std::stoi(record.substr(137, 9));
ergIn = std::stod(record.substr(161, 10));
}
};
class Event
{
private:
/* data */
public:
ulong histNo=0; // history number
int cellNo=0; // cell number
double ergDpst=0; // energy of incoming photon, MeV
Vector3D pos=Vector3D(0,0,0);
Event() {}
Event(std::string record) {
histNo = std::stoi(record.substr(0, 10));
cellNo = std::stoi(record.substr(35, 6));
ergDpst = std::stod(record.substr(44, 8));
pos = Vector3D(stod(record.substr(72,6)),
stod(record.substr(80,6)),
stod(record.substr(88,6)));
}
};
Vector3D getPos(const std::vector<Event>& events, const int& cellNo);
\ No newline at end of file
Postprocessing/Sources/CMakeLists.txt 0 → 100644
+ 4
− 0
View file @ f60c174d
add_library(reader SHARED reader.cpp)
add_executable(main main.cpp)
target_link_libraries(main PUBLIC reader)
\ No newline at end of file
Postprocessing/Sources/main.cpp 0 → 100644
+ 175
− 0
View file @ f60c174d
#include <iostream>
#include <fstream>
#include <vector>
#include <string>
#include <algorithm>
#include <iomanip>
#include <chrono>
#include <set>
#include "reader.h"
int main(int argc, char** argv)
{
// std::string dumnFile("/media/ming/DATA/projects/Imager/polimi/dumn1");
// std::string pulseFile("/media/ming/DATA/projects/Imager/polimi/imager_All_pulses.o");
// std::string outpath("/media/ming/DATA/projects/Imager/polimi/cones.txt");
// test run
std::string dumnFile("/media/ming/DATA/projects/Imager/Test/collision.txt");
std::string pulseFile("/media/ming/DATA/projects/Imager/Test/pulses.txt");
std::string outpath("/media/ming/DATA/projects/Imager/Test/cones.txt");
std::set<int> channels = {13, 15, 17, 19, 21,
53, 55, 57, 59, 61,
73, 75, 77, 79, 81,
93, 95, 97, 99, 101,
303, 305, 307, 309, 311,
503, 505, 507, 509, 511};
std::ifstream fdumn;
fdumn.open(dumnFile, std::ios::in);
if (!fdumn.good())
{
throw std::invalid_argument("Cannot open collision file");
}
std::ifstream fpulse;
fpulse.open(pulseFile, std::ios::in);
if (!fpulse.good())
{
throw std::invalid_argument("Cannot open pulse file");
}
std::ofstream outfile;
outfile.open(outpath, std::ios::out);
if (!outfile.good())
{
throw std::invalid_argument("Cannot create file: " + outpath);
}
auto startTime = std::chrono::high_resolution_clock::now();
// reconstruct events
std::vector<Cone> cones;
std::vector<ulong> histryNumbers;
double cosHalfangle;
Vector3D apex;
Vector3D axis;
std::string pulseLine;
std::string dumnLine;
std::string prevDumnLine("");
std::vector<Pulse> pulses; // pulses with same history number
u_long histNo(0);
while (std::getline(fpulse, pulseLine))
{
// parse line
Pulse newPulse(pulseLine);
// if line.histNo == histNo, add this pulse, go to next line
if (newPulse.histNo == histNo || histNo == 0 || pulses.size() == 0)
{
// // pick up channels of interest
// if (channels.find(newPulse.cellNo) != channels.end())
{
// update histNo;
histNo = newPulse.histNo;
pulses.push_back(std::move(newPulse));
continue;
}
}
// else, process vector of pulses;
if (pulses.size() > 1)
{
// sort the pulses based on incoming photon energy
std::sort(pulses.begin(), pulses.end(),
[](const Pulse & a, const Pulse & b) -> bool
{
return a.ergIn > b.ergIn;
});
// remove anomalies
if (std::abs(pulses[0].ergIn - 0.661) < 0.001 && pulses[0].height < pulses[0].ergIn)
{
// compute scattering angle
cosHalfangle = 1 - pulses[0].height * 0.511 / (pulses[0].ergIn * (pulses[0].ergIn - pulses[0].height));
if (std::abs(cosHalfangle) < 1)
{
// discard cost >= 1, i.e., multiple compton scattering
// get position from collision file
std::vector<Event> interactions;
if (prevDumnLine.length() > 0)
{
u_long histNotmp = stoi(prevDumnLine.substr(0, 10));
if (histNotmp == histNo)
{
interactions.push_back(Event(prevDumnLine));
}
}
while (std::getline(fdumn, dumnLine))
{
u_long histNotmp = stoi(dumnLine.substr(0, 10));
if (histNotmp < histNo)
continue;
else if (histNotmp > histNo)
{
// overshoot
prevDumnLine = dumnLine;
// process vector interactions
// find the interaction that deposits most energy in given cell
// apex = site 1
// axis = site 2 -site 1
apex = getPos(interactions, pulses[0].cellNo);
axis = getPos(interactions, pulses[1].cellNo) - apex;
break;
}
else
{
interactions.push_back(Event(dumnLine));
}
}
if (fdumn.eof())
{
break;
}
if (axis * axis > 0)
{
// axis could be a zero vector if interactions occur at the edge
cones.push_back(Cone(apex, axis, cosHalfangle, pulses[0].ergIn, pulses[0].height));
histryNumbers.push_back(histNo);
}
}
}
}
// empty vector pulses;
pulses.clear();
// update histNo;
histNo = newPulse.histNo;
// add current line to vector;
pulses.push_back(Pulse(std::move(newPulse)));
//go to next line
}
fdumn.close();
fpulse.close();
// save cones to text file
outfile << " apex.x apex.y apex.z axis.x axis.y axis.z cosHalfAngle E0 Edpst history_number\n";
for (int i = 0; i < cones.size(); i++)
{
outfile << std::fixed << std::setprecision(2)
<< std::setw(8) << cones[i].apex.X // << " "
<< std::setw(8) << cones[i].apex.Y // << " "
<< std::setw(8) << cones[i].apex.Z // << " "
<< std::setw(8) << cones[i].axis.X // << " "
<< std::setw(8) << cones[i].axis.Y // << " "
<< std::setw(8) << cones[i].axis.Z // << " "
<< std::fixed << std::setprecision(8)
<< std::setw(13) << cones[i].cosHalfAngle
<< std::setw(13) << cones[i].E0
<< std::setw(13) << cones[i].Edpst
<< std::setw(13) << histryNumbers[i] << '\n';
}
outfile.close();
auto endTime = std::chrono::high_resolution_clock::now();
std::cout << std::chrono::duration_cast<std::chrono::milliseconds>(endTime - startTime).count() << "ms" << std::endl;
return 0;
}
\ No newline at end of file
Postprocessing/Sources/reader.cpp 0 → 100644
+ 18
− 0
View file @ f60c174d
#include "reader.h"
Vector3D getPos(const std::vector<Event>& events, const int& cellNo)
{
// the interaction that deposits most energy in the target cell
Vector3D pos;
double maxDpstErg(0);
for (int i = 0; i < events.size(); i++)
{
if (events[i].cellNo == cellNo && events[i].ergDpst > maxDpstErg)
{
maxDpstErg = events[i].ergDpst;
pos = events[i].pos;
}
}
return pos;
}
\ No newline at end of file
Postprocessing/Test/CMakeLists.txt 0 → 100644
+ 8
− 0
View file @ f60c174d
add_executable(readerTest test.cpp)
target_link_libraries(readerTest PUBLIC gtest_main reader)
add_test(
NAME readerTest
COMMAND readerTest
)
\ No newline at end of file
Postprocessing/Test/collision.txt 0 → 100644
+ 200
− 0
View file @ f60c174d
This diff is collapsed.
Postprocessing/Test/cones.txt 0 → 100644
+ 25
− 0
View file @ f60c174d
apex.x apex.y apex.z axis.x axis.y axis.z cosHalfAngle E0 Edpst history_number
1.76 2.62 -0.05 0.20 -0.16 0.12 -0.68048645 0.66100001 0.45273100 48
5.11 3.38 0.52 -1.80 0.44 -0.80 0.06025452 0.66100001 0.36266100 54
6.30 1.52 -0.90 0.88 0.58 1.27 0.51782044 0.66100001 0.25391000 69
4.51 4.26 0.82 0.02 0.19 -0.04 0.97239547 0.66100001 0.02278900 95
1.57 3.96 -0.64 -0.28 0.92 -0.05 0.99492261 0.66100001 0.00431300 97
2.84 4.49 1.03 -0.18 -0.28 -0.15 -0.61549360 0.66100001 0.44706400 100
6.49 2.16 -0.46 -1.37 1.49 0.59 0.51941873 0.66100001 0.25339100 109
3.40 3.81 0.56 -0.70 -0.53 -0.74 -0.00225151 0.66100001 0.37316500 122
5.25 1.97 1.37 0.14 0.69 0.38 0.92122555 0.66100001 0.06112600 164
1.54 3.17 -0.58 -0.14 -0.34 -0.04 -0.70780034 0.66100001 0.45502400 176
4.15 2.07 0.29 -0.88 1.02 -0.98 0.58314659 0.66100001 0.23156100 178
3.42 2.48 -1.00 -0.31 -0.39 -0.27 -0.52020828 0.66100001 0.43817500 230
3.75 1.40 0.16 -1.40 2.66 -0.64 0.84702101 0.66100001 0.10919400 235
6.69 3.43 -0.32 0.73 1.52 -1.39 0.78720127 0.66100001 0.14267600 301
1.43 1.43 0.74 -0.59 2.19 0.48 0.97727465 0.66100001 0.01887600 318
1.77 1.41 -0.38 0.32 0.78 0.01 0.85542739 0.66100001 0.10413900 320
2.42 2.08 1.26 0.04 0.37 0.43 0.69093463 0.66100001 0.18878600 342
5.95 1.91 -1.28 -0.02 0.28 -0.07 0.95282269 0.66100001 0.03801800 346
2.82 2.35 -0.81 1.70 2.96 -0.02 0.84884332 0.66100001 0.10810600 363
0.48 1.60 1.44 0.66 2.84 -0.82 0.79464811 0.66100001 0.13873100 388
5.47 4.61 -0.03 0.46 -0.02 -0.88 -0.78130296 0.66100001 0.46095100 403
2.82 4.66 1.06 0.53 0.84 0.64 0.72948559 0.66100001 0.17134200 411
4.09 1.57 0.12 -0.97 1.72 -0.04 0.85417659 0.66100001 0.10489700 440
5.67 1.50 -0.83 -0.90 0.16 0.27 -0.03122455 0.66100001 0.37778700 464
Postprocessing/Test/pulses.txt 0 → 100644
+ 100
− 0
View file @ f60c174d
13. 0.40000000 2. 21. 0.737400 1. 1. -999. 0.10000000E+01 35 0.66100001
23. 0.40000000 2. 48. 0.215864 3. 1. -999. 0.10000000E+01 55 0.03062000
24. 0.40000000 2. 48. 0.452731 1. 1. -999. 0.10000000E+01 56 0.66100001
32. 0.40000000 2. 51. 0.046661 1. 1. -999. 0.10000000E+01 74 0.62809998
65. 0.40000000 2. 54. 0.362661 1. 1. -999. 0.10000000E+01 307 0.66100001
46. 0.50000000 2. 54. 0.083879 1. 1. -999. 0.10000000E+01 98 0.30280000
35. 0.50000000 2. 54. 0.207987 3. 1. -999. 0.10000000E+01 77 0.20850000
92. 0.40000000 2. 59. 0.103136 1. 1. -999. 0.10000000E+01 604 0.66100001
99. 0.50000000 2. 64. 0.102308 1. 1. -999. 0.10000000E+01 611 0.65880001
21. 0.40000000 2. 67. 0.503015 1. 1. -999. 0.10000000E+01 53 0.66100001
71. 0.40000000 2. 69. 0.253910 1. 1. -999. 0.10000000E+01 403 0.66100001
82. 0.40000000 2. 69. 0.390592 1. 1. -999. 0.10000000E+01 504 0.39390001
72. 0.40000000 2. 86. 0.172977 1. 1. -999. 0.10000000E+01 404 0.50120002
62. 0.50000000 2. 86. 0.341322 3. 1. -999. 0.10000000E+01 304 0.33640000
57. 0.50000000 2. 95. 0.022789 1. 1. -999. 0.10000000E+01 209 0.66100001
58. 0.50000000 2. 95. 0.631979 3. 1. -999. 0.10000000E+01 210 0.02831000
17. 0.50000000 2. 97. 0.004313 1. 1. -999. 0.10000000E+01 39 0.66100001
19. 0.50000000 2. 97. 0.351328 1. 1. -999. 0.10000000E+01 41 0.65170002
20. 0.50000000 2. 97. 0.295631 1. 1. -999. 0.10000000E+01 42 0.47999999
37. 0.50000000 2. 100. 0.207598 3. 2. -999. 0.10000000E+01 79 0.00400000
38. 0.50000000 2. 100. 0.447064 1. 1. -999. 0.10000000E+01 80 0.66100001
83. 0.40000000 2. 109. 0.253391 1. 1. -999. 0.10000000E+01 505 0.66100001
66. 0.50000000 2. 109. 0.117387 1. 1. -999. 0.10000000E+01 308 0.41060001
26. 0.60000000 2. 109. 0.061210 1. 1. -999. 0.10000000E+01 58 0.29769999
27. 0.60000000 2. 109. 0.240437 3. 1. -999. 0.10000000E+01 59 0.23890001
35. 0.50000000 2. 122. 0.286155 3. 1. -999. 0.10000000E+01 77 0.13630000
46. 0.50000000 2. 122. 0.373165 1. 1. -999. 0.10000000E+01 98 0.66100001
20. 0.50000000 2. 135. 0.253535 1. 1. -999. 0.10000000E+01 42 0.66100001
15. 0.40000000 2. 146. 0.523320 1. 1. -999. 0.10000000E+01 37 0.66100001
14. 0.50000000 2. 146. 0.137591 3. 1. -999. 0.10000000E+01 36 0.03062000
95. 0.50000000 2. 147. 0.456155 1. 1. -999. 0.10000000E+01 607 0.66100001
51. 0.40000000 2. 163. 0.425994 1. 1. -999. 0.10000000E+01 203 0.66100001
62. 0.40000000 2. 164. 0.061126 1. 1. -999. 0.10000000E+01 304 0.66100001
64. 0.40000000 2. 164. 0.139104 1. 1. -999. 0.10000000E+01 306 0.60450000
87. 0.50000000 2. 164. 0.489163 3. 1. -999. 0.10000000E+01 509 0.02861000
15. 0.40000000 2. 176. 0.455024 1. 1. -999. 0.10000000E+01 37 0.66100001
14. 0.50000000 2. 176. 0.190237 3. 1. -999. 0.10000000E+01 36 0.19910000
52. 0.40000000 2. 178. 0.231561 1. 1. -999. 0.10000000E+01 204 0.66100001
45. 0.50000000 2. 178. 0.452133 3. 1. -999. 0.10000000E+01 97 0.16550000
21. 0.40000000 2. 189. 0.253628 1. 1. -999. 0.10000000E+01 53 0.45820001
3. 0.40000000 2. 223. 0.167791 1. 1. -999. 0.10000000E+01 15 0.66100001
49. 0.50000000 2. 227. 0.020326 1. 1. -999. 0.10000000E+01 101 0.66100001
32. 0.40000000 2. 230. 0.207351 3. 1. -999. 0.10000000E+01 74 0.03098000
43. 0.40000000 2. 230. 0.438175 1. 1. -999. 0.10000000E+01 95 0.66100001
63. 0.40000000 2. 231. 0.311037 1. 1. -999. 0.10000000E+01 305 0.66100001
31. 0.40000000 2. 233. 0.091918 1. 1. -999. 0.10000000E+01 73 0.30570000
33. 0.40000000 2. 233. 0.013908 1. 1. -999. 0.10000000E+01 75 0.21089999
41. 0.40000000 2. 235. 0.109194 1. 1. -999. 0.10000000E+01 93 0.66100001
27. 0.50000000 2. 235. 0.553914 3. 1. -999. 0.10000000E+01 59 0.55779999
31. 0.40000000 2. 252. 0.407080 3. 1. -999. 0.10000000E+01 73 0.41530001
24. 0.40000000 2. 255. 0.683741 3. 1. -999. 0.10000000E+01 56 0.66100001
86. 0.50000000 2. 264. 0.701056 3. 1. -999. 0.10000000E+01 508 0.66100001
1. 0.40000000 2. 273. 0.363796 1. 1. -999. 0.10000000E+01 13 0.62970001
98. 0.50000000 2. 279. 0.183004 1. 1. -999. 0.10000000E+01 610 0.66100001
65. 0.40000000 2. 293. 0.289512 1. 1. -999. 0.10000000E+01 307 0.66100001
66. 0.50000000 2. 295. 0.671339 3. 1. -999. 0.10000000E+01 308 0.66100001
85. 0.50000000 2. 301. 0.142676 1. 1. -999. 0.10000000E+01 507 0.66100001
99. 0.50000000 2. 301. 0.069063 1. 1. -999. 0.10000000E+01 611 0.50779998
11. 0.40000000 2. 318. 0.018876 1. 1. -999. 0.10000000E+01 33 0.66100001
15. 0.50000000 2. 318. 0.021255 1. 2. -999. 0.10000000E+01 37 0.17160000
14. 0.50000000 2. 318. 0.126043 3. 2. -999. 0.10000000E+01 36 0.15040000
5. 0.50000000 2. 318. 0.324207 3. 1. -999. 0.10000000E+01 17 0.03098000
16. 0.50000000 2. 318. 0.140536 1. 2. -999. 0.10000000E+01 38 0.30829999
23. 0.40000000 2. 320. 0.363872 1. 1. -999. 0.10000000E+01 55 0.55489999
21. 0.40000000 2. 320. 0.104139 1. 1. -999. 0.10000000E+01 53 0.66100001
11. 0.50000000 2. 320. 0.183944 3. 1. -999. 0.10000000E+01 33 0.18380000
71. 0.40000000 2. 329. 0.669109 1. 1. -999. 0.10000000E+01 403 0.66100001
14. 0.40000000 2. 334. 0.273335 1. 1. -999. 0.10000000E+01 36 0.66100001
12. 0.40000000 2. 339. 0.186464 3. 1. -999. 0.10000000E+01 34 0.18449999
13. 0.40000000 2. 339. 0.418096 1. 1. -999. 0.10000000E+01 35 0.60009998
32. 0.40000000 2. 342. 0.188786 1. 1. -999. 0.10000000E+01 74 0.66100001
33. 0.40000000 2. 342. 0.451286 3. 1. -999. 0.10000000E+01 75 0.47990000
72. 0.40000000 2. 346. 0.038018 1. 1. -999. 0.10000000E+01 404 0.66100001
73. 0.40000000 2. 346. 0.605070 3. 1. -999. 0.10000000E+01 405 0.62220001
2. 0.40000000 2. 349. 0.649951 3. 1. -999. 0.10000000E+01 14 0.66100001
45. 0.40000000 2. 351. 0.435611 1. 1. -999. 0.10000000E+01 97 0.66100001
71. 0.40000000 2. 360. 0.570391 1. 1. -999. 0.10000000E+01 403 0.59140003
72. 0.40000000 2. 360. 0.040621 3. 1. -999. 0.10000000E+01 404 0.03228000
33. 0.40000000 2. 363. 0.108106 1. 1. -999. 0.10000000E+01 75 0.66100001
60. 0.50000000 2. 363. 0.139398 1. 1. -999. 0.10000000E+01 212 0.54360002
1. 0.40000000 2. 365. 0.212371 1. 1. -999. 0.10000000E+01 13 0.66100001
65. 0.40000000 2. 371. 0.671448 1. 1. -999. 0.10000000E+01 307 0.66100001
1. 0.40000000 2. 388. 0.138731 1. 1. -999. 0.10000000E+01 13 0.66100001
18. 0.50000000 2. 388. 0.060128 1. 1. -999. 0.10000000E+01 40 0.52999997
68. 0.50000000 2. 403. 0.460951 1. 1. -999. 0.10000000E+01 310 0.66100001
78. 0.50000000 2. 403. 0.207285 3. 1. -999. 0.10000000E+01 410 0.16700000
38. 0.50000000 2. 411. 0.171342 1. 1. -999. 0.10000000E+01 80 0.66100001
49. 0.50000000 2. 411. 0.164831 3. 1. -999. 0.10000000E+01 101 0.17810000
50. 0.50000000 2. 411. 0.306314 1. 1. -999. 0.10000000E+01 102 0.48740000
21. 0.40000000 2. 425. 0.400719 1. 1. -999. 0.10000000E+01 53 0.66100001
91. 0.40000000 2. 426. 0.494278 1. 1. -999. 0.10000000E+01 603 0.66100001
51. 0.40000000 2. 440. 0.104897 1. 1. -999. 0.10000000E+01 203 0.66100001
39. 0.50000000 2. 440. 0.561963 3. 1. -999. 0.10000000E+01 81 0.19310001
35. 0.50000000 2. 440. 0.028253 1. 1. -999. 0.10000000E+01 77 0.56860000
51. 0.40000000 2. 464. 0.280872 3. 1. -999. 0.10000000E+01 203 0.03098000
71. 0.40000000 2. 464. 0.377787 1. 1. -999. 0.10000000E+01 403 0.66100001
12. 0.40000000 2. 468. 0.185954 1. 1. -999. 0.10000000E+01 34 0.56360000
25. 0.50000000 2. 468. 0.149822 1. 1. -999. 0.10000000E+01 57 0.38400000
35. 0.50000000 2. 468. 0.248850 3. 1. -999. 0.10000000E+01 77 0.24680001
14. 0.40000000 2. 470. 0.270426 1. 1. -999. 0.10000000E+01 36 0.66100001
Postprocessing/Test/test.cpp 0 → 100644
+ 116
− 0
View file @ f60c174d
#include <gtest/gtest.h>
#include <string>
#include "reader.h"
TEST(PulseTest, constructor)
{
std::string record(" 13. 0.40000000 2. 21. 0.737400 1. 1. -999. 0.10000000E+01 35 0.66100001");
Pulse pulse(record);
EXPECT_EQ(pulse.histNo, 21);
EXPECT_DOUBLE_EQ(pulse.height, 0.737400);
EXPECT_EQ(pulse.reaction, 1);
EXPECT_EQ(pulse.cellNo, 35);
EXPECT_DOUBLE_EQ(pulse.ergIn, 0.66100001);
record = " 22. 0.40000000 2. 89765585. 0.114611 3. 1. -999. 0.10000000E+01 54 0.11270000";
pulse = Pulse(record);
EXPECT_EQ(pulse.histNo, 89765585);
EXPECT_DOUBLE_EQ(pulse.height, 0.114611);
EXPECT_EQ(pulse.reaction, 3);
EXPECT_EQ(pulse.cellNo, 54);
EXPECT_DOUBLE_EQ(pulse.ergIn, 0.11270000);
}
TEST(EventTest, constructor)
{
std::string record(" 21 1 2 1 55 35 0.228927 0.04 0.84 2.26 -0.16 1.000E+00 0 0 0 6.610E-01");
Event event(record);
EXPECT_EQ(event.histNo, 21);
EXPECT_EQ(event.cellNo, 35);
EXPECT_DOUBLE_EQ(event.ergDpst, 0.228927);
Vector3D pos(0.84, 2.26, -0.16);
EXPECT_DOUBLE_EQ(event.pos.X, pos.X);
EXPECT_DOUBLE_EQ(event.pos.Y, pos.Y);
EXPECT_DOUBLE_EQ(event.pos.Z, pos.Z);
record = " 99999990 1 2 1 55 604 0.202758 0.04 8.00 1.72 -0.13 1.000E+00 0 0 0 6.610E-01";
event = Event(record);
EXPECT_EQ(event.histNo, 99999990);
EXPECT_EQ(event.cellNo, 604);
EXPECT_DOUBLE_EQ(event.ergDpst, 0.202758 );
pos = Vector3D(8.00, 1.72, -0.13);
EXPECT_DOUBLE_EQ(event.pos.X, pos.X);
EXPECT_DOUBLE_EQ(event.pos.Y, pos.Y);
EXPECT_DOUBLE_EQ(event.pos.Z, pos.Z);
}
TEST(ConeTest, constructor)
{
std::string record(" 1.76 2.62 -0.05 0.20 -0.16 0.12 -0.68048645 0.66100001 0.45273100 48");
Cone cone(record);
EXPECT_DOUBLE_EQ(cone.apex.X, 1.76);
EXPECT_DOUBLE_EQ(cone.apex.Y, 2.62);
EXPECT_DOUBLE_EQ(cone.apex.Z, -0.05);
EXPECT_DOUBLE_EQ(cone.axis.X, 0.20);
EXPECT_DOUBLE_EQ(cone.axis.Y, -0.16);
EXPECT_DOUBLE_EQ(cone.axis.Z, 0.12);
EXPECT_DOUBLE_EQ(cone.cosHalfAngle, -0.68048645);
EXPECT_DOUBLE_EQ(cone.E0, 0.66100001);
EXPECT_DOUBLE_EQ(cone.Edpst, 0.45273100);
record = " 1.43 1.43 0.74 -0.59 2.19 0.48 0.97727465 0.66100001 0.01887600 318";
cone = Cone(record);
EXPECT_DOUBLE_EQ(cone.apex.X, 1.43);
EXPECT_DOUBLE_EQ(cone.apex.Y, 1.43);
EXPECT_DOUBLE_EQ(cone.apex.Z, 0.74);
EXPECT_DOUBLE_EQ(cone.axis.X, -0.59);
EXPECT_DOUBLE_EQ(cone.axis.Y, 2.19);
EXPECT_DOUBLE_EQ(cone.axis.Z, 0.48);
EXPECT_DOUBLE_EQ(cone.cosHalfAngle, 0.97727465);
EXPECT_DOUBLE_EQ(cone.E0, 0.66100001);
EXPECT_DOUBLE_EQ(cone.Edpst, 0.01887600);
}
TEST(getPosTest, getPosTest){
std::vector<std::string> records = {
" 48 1 2 1 55 56 0.429851 0.04 1.76 2.62 -0.05 1.000E+00 0 0 0 6.610E-01",
" 48 1 2 1 55 55 0.073002 0.04 1.95 2.45 0.11 1.000E+00 0 1 0 2.096E-01",
" 48 1 2 3 55 55 0.105975 0.04 1.96 2.46 0.07 1.000E+00 0 2 0 1.366E-01",
" 48 1 2 3 55 55 0.030622 0.04 1.97 2.46 0.08 1.000E+00 0 2 0 3.062E-02",
" 48 3 2 3 55 56 0.013385 0.04 1.76 2.62 -0.05 1.000E+00 0 1 0 1.338E-02",
" 48 2 2 3 55 56 0.008165 0.04 1.76 2.62 -0.05 1.000E+00 0 1 0 8.165E-03"};
std::vector<Event> interactions;
for (int i = 0; i < records.size(); i++)
{
interactions.push_back(Event(records[i]));
}
Vector3D apex = getPos(interactions, 56);
EXPECT_DOUBLE_EQ(apex.X, 1.76);
EXPECT_DOUBLE_EQ(apex.Y, 2.62);
EXPECT_DOUBLE_EQ(apex.Z, -0.05);
Vector3D site2 = getPos(interactions, 55);
EXPECT_DOUBLE_EQ(site2.X, 1.96);
EXPECT_DOUBLE_EQ(site2.Y, 2.46);
EXPECT_DOUBLE_EQ(site2.Z, 0.07);
}
TEST(Vector3DTest, operators)
{
Vector3D v(1,2,3);
EXPECT_EQ(v*v, 14);
Vector3D expected(0.1,0.2,0.3);
Vector3D leftMulti = 0.1 * v;
EXPECT_DOUBLE_EQ(expected.X, leftMulti.X);
EXPECT_DOUBLE_EQ(expected.Y, leftMulti.Y);
EXPECT_DOUBLE_EQ(expected.Z, leftMulti.Z);
Vector3D rightMulti = v * 0.1;
EXPECT_DOUBLE_EQ(expected.X, rightMulti.X);
EXPECT_DOUBLE_EQ(expected.Y, rightMulti.Y);
EXPECT_DOUBLE_EQ(expected.Z, rightMulti.Z);
Vector3D division = v / 10;
EXPECT_DOUBLE_EQ(expected.X, division.X);
EXPECT_DOUBLE_EQ(expected.Y, division.Y);
EXPECT_DOUBLE_EQ(expected.Z, division.Z);
}
\ No newline at end of file
Postprocessing/format.txt 0 → 100644
+ 103
− 0
View file @ f60c174d
dum1 file:
history_No pid gamma reac target cell E_deposited time x y z weight / nsca / erg_par_in(MeV)
21 1 2 1 55 35 0.228927 0.04 0.84 2.26 -0.16 1.000E+00 0 0 0 6.610E-01
21 1 2 1 55 35 0.199572 0.04 1.11 2.57 -0.04 1.000E+00 0 1 0 4.220E-01
21 1 2 3 55 35 0.186618 0.04 1.42 2.27 -0.11 1.000E+00 0 2 0 2.224E-01
21 1 2 3 53 35 0.035813 0.04 1.42 2.27 -0.11 1.000E+00 0 2 0 3.581E-02
21 2 2 3 53 35 0.010071 0.04 0.84 2.26 -0.16 1.000E+00 0 1 0 1.007E-02
48 1 2 1 55 56 0.429851 0.04 1.76 2.62 -0.05 1.000E+00 0 0 0 6.610E-01
48 1 2 1 55 55 0.073002 0.04 1.95 2.45 0.11 1.000E+00 0 1 0 2.096E-01
48 1 2 3 55 55 0.105975 0.04 1.96 2.46 0.07 1.000E+00 0 2 0 1.366E-01
48 1 2 3 55 55 0.030622 0.04 1.97 2.46 0.08 1.000E+00 0 2 0 3.062E-02
48 3 2 3 55 56 0.013385 0.04 1.76 2.62 -0.05 1.000E+00 0 1 0 1.338E-02
48 2 2 3 55 56 0.008165 0.04 1.76 2.62 -0.05 1.000E+00 0 1 0 8.165E-03
51 1 2 1 53 74 0.050400 0.04 2.86 1.86 0.28 1.000E+00 0 1 0 6.281E-01
54 1 2 1 55 307 0.358192 0.04 5.11 3.38 0.52 1.000E+00 0 0 0 6.610E-01
54 1 2 1 53 98 0.094290 0.05 3.31 3.82 -0.28 1.000E+00 0 1 0 3.028E-01
54 1 2 3 53 77 0.179909 0.05 2.91 3.41 0.01 1.000E+00 0 2 0 2.085E-01
54 1 2 3 55 77 0.028609 0.05 2.90 3.40 0.00 1.000E+00 0 2 0 2.861E-02
59 1 2 1 53 604 0.108711 0.04 7.30 2.09 1.47 1.000E+00 0 0 0 6.610E-01
64 1 2 1 55 611 0.099624 0.05 7.21 5.02 -0.33 1.000E+00 0 1 0 6.588E-01
67 1 2 1 55 53 0.186788 0.04 2.16 1.40 -0.55 1.000E+00 0 0 0 6.610E-01
67 1 2 1 55 53 0.304438 0.04 2.18 1.67 -0.86 1.000E+00 0 1 0 4.742E-01
99999939 2 2 3 53 14 0.001967 0.04 0.23 1.79 -0.98 1.000E+00 0 1 0 1.967E-03
99999970 1 2 1 55 33 0.038245 0.04 0.83 1.69 -0.24 1.000E+00 0 0 0 6.610E-01
99999970 1 2 3 53 16 0.591504 0.04 0.49 2.65 0.05 1.000E+00 0 1 0 6.228E-01
99999970 1 2 3 55 16 0.028609 0.04 0.50 2.65 0.06 1.000E+00 0 1 0 2.861E-02
99999970 2 2 3 53 16 0.002642 0.04 0.49 2.65 0.05 1.000E+00 0 1 0 2.642E-03
99999988 1 2 1 55 39 0.472426 0.05 1.01 3.91 -1.14 1.000E+00 0 0 0 6.610E-01
99999988 1 2 3 55 37 0.155694 0.05 1.06 3.28 -0.93 1.000E+00 0 1 0 1.867E-01
99999988 1 2 3 53 37 0.030979 0.05 1.08 3.26 -0.92 1.000E+00 0 1 0 3.098E-02
99999988 2 2 3 53 39 0.001901 0.05 1.01 3.91 -1.14 1.000E+00 0 1 0 1.901E-03
99999990 1 2 1 55 604 0.202758 0.04 8.00 1.72 -0.13 1.000E+00 0 0 0 6.610E-01
119449 1 2 1 55 78 0.165062 0.05 2.40 3.58 0.68 1.000E+00 0 0 0 6.610E-01
119449 2 2 3 53 58 0.034704 0.05 2.40 3.58 0.68 1.000E+00 0 1 0 3.470E-02
26. 0.50000000 2. 119449. 0.030964 3. 2. -999. 0.10000000E+01 58 0.03470000
36. 0.50000000 2. 119449. 0.157417 1. 1. -999. 0.10000000E+01 78 0.66100001
All_pulses.o
time(ns) pulse_type history_No pulse_height(MeV) reac. weight cell erg_par_in(MeV)
13. 0.40000000 2. 21. 0.737400 1. 1. -999. 0.10000000E+01 35 0.66100001
23. 0.40000000 2. 48. 0.215864 3. 1. -999. 0.10000000E+01 55 0.03062000
24. 0.40000000 2. 48. 0.452731 1. 1. -999. 0.10000000E+01 56 0.66100001
32. 0.40000000 2. 51. 0.046661 1. 1. -999. 0.10000000E+01 74 0.62809998
65. 0.40000000 2. 54. 0.362661 1. 1. -999. 0.10000000E+01 307 0.66100001
46. 0.50000000 2. 54. 0.083879 1. 1. -999. 0.10000000E+01 98 0.30280000
35. 0.50000000 2. 54. 0.207987 3. 1. -999. 0.10000000E+01 77 0.20850000
92. 0.40000000 2. 59. 0.103136 1. 1. -999. 0.10000000E+01 604 0.66100001
99. 0.50000000 2. 64. 0.102308 1. 1. -999. 0.10000000E+01 611 0.65880001
21. 0.40000000 2. 67. 0.503015 1. 1. -999. 0.10000000E+01 53 0.66100001
71. 0.40000000 2. 69. 0.253910 1. 1. -999. 0.10000000E+01 403 0.66100001
82. 0.40000000 2. 69. 0.390592 1. 1. -999. 0.10000000E+01 504 0.39390001
72. 0.40000000 2. 86. 0.172977 1. 1. -999. 0.10000000E+01 404 0.50120002
62. 0.50000000 2. 86. 0.341322 3. 1. -999. 0.10000000E+01 304 0.33640000
57. 0.50000000 2. 95. 0.022789 1. 1. -999. 0.10000000E+01 209 0.66100001
58. 0.50000000 2. 95. 0.631979 3. 1. -999. 0.10000000E+01 210 0.02831000
17. 0.50000000 2. 97. 0.004313 1. 1. -999. 0.10000000E+01 39 0.66100001
19. 0.50000000 2. 97. 0.351328 1. 1. -999. 0.10000000E+01 41 0.65170002
20. 0.50000000 2. 97. 0.295631 1. 1. -999. 0.10000000E+01 42 0.47999999
37. 0.50000000 2. 100. 0.207598 3. 2. -999. 0.10000000E+01 79 0.00400000
31. 0.40000000 2. 89765558. 0.134197 1. 1. -999. 0.10000000E+01 73 0.27219999
41. 0.40000000 2. 89765558. 0.385267 1. 1. -999. 0.10000000E+01 93 0.66100001
42. 0.40000000 2. 89765558. 0.143121 3. 1. -999. 0.10000000E+01 94 0.13689999
79. 0.50000000 2. 89765559. 0.133292 1. 1. -999. 0.10000000E+01 411 0.13230000
87. 0.50000000 2. 89765559. 0.540202 1. 1. -999. 0.10000000E+01 509 0.66100001
97. 0.50000000 2. 89765569. 0.414940 1. 1. -999. 0.10000000E+01 609 0.66100001
21. 0.40000000 2. 89765585. 0.466388 1. 1. -999. 0.10000000E+01 53 0.60430002
22. 0.40000000 2. 89765585. 0.114611 3. 1. -999. 0.10000000E+01 54 0.11270000
50. 0.50000000 2. 89765590. 0.098971 1. 1. -999. 0.10000000E+01 102 0.66100001
92. 0.40000000 2. 89765608. 0.675097 1. 1. -999. 0.10000000E+01 604 0.66100001
to reconstruct an event, we need first and second interaction sites, and scattering angle
1. assume initial energy is known as 0.661MeV. we only need to know the energy deposited in the first cell; thus we can select 1-->3 or 1-->1 event pairs. the scattering angle is:
cos(\theta) = 1 - E_d1 * mc^2 / (0.661 * (0.661-E_d1))
2. if the initial energy is unknown, we select 1-->3 event pairs, i.e., photon deposits all its energt in second cell. the scattering angle is:
cos(\theta) = 1 - E_d1 * mc^2 / ((E_d1+E_d2) * E_d2)
while histNo < max:
get pulses with same history number
if number of pulse != 2:
go to next hist number
else:
parse line
sort pulses by erg_par_in
if interaction of second pulse is not 3:
go to next hist number
calculate scattering angle
get interaction sites from dum1 file, based on history number and cell number:
select all interactions with the same history number
sort by cell number
pick up the one with the mathching cell number and max energy of all
how do we know whether the photon deposits all its energy in the second interaction?
Simulation shows that there are many multiple scattering events.
\ No newline at end of file
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment