diff --git a/.vscode/c_cpp_properties.json b/.vscode/c_cpp_properties.json
index fa67db65a4e01f1df7d0fed706ca76b0419ab0cf..4a3ed5d1fc58f151be123195aa1b3da79524ab93 100644
--- a/.vscode/c_cpp_properties.json
+++ b/.vscode/c_cpp_properties.json
@@ -4,8 +4,8 @@
"name": "Linux",
"includePath": [
"${workspaceFolder}/**",
- "/home/ming/root/include",
- "/home/ming/root/lib"
+ "/home/mingf2/root/include",
+ "/home/mingf2/root/lib"
],
"defines": [],
"compilerPath": "/usr/bin/gcc",
diff --git a/Headers/gui.h b/Headers/gui.h
index 4293573aa3d7fbc53b04768f8d6b34889299ebe4..976052ff65d3ea0f9df3755fc490fbabf0eb324b 100644
--- a/Headers/gui.h
+++ b/Headers/gui.h
@@ -19,6 +19,7 @@
#include <TGComboBox.h>
#include <TGFileDialog.h>
#include <TText.h>
+#include <TPaletteAxis.h>
#include "setup.h"
@@ -63,6 +64,8 @@ public:
MyMainFrame(const Setup* config_, const TGWindow *p,UInt_t w,UInt_t h);
virtual ~MyMainFrame();
void addCone(std::vector<Cone>::const_iterator first, std::vector<Cone>::const_iterator last);
+ void addConeNormalized(std::vector<Cone>::const_iterator first,
+ std::vector<Cone>::const_iterator last);
// slots
void Start();
void Stop();
diff --git a/Headers/setup.h b/Headers/setup.h
index ee3310b55ec8e0266d9a48c4b8fed30f62f7e558..6cb83ebed2cabcc217c046568c896cc38be618f6 100644
--- a/Headers/setup.h
+++ b/Headers/setup.h
@@ -63,7 +63,7 @@ public:
}
Setup(const Vector3D v):
R(std::sqrt(v.X*v.X + v.Y*v.Y + v.Z * v.Z)),
- trueSource(v)
+ trueSource(v)
{
getPixelPos();
pixelateSphere();
@@ -78,10 +78,12 @@ public:
const double dtheta = 180 / double(thetaBins) * M_PI / 180;
const double thetaMin = -90 * M_PI / 180 + dtheta / 2;
const double thetaMax = 90 * M_PI / 180 - dtheta / 2;
+ std::vector<double> thetaBinCenters;
// Azimuthal angle, rad
const double dphi = 2 * M_PI /double(phiBins);
const double phiMin = -M_PI + dphi / 2;
const double phiMax = M_PI - dphi / 2;
+ std::vector<double> phiBinCenters;
// pixelate the projection sphere
std::vector<std::vector<Vector3D>> xbs;
int pixelateSphere();
@@ -98,7 +100,7 @@ public:
const int order=3;
// refresh every 10 events
- const int chuckSize = 100;
+ const int chuckSize = 10;
// queue capacity
ulong capacity=1024;
diff --git a/Sources/gui.cpp b/Sources/gui.cpp
index ada7ddeb17baf93fbcae72fd993284502a4840e9..138b5bea33a322afd6df7bc660f4abad061ae69e 100644
--- a/Sources/gui.cpp
+++ b/Sources/gui.cpp
@@ -233,6 +233,64 @@ void MyMainFrame::addCone(std::vector<Cone>::const_iterator first,
}
+void MyMainFrame::addConeNormalized(std::vector<Cone>::const_iterator first,
+ std::vector<Cone>::const_iterator last)
+{
+ if (first == last)
+ {
+ return;
+ }
+
+ // project cones onto the spherical surface
+ // std::cout << "Projecting cones onto the designated spherical surface..."<<std::endl;
+ double alpha(0);
+ double beta(0);
+ double sgma2(0);
+ double sgmb2(0);
+ double summ(0);
+ Vector3D ray;
+ std::vector<std::vector<double>> probDist(config->thetaBins, std::vector<double>(config->phiBins, 0));
+ for (auto k = first; k < last; k++)
+ {
+ counts++;
+ // update the image
+ alpha = k->cosHalfAngle;
+ sgma2 = std::pow(0.511*k->E0/std::pow((k->E0-k->Edpst),2)*config->sgmE, 2);
+ summ = 0;
+ #pragma omp parallel for private(ray, beta, sgmb2) shared(probDist) reduction(+:summ)
+ for (int i = 0; i < config->thetaBins; i++)
+ {
+ for (int j = 0; j < config->phiBins; j++)
+ {
+ ray = k->apex - config->xbs[i][j];
+ beta = getCosAngle(ray, k->axis);
+ sgmb2 = std::pow((ray/(ray*ray) + k->axis/(k->axis*k->axis)-
+ (ray+k->axis)/(ray*k->axis))*config->sgmpos * beta, 2);
+ sgmb2+= std::pow((ray/(ray*k->axis)-k->axis/(k->axis*k->axis))*config->sgmpos * beta, 2);
+ probDist[i][j] = std::exp(-std::pow((std::pow(beta, config->order)-std::pow(alpha, config->order)), 2)/
+ (2*std::pow(config->order, 2)*
+ (std::pow(alpha,2*config->order-2)*sgma2 +std::pow(beta, 2*config->order-2)*sgmb2)));
+ // use proability density
+ // uncomment to use proability counts (integral)
+ // probDist[i][j] = probDist[i][j] * config->dtheta * config->dphi * std::cos(config->thetaBinCenters[i]);
+ // summ += probDist[i][j];
+ summ += probDist[i][j] * config->dtheta * config->dphi * std::cos(config->thetaBinCenters[i]);
+ }
+ }
+ // #pragma omp parallel for
+ for (int i = 0; i < config->thetaBins; i++)
+ {
+ for (int j = 0; j < config->phiBins; j++)
+ {
+ histo->SetBinContent(j+1, i+1, (histo->GetBinContent(j+1, i+1)*(counts - 1) + probDist[i][j] / summ) / counts);
+ }
+ }
+ }
+ // counts += (last - first);
+ // redraw
+ redraw();
+}
+
void MyMainFrame::redraw() {
std::string strtmp = "Total counts: " + std::to_string(counts);
countsText->SetText(0.7, 0.92, strtmp.c_str());
diff --git a/Sources/main.cpp b/Sources/main.cpp
index 93731ed27efc259e3811a3ca5baf175510399f22..43a228f93d588c11bee6f43895541fd02264b8f2 100644
--- a/Sources/main.cpp
+++ b/Sources/main.cpp
@@ -38,7 +38,8 @@ int main(int argc, char** argv) {
i++;
}
// update image
- myMainFrame->addCone(cones.cbegin(), cones.cbegin() + i);
+ // myMainFrame->addCone(cones.cbegin(), cones.cbegin() + i);
+ myMainFrame->addConeNormalized(cones.cbegin(), cones.cbegin() + i);
if (!unfinished)
{
// stop updating after this final run
diff --git a/Sources/setup.cpp b/Sources/setup.cpp
index ea628cf4bcdc823f703c04e64dc2ad539362d701..18f66453815ce86661fa6f8d7eea75f822f12e38 100644
--- a/Sources/setup.cpp
+++ b/Sources/setup.cpp
@@ -21,12 +21,19 @@ int Setup::getPixelPos()
int Setup::pixelateSphere() {
xbs.resize(thetaBins, std::vector<Vector3D>(phiBins, Vector3D(0, 0, 0)));
+ thetaBinCenters.resize(thetaBins, 0);
+ phiBinCenters.resize(phiBins, 0);
double phi = phiMin;
double theta = thetaMin;
for (int i = 0; i < thetaBins; i++)
{
+ thetaBinCenters[i] = theta;
+ phi = phiMin;
for (int j = 0; j < phiBins; j++)
{
+ if (i==0)
+ phiBinCenters[j] = phi;
+
xbs[i][j].X = R * std::cos(theta) * std::cos(phi);
xbs[i][j].Y = R * std::cos(theta) * std::sin(phi);
xbs[i][j].Z = R * std::sin(theta);