diff --git a/Data/README.md b/Data/README.md
new file mode 100644
index 0000000000000000000000000000000000000000..00f2a3bcce167b8670189c7a3c8e8ec2d2c93cde
--- /dev/null
+++ b/Data/README.md
@@ -0,0 +1,14 @@
+# Cones and energy deposited simulated by MCNP
+
+## Files
+- `cones_ideal.txt`: all channels used; exact initial energy used. For best result, use this one.
+- `cones_all_channels.txt`: all channels used; initial energy estimated as the sum of two deposited energies
+- `cones_selected_channels.txt`: channels chosen to match our PCB layout; initial energy estimated as the sum of two deposited energies. For most realistic, use this one.
+
+## Fromat
+- `apex`: x_0, cm
+- `axis`: x_0 - x_1, cm
+- `cosHalfAngle`: cos(Theta)
+- `E1+E2`: E_{D1} + E_{D2}, MeV
+- `E_1`: E_{D1}, MeV
+- `history number`: history number in MCNP simulation. Useful only for debugging.
\ No newline at end of file
diff --git a/Headers/reconstruction.h b/Headers/reconstruction.h
index cf082b952e5cc02afeba00d8504a67bbc5197d47..ada2e0f989697f380a5066ffc17a0df962e64a1a 100644
--- a/Headers/reconstruction.h
+++ b/Headers/reconstruction.h
@@ -12,16 +12,23 @@ public:
     RecoImage(const Setup* config_);
     ~RecoImage();
     std::mutex mMutex;
-    ULong64_t counts=0;
+    ULong64_t counts=0; // total number of cones (events)
 
     void clear();
     bool saveImage(const std::string& fileName);
+    /**
+     * @brief add cones to current image and update.
+     * 
+     * @param first Input. First iterator to the vector of cones.
+     * @param last Input. Last iterator to the vector of cones.
+     * @param normailzed Input. Normalize the image if true.
+     */
     void updateImage(std::vector<Cone>::const_iterator first,
                      std::vector<Cone>::const_iterator last,
                      const bool normailzed=true);
 private:
-    const Setup* config;
-    TH2D* hist;
+    const Setup* config; // Global settings.
+    TH2D* hist; // Image to be updated by worker thread and read by main thread.
     void createHist();
     bool hist2txt(const std::string& fileName);
 
@@ -30,9 +37,6 @@ private:
      *        Implementation based on
      *        https://www.nature.com/articles/s41598-020-58857-z
      *
-     * @param config Input. Global settings.
-     * @param histo Input and output. Image to be updated.
-     * @param counts Input and output. Number of cones (events).
      * @param first Input. First iterator to the vector of cones.
      * @param last Input. Last iterator to the vector of cones.
      * @return 0
@@ -45,9 +49,6 @@ private:
      *        Imeplemention based on
      *        https://www.overleaf.com/read/hjdvrcrjcvpx
      *
-     * @param config Input. Global settings.
-     * @param histo Input and output. Image to be updated.
-     * @param counts Input and output. Number of cones (events).
      * @param first Input. First iterator of the vector of cones.
      * @param last Input. Last iterator of the vector of cones.
      * @return 0
diff --git a/Headers/setup.h b/Headers/setup.h
index 6c7c8f0219d3ef9a26da8dff0f23094f32cae374..ed050a40b796aee5a405a1c1b9903d5054e5b4b7 100644
--- a/Headers/setup.h
+++ b/Headers/setup.h
@@ -107,8 +107,6 @@ public:
     // refresh every 10 events
     const int chuckSize = 10;
 
-    // queue capacity
-    ulong capacity=1024;
     // max number of cones to process
     ulong maxN=10000;
 };
@@ -118,7 +116,7 @@ class Cone
 public:
     // interaction site 1 is the apex
     Vector3D apex;
-    // site 2 - site 1 is the axis
+    // site 1 - site 2 is the axis
     Vector3D axis;
     // half angle = abs (scattering angle)
     // for neutron, 0 < half angle < pi /2
diff --git a/README.md b/README.md
index aad3e7e90bf27eeb77a522a2a7ccb44bbdd65bb9..41e20a46ec90cd780493a5ee142f3cac929af6c9 100644
--- a/README.md
+++ b/README.md
@@ -21,3 +21,4 @@ sudo apt-get install qt5-default qtdeclarative5-dev
 ```
 ## Demo
 <img src="Demo/Demo.gif" width="600" height="400" />
+
diff --git a/Sources/reconstruction.cpp b/Sources/reconstruction.cpp
index fc19ee5454c2af6aa70d93b6e4bf71852f8f1ac7..617a961e357be3921d1a65a605f1c328084c0b7b 100644
--- a/Sources/reconstruction.cpp
+++ b/Sources/reconstruction.cpp
@@ -164,10 +164,6 @@ int RecoImage::addConesNormalized(std::vector<Cone>::const_iterator first,
         {
             for (int j = 0; j < config->phiBins; j++)
             {
-		// if (i==0&&j==0)
-		// {
-		//     std::cout << omp_get_num_threads() << '\n';
-		// }
                 ray = k->apex - config->xbs[i][j];
                 beta = getCosAngle(ray, k->axis);
                 sgmb2 = std::pow((ray/(ray*ray) + k->axis/(k->axis*k->axis)-
@@ -177,10 +173,10 @@ int RecoImage::addConesNormalized(std::vector<Cone>::const_iterator first,
                                 (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
+                summ += probDist[i][j] * config->dtheta * config->dphi * std::cos(config->thetaBinCenters[i]);
                 // uncomment to use integral in the bin
                 // 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