normalize image

.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",

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();

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;

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());

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

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);