Skip to content
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.

Sign up for GitHub

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

Jump to bottom

[WIP] Updated version of tripple coincidences code, accoring to Pawel remarks - left for reference #24

Open
wants to merge 8 commits into
base: develop
Choose a base branch
from
Open

[WIP] Updated version of tripple coincidences code, accoring to Pawel remarks - left for reference #24

Changes from 4 commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
c497707
Updated version of tripple coincidences code, accoring to Pawel remarks
dominikpanek Aug 31, 2022
d2a572e
New verison moved to the proper place
dominikpanek Aug 31, 2022
0c7e27a
Event analysis code from the develop brnach updated with tripple coin…
dominikpanek Aug 31, 2022
d2d0f2a
Event analysis code from develop branch, updated with tripple coincid…
dominikpanek Aug 31, 2022
21b2fce
Tripple coincidence code updated with the remarks from mateusz. The c…
dominikpanek Oct 11, 2022
a7eec1a
Functions EventAnalysis.cpp and EventAnalysis.h corrected accoring to…
dominikpanek Oct 14, 2022
549e819
Updated version of the tripple coincidence code. Error were fixed. Co…
dominikpanek Nov 6, 2022
cf2b60f
Updated version of tripple coincidence code
dominikpanek Dec 14, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
116 changes: 113 additions & 3 deletions goja/EventAnalysis.cpp
100755 → 100644
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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) {
Copy link
Member

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:

bool all_same_event(const vector<Hit>& hits) {
 unsigned int nhits = hits.size();
 if ( nhits == 0 ) {
  return false;
 }
 const int ref_id = hits.at(0).eventID;
 return std::all_of(hits.cbegin(), hits.cend(), [ref_id](Hit& hit){ return hit.eventID == ref_id; });
}

bool not_phantom_compton(const vector<Hit>& hits) {
 return std::all_of(hits.cbegin(), hits.cend(), [](Hit& hit){ return hit.nPhantomCompton == 0; });
}

bool all_crystal_compton(const vector<Hit>& hits) {
 return std::all_of(hits.cbegin(), hits.cend(), [](Hit& hit){ return hit.nCrystalCompton == 1; });
}

EventType verify_type_of_triple_coincidence(const std::vector<Hit>& hits) {
 if ( all_same_event(hits) ) {
  //true, phantom-scattered and detector-scattered
  if ( not_phantom_compton(hits) ) {
   if ( all_crystal_compton(hits) ) {
    return kTrue;
   }
   return kDetectorScattered;
  }
  return kPhantomScattered;
 }
 //otherwise --> accidental
 return kAccidental;
}

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)


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
//================================================================================
Expand Down Expand Up @@ -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)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

As I understand these variables are not modified by by function, hence:
void analyze_event(const vector<Hit> &hits, const bool hits_are_singles, const bool triple_coincidence)

{

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;
Expand All @@ -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);
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is error fragile usage of the if().
You should use brackets: if (N==MAX_N and N0<=MAX_N0) {print_coincidences(selected_hits);}

}
else {
if (N==MAX_N and N0<=MAX_N0) print_triple_coincidences(selected_hits);
}

if (DEBUG) {
cout.setf(ios::fixed);
cout.precision(1);
Expand Down Expand Up @@ -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
}
}
}
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Missing end of line