leverington/hit_analyse_v2
1
0
Fork 0
Code Issues Pull Requests repo.project_board Releases Wiki Activity

added the missing Rsq code

Browse Source
This commit is contained in:
Blake Leverington 2021-02-10 17:07:29 +01:00
parent bbb9e3289c
commit 8d47a11394
1 changed files with 8 additions and 3 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

11
Scripts_20210119/hit_analyse_v2.c
View File

@ -315,7 +315,7 @@ beamRecon beamreconstruction(bpm_frame_v2 frametoanalyse, double threshold = 30.
TVectorD M2(3); TVectorD M2(3);
vector<double> signal_list; vector<double> signal_list;
vector<double> channel_list; vector<double> channel_list;
channel_list.clear();
SumY = 0.; SumY = 0.;
SumS = 0.; SumS = 0.;
SumT = 0.; SumT = 0.;
@ -358,7 +358,7 @@ beamRecon beamreconstruction(bpm_frame_v2 frametoanalyse, double threshold = 30.
int temp_lastneighbour= -128; int temp_lastneighbour= -128;
for (int i = 0; i< array_length; i++){ for (int i = 0; i< array_length; i++){
if (std::count( masked_channels[frametoanalyse.board_number].begin(), masked_channels[frametoanalyse.board_number].end(), i)) continue; //check masked pixel list, and do not add it to the list of channels for analysis. if (count( masked_channels[frametoanalyse.board_number].begin(), masked_channels[frametoanalyse.board_number].end(), i)) continue; //check masked pixel list, and do not add it to the list of channels for analysis.
if (frametoanalyse.channel_amp[i]>=threshold) { if (frametoanalyse.channel_amp[i]>=threshold) {
// cout << "ch: " << i << endl; // cout << "ch: " << i << endl;
// signal_list.push_back(frametoanalyse.channel_amp[i]); // signal_list.push_back(frametoanalyse.channel_amp[i]);
@ -452,10 +452,15 @@ beamRecon beamreconstruction(bpm_frame_v2 frametoanalyse, double threshold = 30.
for(int k=0; k<vector_length;k++){ for(int k=0; k<vector_length;k++){
b_num += exp(-p*(channel_list[k]-c)*(channel_list[k]-c)) * signal_list[k]; b_num += exp(-p*(channel_list[k]-c)*(channel_list[k]-c)) * signal_list[k];
b_den += exp(-2*p*(channel_list[k]-c)*(channel_list[k]-c)); b_den += exp(-2*p*(channel_list[k]-c)*(channel_list[k]-c));
} }
b = b_num/b_den; b = b_num/b_den;
for(int k=0; k<vector_length;k++){
SumYYM+= (signal_list[k]-MeanY)*(signal_list[k]-MeanY);
SumYYP+= (b*exp(-p*(channel_list[k]-c)*(channel_list[k]-c)) - MeanY )*(b*exp(-p*(channel_list[k]-c)*( channel_list[k]-c)) - MeanY );
}
// cout << "R-squared = " << SumYYP/SumYYM << endl;
beam.Position = -ABC(1)/ ABC(0); beam.Position = -ABC(1)/ ABC(0);
beam.Focus = 2.3548/sqrt(2*p); beam.Focus = 2.3548/sqrt(2*p);
beam.Peak = b; beam.Peak = b;