add uncertainty

Headers/setup.h

@@ -23,9 +23,17 @@ public:
     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));
@@ -86,7 +94,7 @@ public:
     int randSource();
 
     // refresh every 10 events
-    const int chuckSize = 10;
+    const int chuckSize = 100;
 
     // queue capacity
     ulong capacity=1024;
@@ -103,16 +111,23 @@ public:
     // for neutron, 0 < half angle < pi /2
     // for gamma,   0 < half angle < pi
     double cosHalfAngle;
-    Cone(): Cone(Vector3D(0,0,0), Vector3D(0,0,0), 0) {}
-    Cone(const Vector3D point, const Vector3D line, const double t): apex(point), axis(line), cosHalfAngle(t) {}
+    // 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, 16)), 
-                        std::stod(record.substr(16, 24)));
-        axis = Vector3D(std::stod(record.substr(24, 32)), 
-                        std::stod(record.substr(32, 40)), 
-                        std::stod(record.substr(40, 48)));
-        cosHalfAngle = stod(record.substr(48, 61));
+                        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));
     }
 };
 

Sources/gui.cpp

@@ -213,19 +213,34 @@ void MyMainFrame::addCone(std::vector<Cone>::const_iterator first,
     // std::cout << "Projecting cones onto the designated spherical surface..."<<std::endl;
     double alpha(0);
     double beta(0);
-    // for (const Cone& event : cones)
+    double sgma2(0);
+    double sgmb2(0);
+    const double sgmE(0.06); // 6% energy resolution
+    const Vector3D sgmpos(0.8 / 2, 0.43 / 2, 1.0 / 2); // mm
+    // const Vector3D sgmpos(0, 0, 0); // mm
+    const int order(3);
+    Vector3D ray;
     for (auto k = first; k < last; k++)    
     {
         // update the image
-        //TODO: add the uncertainty
         alpha = k->cosHalfAngle;
-        #pragma omp parallel for
+        sgma2 = std::pow(0.511*k->E0/std::pow((k->E0-k->Edpst),2)*sgmE, 2);
+        #pragma omp parallel for private(ray, beta, sgmb2)
         for (int i = 0; i < config->thetaBins; i++)
         {
             for (int j = 0; j < config->phiBins; j++)
             {
-                double beta = getCosAngle(k->apex - config->xbs[i][j], k->axis);
-                histo->SetBinContent(j+1, i+1, histo->GetBinContent(j+1, i+1) + std::exp(-(beta-alpha)*(beta-alpha) / 0.002));
+                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))*sgmpos * beta, 2);
+                sgmb2+= std::pow((ray/(ray*k->axis)-k->axis/(k->axis*k->axis))*sgmpos * beta, 2);
+                histo->SetBinContent(j+1, i+1, 
+                    histo->GetBinContent(j+1, i+1) + 
+                    std::exp(-std::pow((std::pow(beta, order)-std::pow(alpha, order)), 2)/
+                             (2*order*order*(std::pow(alpha,2*order-2)*sgma2 +std::pow(beta, 2*order-2)*sgmb2))
+                            )
+                );              
             }
         }
     }

Sources/simulation.cpp

@@ -23,7 +23,7 @@ Cone createCone(const Setup* const config){
     // calculate angle
     Vector3D axis = site2 - site1;
     double cosAngle = getCosAngle(site1 - config->trueSource, axis);
-    return Cone(site1, axis, cosAngle);
+    return Cone(site1, axis, cosAngle, 0, 0);
 }
 
 void * Worker::handle(void *ptr) {