-
Notifications
You must be signed in to change notification settings - Fork 13
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
[WIP] Updated version of tripple coincidences code, accoring to Pawel remarks - left for reference #24
base: develop
Are you sure you want to change the base?
[WIP] Updated version of tripple coincidences code, accoring to Pawel remarks - left for reference #24
Changes from 4 commits
c497707
d2a572e
0c7e27a
d2d0f2a
21b2fce
a7eec1a
549e819
cf2b60f
File filter
Filter by extension
Conversations
Jump to
Diff view
Diff view
There are no files selected for viewing
Original file line number | Diff line number | Diff line change |
---|---|---|
|
@@ -232,6 +232,56 @@ EventType verify_type_of_coincidence_castor(const Hit &h1, const Hit &h2) { | |
return kUnspecified; | ||
} | ||
|
||
EventType verify_type_of_triple_coincidence(const Hit &h1, const Hit &h2, const Hit &h3) { | ||
|
||
EventType t = kUnspecified; | ||
|
||
if (h1.eventID==h2.eventID and h1.eventID==h3.eventID) { //true, phantom-scattered and detector-scattered | ||
if (h1.nPhantomCompton==0 and h2.nPhantomCompton==0 and h3.nPhantomCompton==0) { | ||
if (h1.nCrystalCompton==1 and h2.nCrystalCompton==1 and h3.nCrystalCompton==1) { //true | ||
t = kTrue; | ||
} | ||
else { //detector-scattered | ||
t = kDetectorScattered; | ||
} | ||
} | ||
else { //phantom-scattered | ||
t = kPhantomScattered; | ||
} | ||
} | ||
else { //accidental | ||
t = kAccidental; | ||
} | ||
|
||
return t; | ||
|
||
} | ||
|
||
std::tuple<bool, std::vector<Hit>> check_gamma_type(const vector<Hit> &hits, double prompt_energy_threshold) | ||
{ | ||
std::vector<Hit> triple_hits; | ||
|
||
for (unsigned int i=0; i<hits.size(); i++) { | ||
if (hits[i].edep <= prompt_energy_threshold) triple_hits.push_back(hits[i]); | ||
} | ||
int prompt_multiplicity = 0; | ||
for (unsigned int i=0; i<hits.size(); i++) { | ||
if (hits[i].edep > prompt_energy_threshold) { | ||
triple_hits.push_back(hits[i]); | ||
prompt_multiplicity++; | ||
} | ||
} | ||
if(prompt_multiplicity == 1) { | ||
bool isGood = true; | ||
} | ||
else { | ||
bool isGood = false; | ||
} | ||
return std::make_tuple(isGood, triple_hits); | ||
} | ||
|
||
} | ||
|
||
//================================================================================ | ||
// PRINTING | ||
//================================================================================ | ||
|
@@ -295,14 +345,68 @@ void print_coincidences(const std::vector<Hit>& hits) { | |
|
||
} | ||
|
||
void print_triple_coincidences(const std::vector<Hit>& hits, PROMPT_E_TH) { | ||
|
||
assert(hits.size() ==3); | ||
|
||
std::vector<Hit> triple_hits; | ||
bool isGood; | ||
std::tie(isGood, triple_hits) = check_gamma_type(hits, PROMPT_E_TH); | ||
|
||
if(isGood){ | ||
|
||
cout.setf(ios::fixed); | ||
|
||
Hit h1 = triple_hits[0]; | ||
Hit h2 = triple_hits[1]; | ||
Hit h3 = triple_hits[2]; | ||
|
||
cout.precision(2); | ||
cout << h1.posX/10. << "\t" << h1.posY/10. << "\t" << h1.posZ/10. << "\t"; | ||
cout.precision(1); | ||
cout << h1.time << "\t"; | ||
|
||
cout.precision(2); | ||
cout << h2.posX/10. << "\t" << h2.posY/10. << "\t" << h2.posZ/10. << "\t"; | ||
cout.precision(1); | ||
cout << h2.time << "\t"; | ||
|
||
cout << h1.volumeID << "\t"; | ||
cout << h2.volumeID << "\t"; | ||
|
||
cout.precision(2); | ||
cout << h1.edep << "\t"; | ||
cout << h2.edep << "\t"; | ||
|
||
cout << verify_type_of_triple_coincidence(h1, h2, h3) << "\t"; | ||
|
||
cout.precision(2); | ||
cout << h1.sourcePosX/10. << "\t" << h1.sourcePosY/10. << "\t" << h1.sourcePosZ/10. << "\t"; | ||
cout << h2.sourcePosX/10. << "\t" << h2.sourcePosY/10. << "\t" << h2.sourcePosZ/10. << endl; | ||
|
||
cout.precision(2); | ||
cout << h3.posX/10. << "\t" << h3.posY/10. << "\t" << h3.posZ/10. << "\t"; | ||
cout.precision(1); | ||
cout << h3.time << "\t"; | ||
cout << h3.volumeID << "\t"; | ||
cout.precision(2); | ||
cout << h3.edep << "\t"; | ||
cout.precision(2); | ||
cout << h3.sourcePosX/10. << "\t" << h3.sourcePosY/10. << "\t" << h3.sourcePosZ/10. << "\t"; | ||
|
||
} | ||
|
||
} | ||
|
||
//================================================================================ | ||
// MAIN ANALYSIS FUNCTION | ||
//================================================================================ | ||
|
||
void analyze_event(vector<Hit> &hits, bool hits_are_singles) | ||
void analyze_event(vector<Hit> &hits, bool hits_are_singles, bool triple_coincidence) | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. As I understand these variables are not modified by by function, hence: |
||
{ | ||
|
||
double COMPTON_E_TH = atof(getenv("GOJA_COMPTON_E_TH"))*1e3; | ||
double PROMPT_E_TH = atof(getenv(“GOJA_PROMPT_E_TH”))*1e3; | ||
int MAX_N = int(atof(getenv("GOJA_MAX_N"))); | ||
int MAX_N0 = int(atof(getenv("GOJA_MAX_N0"))); | ||
int N = 0; | ||
|
@@ -326,6 +430,13 @@ void analyze_event(vector<Hit> &hits, bool hits_are_singles) | |
if (N==MAX_N and N0<=MAX_N0) print_coincidences(selected_hits); | ||
//if (N>=2 and N<=MAX_N and N0<=MAX_N0) print_coincidences(selected_hits); | ||
|
||
if (triple_coincidences) { | ||
if (N==MAX_N and N0<=MAX_N0) print_coincidences(selected_hits); | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. This is error fragile usage of the if(). |
||
} | ||
else { | ||
if (N==MAX_N and N0<=MAX_N0) print_triple_coincidences(selected_hits); | ||
} | ||
|
||
if (DEBUG) { | ||
cout.setf(ios::fixed); | ||
cout.precision(1); | ||
|
@@ -399,5 +510,4 @@ void RunTests() | |
assert(single4.eventID == 10); // 10 is a value set manually | ||
assert(single4.volumeID == 9); // 9 is a value set manually | ||
assert(std::abs(single4.edep - 140) < epsilon); // Energy of single == sum of energies of hits: 10+80+50 = 140 | ||
} | ||
} | ||
} | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Missing end of line |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Too many if{}else{} plus this implantation is easy to cause some potential errors.
In my opinion you should replace this function with following set of functions:
You should use these 3 new function in other part of code too.
For example, function
EventType verify_type_of_coincidence(const Hit &h1, const Hit &h2)
can be replaced by
EventType verify_type_of_coincidence(const vector<Hit>& hits)