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

demo data

parent e58b438a
No related branches found
No related tags found
No related merge requests found
This diff is collapsed.
...@@ -86,12 +86,16 @@ public: ...@@ -86,12 +86,16 @@ public:
std::vector<std::vector<Vector3D>> xbs; std::vector<std::vector<Vector3D>> xbs;
int pixelateSphere(); int pixelateSphere();
// dummy simulation setup // simulation setup
const int nCones= 100; const int nCones= 100;
double trueTheta; double trueTheta;
double truePhi; double truePhi;
Vector3D trueSource; Vector3D trueSource;
int randSource(); int randSource();
const double sgmE = 0.06; // 6% energy resolution
const Vector3D sgmpos = Vector3D(0.8 / 2, 0.43 / 2, 1.0 / 2); // mm
// const Vector3D sgmpos(0, 0, 0); // mm
const int order=3;
// refresh every 10 events // refresh every 10 events
const int chuckSize = 100; const int chuckSize = 100;
......
...@@ -59,22 +59,12 @@ int MyMainFrame::Init(){ ...@@ -59,22 +59,12 @@ int MyMainFrame::Init(){
fClearB->Connect("Clicked()", "MyMainFrame", this, "Clear()"); fClearB->Connect("Clicked()", "MyMainFrame", this, "Clear()");
fStopB = new TGTextButton(fF3, "Stop"); fStopB = new TGTextButton(fF3, "Stop");
fStopB->Connect("Clicked()", "MyMainFrame", this, "Stop()"); fStopB->Connect("Clicked()", "MyMainFrame", this, "Stop()");
// // counts
// fCountsLabel = new TGLabel(fF3, "Total counts");
// fCounts = new TGNumberEntryField(fF3, 1, 0, TGNumberFormat::kNESInteger,
// TGNumberEntryField::kNEANonNegative);
// fCounts->SetAlignment(kTextLeft);
// fCounts->SetLimits(TGNumberEntryField::kNELLimitMin, 0);
// fCounts->SetEditable(kFALSE);
// // fCounts->Resize(100,20);
fL3 = new TGLayoutHints(kLHintsTop | kLHintsLeft, fL3 = new TGLayoutHints(kLHintsTop | kLHintsLeft,
5, 5, 5, 5); 5, 5, 5, 5);
fF3->AddFrame(fStartB, fL3); fF3->AddFrame(fStartB, fL3);
fF3->AddFrame(fClearB, fL3); fF3->AddFrame(fClearB, fL3);
fF3->AddFrame(fStopB, fL3); fF3->AddFrame(fStopB, fL3);
// fF3->AddFrame(fCountsLabel, new TGLayoutHints(kLHintsCenterY | kLHintsLeft, 2, 2, 2, 2));
// fF3->AddFrame(fCounts, fL3);
// embedded canvas // embedded canvas
fF5 = new TGCompositeFrame(fMain, 60, 60, kHorizontalFrame); fF5 = new TGCompositeFrame(fMain, 60, 60, kHorizontalFrame);
...@@ -101,8 +91,6 @@ void MyMainFrame::initHist() { ...@@ -101,8 +91,6 @@ void MyMainFrame::initHist() {
histo = new TH2D("ROI", " ; Azimuth; Elevation", histo = new TH2D("ROI", " ; Azimuth; Elevation",
config->phiBins, -180, 180, config->phiBins, -180, 180,
config->thetaBins, -90, 90); config->thetaBins, -90, 90);
// canvas = new TCanvas("THCanvas","Sinocanvas", 1000, 500);
// canvas->Connect("Closed()", "Image", this, "closed()");
// init image // init image
for (int i = 0; i < config->phiBins; i++) for (int i = 0; i < config->phiBins; i++)
{ {
...@@ -215,16 +203,12 @@ void MyMainFrame::addCone(std::vector<Cone>::const_iterator first, ...@@ -215,16 +203,12 @@ void MyMainFrame::addCone(std::vector<Cone>::const_iterator first,
double beta(0); double beta(0);
double sgma2(0); double sgma2(0);
double sgmb2(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; Vector3D ray;
for (auto k = first; k < last; k++) for (auto k = first; k < last; k++)
{ {
// update the image // update the image
alpha = k->cosHalfAngle; alpha = k->cosHalfAngle;
sgma2 = std::pow(0.511*k->E0/std::pow((k->E0-k->Edpst),2)*sgmE, 2); sgma2 = std::pow(0.511*k->E0/std::pow((k->E0-k->Edpst),2)*config->sgmE, 2);
#pragma omp parallel for private(ray, beta, sgmb2) #pragma omp parallel for private(ray, beta, sgmb2)
for (int i = 0; i < config->thetaBins; i++) for (int i = 0; i < config->thetaBins; i++)
{ {
...@@ -233,14 +217,13 @@ void MyMainFrame::addCone(std::vector<Cone>::const_iterator first, ...@@ -233,14 +217,13 @@ void MyMainFrame::addCone(std::vector<Cone>::const_iterator first,
ray = k->apex - config->xbs[i][j]; ray = k->apex - config->xbs[i][j];
beta = getCosAngle(ray, k->axis); beta = getCosAngle(ray, k->axis);
sgmb2 = std::pow((ray/(ray*ray) + k->axis/(k->axis*k->axis)- sgmb2 = std::pow((ray/(ray*ray) + k->axis/(k->axis*k->axis)-
(ray+k->axis)/(ray*k->axis))*sgmpos * beta, 2); (ray+k->axis)/(ray*k->axis))*config->sgmpos * beta, 2);
sgmb2+= std::pow((ray/(ray*k->axis)-k->axis/(k->axis*k->axis))*sgmpos * beta, 2); sgmb2+= std::pow((ray/(ray*k->axis)-k->axis/(k->axis*k->axis))*config->sgmpos * beta, 2);
histo->SetBinContent(j+1, i+1, histo->SetBinContent(j+1, i+1,
histo->GetBinContent(j+1, i+1) + histo->GetBinContent(j+1, i+1) +
std::exp(-std::pow((std::pow(beta, order)-std::pow(alpha, order)), 2)/ std::exp(-std::pow((std::pow(beta, config->order)-std::pow(alpha, config->order)), 2)/
(2*order*order*(std::pow(alpha,2*order-2)*sgma2 +std::pow(beta, 2*order-2)*sgmb2)) (2*std::pow(config->order, 2)*
) (std::pow(alpha,2*config->order-2)*sgma2 +std::pow(beta, 2*config->order-2)*sgmb2))));
);
} }
} }
} }
......
...@@ -45,8 +45,13 @@ void * Worker::handle(void *ptr) { ...@@ -45,8 +45,13 @@ void * Worker::handle(void *ptr) {
void * Worker::reader(void *ptr) { void * Worker::reader(void *ptr) {
std::ifstream conefile; std::ifstream conefile;
std::string fpath("/media/ming/DATA/projects/Imager/polimi/cones.txt"); std::string fpath("Data/cones.txt");
conefile.open(fpath, std::ios::in); conefile.open(fpath, std::ios::in);
if (!conefile.good())
{
throw std::invalid_argument("Cannot find file: " + fpath);
}
std::string line; std::string line;
// skip header (first line) // skip header (first line)
std::getline(conefile, line); std::getline(conefile, line);
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment