diff --git a/Scripts_20210119/hit_analyse_v2.c b/Scripts_20210119/hit_analyse_v2.c
index 376dd84..30fbdb4 100644
--- a/Scripts_20210119/hit_analyse_v2.c
+++ b/Scripts_20210119/hit_analyse_v2.c
@@ -315,7 +315,7 @@ beamRecon beamreconstruction(bpm_frame_v2 frametoanalyse, double threshold = 30.
   TVectorD M2(3);
   vector<double> signal_list;
   vector<double> channel_list;
-
+  channel_list.clear();
   SumY = 0.;
   SumS = 0.;
   SumT = 0.;
@@ -358,7 +358,7 @@ beamRecon beamreconstruction(bpm_frame_v2 frametoanalyse, double threshold = 30.
   
   int temp_lastneighbour= -128;
   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) {
       //  cout << "ch: " << i << endl;
       //    signal_list.push_back(frametoanalyse.channel_amp[i]);
@@ -452,9 +452,14 @@ beamRecon beamreconstruction(bpm_frame_v2 frametoanalyse, double threshold = 30.
   for(int k=0; k<vector_length;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 = 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.Focus = 2.3548/sqrt(2*p);