From 8a667cd5d6eee31f47312e05ad402eeb162aa914 Mon Sep 17 00:00:00 2001
From: redford_s <sophie.redford@psi.ch>
Date: Tue, 5 Feb 2019 17:04:44 +0100
Subject: [PATCH] Better filtering of saturation.

Previously, saturated data points were not used, but steps including saturated steps will then still be biased.
Now, if any point within a step saturates the step is not used.
Results look better both per pixel and overall.
Note this requires full rerun (createHistoFile).
---
 CS_analysis.cpp | 33 +++++++++++++++++++++------------
 1 file changed, 21 insertions(+), 12 deletions(-)

diff --git a/CS_analysis.cpp b/CS_analysis.cpp
index 634ec16..cfd1215 100644
--- a/CS_analysis.cpp
+++ b/CS_analysis.cpp
@@ -205,6 +205,9 @@ int main(int argc, char* argv[]) {
   jungfrauPixelMask *pixelMaskObject = new jungfrauPixelMask();
   bool pixel_mask [NCH];
 
+  bool saturation_mask [NCH];
+  std::fill_n(saturation_mask, NCH, true);
+
   jungfrauPedestal *pedestalObject = new jungfrauPedestal();
   pedestalObject->pedestalSetNFrames(640);
 
@@ -481,14 +484,15 @@ int main(int argc, char* argv[]) {
 	      uint16_t adc = imagedptr[i]&0x3fff;
 	      uint16_t gain = (imagedptr[i]&0xc000) >> 14;
 
-	      if (adc > 0) {
-		gain_histos[i]->Fill(gain);
-		if (gain == 0) {
-		  adc_histos_g0[i]->Fill(adc);
-		} else if (gain == 1) {
-		  adc_histos_g1[i]->Fill(adc);
-		} else if (gain == 3) {
-		  adc_histos_g2[i]->Fill(adc);
+	      gain_histos[i]->Fill(gain);
+	      if (gain == 0) {
+		adc_histos_g0[i]->Fill(adc);
+	      } else if (gain == 1) {
+		adc_histos_g1[i]->Fill(adc);
+	      } else if (gain == 3) {
+		adc_histos_g2[i]->Fill(adc);
+		if (adc == 0) {
+		  saturation_mask[i] = false;
 		}
 	      }
 
@@ -560,9 +564,11 @@ int main(int argc, char* argv[]) {
 	      adcermap_avg_g0[j/640-1]->Fill(i%NC,i/NC,adc_histos_g0[i]->GetMeanError());
 	      adcmap_avg_g1[j/640-1]->Fill(i%NC,i/NC,adc_histos_g1[i]->GetMean());
 	      adcermap_avg_g1[j/640-1]->Fill(i%NC,i/NC,adc_histos_g1[i]->GetMeanError());
-	      adcmap_avg_g2[j/640-1]->Fill(i%NC,i/NC,adc_histos_g2[i]->GetMean());
-	      adcermap_avg_g2[j/640-1]->Fill(i%NC,i/NC,adc_histos_g2[i]->GetMeanError());
-	    
+	      if (saturation_mask[i] == true) {
+		adcmap_avg_g2[j/640-1]->Fill(i%NC,i/NC,adc_histos_g2[i]->GetMean());
+		adcermap_avg_g2[j/640-1]->Fill(i%NC,i/NC,adc_histos_g2[i]->GetMeanError());
+	      }
+	      
 	      if (gain_histos[i]->GetEntries() > 0) {
 		gain_avg->Fill(gain_histos[i]->GetMean());
 	      }
@@ -573,7 +579,9 @@ int main(int argc, char* argv[]) {
 		adc_avg_g1->Fill(adc_histos_g1[i]->GetMean());
 	      }
 	      if (adc_histos_g2[i]->GetEntries() > 0) {
-		adc_avg_g2->Fill(adc_histos_g2[i]->GetMean());
+		if (saturation_mask[i] == true) {
+		  adc_avg_g2->Fill(adc_histos_g2[i]->GetMean());
+		}
 	      }
 
 	      gain_histos[i]->Reset();
@@ -672,6 +680,7 @@ int main(int argc, char* argv[]) {
 	  adc_avg_g0->Reset();
 	  adc_avg_g1->Reset();
 	  adc_avg_g2->Reset();
+	  std::fill_n(saturation_mask, NCH, true);
 	}
       } // end of while