Index: trunk/ArchTOIPipe/ProcWSophya/simoffset.cc
===================================================================
--- trunk/ArchTOIPipe/ProcWSophya/simoffset.cc	(revision 2008)
+++ trunk/ArchTOIPipe/ProcWSophya/simoffset.cc	(revision 2014)
@@ -18,15 +18,31 @@
   : poly((degpol > 1)?degpol:1)
 {
+  SetParams(mwsz, nptfit, degpol);
+  totnscount = 0;
+  totnbblock = 0;
+  bmincut = bmaxcut = -9999.;
+  bgalcut = false;
+  ns_flgcut = ns_bgalcut = 0;
+  SavePolyNTuple();
+}
+
+SimpleOffsetEstimator::~SimpleOffsetEstimator()
+{
+}
+
+void SimpleOffsetEstimator::SetParams(int mwsz, int nptfit, int degpol)
+{
   mWSz = (mwsz > 8) ? mwsz : 8;
   nPtFit = (nptfit > degpol+2) ? nptfit : degpol+2;
   if (nPtFit%2 == 0) nPtFit++;
   degPol = (degpol > 1)?degpol:1;
-  totnscount = 0;
-  totnbblock = 0;
-  SavePolyNTuple();
-}
-
-SimpleOffsetEstimator::~SimpleOffsetEstimator()
-{
+  
+}
+
+void SimpleOffsetEstimator::SetBGalCut(double bmin, double bmax)
+{
+  bgalcut = true;
+  bmincut = bmin;
+  bmaxcut = bmax;
 }
 
@@ -36,6 +52,10 @@
      << " SimpleOffsetEstimator::PrintStatus() - MeanWSize= " << mWSz << " NPtFit=" 
      << nPtFit << " DegPoly=" << degPol << " poly.Degre()=" << poly.Degre() << endl;
+  if (bgalcut) 
+    os << " bGalCut = Yes , bGalMin = " << bmincut << " bGalMax= " << bmaxcut << endl;
+  else os << " bGalCut = No " << endl;
   TOIProcessor::PrintStatus(os);
-  os << " ProcessedSampleCount=" << ProcessedSampleCount() << endl;
+  os << " ProcessedSampleCount=" << ProcessedSampleCount() 
+     << " NS_FlgCut= " << ns_flgcut << " NS_BGalCut= " << ns_bgalcut << endl;
   os << " ------------------------------------------------------ " << endl;  
 }
@@ -45,14 +65,12 @@
   cout << "SimpleOffsetEstimator::init" << endl;
   declareInput("in");
+  declareInput("bgal");
   declareOutput("offset");
   declareOutput("out");
   declareOutput("incopie");
-  declareOutput("poly_a0");
-  declareOutput("poly_a1");
-  declareOutput("poly_a2");
-  declareOutput("poly_sn0");
-  declareOutput("mean_y");
-  declareOutput("sig_y");
-  declareOutput("mean_x");
+  declareOutput("bgalcopie");
+  //  declareOutput("mean_y");
+  //  declareOutput("sig_y");
+  //  declareOutput("mean_x");
   name = "SimpleOffsetEstimator";
 }
@@ -66,16 +84,18 @@
   bool fgout = checkOutputTOIIndex(1);
   bool fgincopie = checkOutputTOIIndex(2);
-  bool fga0 = checkOutputTOIIndex(3);
-  bool fga1 = checkOutputTOIIndex(4);
-  bool fga2 = checkOutputTOIIndex(5);
-  bool fgsn0 = checkOutputTOIIndex(6);
-  bool fgmeany = checkOutputTOIIndex(7);
-  bool fgsigy = checkOutputTOIIndex(8);
-  bool fgmeanx = checkOutputTOIIndex(9);
-  
+  bool fgbgalcopie = checkOutputTOIIndex(3);
+  //  bool fgmeany = checkOutputTOIIndex(4);
+  //  bool fgsigy = checkOutputTOIIndex(5);
+  //  bool fgmeanx = checkOutputTOIIndex(6);
+
   if (!checkInputTOIIndex(0)) {
     cerr << " SimpleOffsetEstimator::run() - Input TOI (in) not connected! "
 	 << endl;
     throw ParmError("SimpleOffsetEstimator::run() Input TOI (in) not connected!");
+  }
+  if (bgalcut && !checkInputTOIIndex(1)) {
+    cerr << " SimpleOffsetEstimator::run() - Input TOI bgal not connected! "
+	 << endl;
+    throw ParmError("SimpleOffsetEstimator::run() Input TOI bgal not connected!");
   }
   if (!fgoffset && !fgout) {
@@ -88,6 +108,7 @@
 
   // NTuple pour sauvegarde des coeff de poly  
-  char * nomsnt[] = {"sncur", "sn0", "meanx", "meany", "sigy", "a0", "a1", "a2", "ycur"};
-  XNTuple xntp(0, 9, 0, 0, nomsnt);
+  char * nomsnt[] = {"sncur", "sn0", "meanx", "meany", "sigy", 
+		       "a0", "a1", "a2", "ycur", "nok", "nbblkok"};
+  XNTuple xntp(0, 11, 0, 0, nomsnt);
 
   try {
@@ -103,4 +124,5 @@
     Vector vinc(wsize);
     TVector<uint_8> vfgc(wsize);
+    Vector bgal(wsize);
 
     
@@ -131,4 +153,9 @@
     double y_last = 0.;
 
+    double last_meanok = 0.;
+    double last_sigok = 1.;
+    double meanx_forlast = 0.;
+    bool fg_last_meansig = false;
+ 
     // Boucle sur les sampleNum
     // 1ere partie, on traite par paquets de wsize
@@ -142,11 +169,25 @@
 
       getData(0, k, wsize, vin.Data(), vfg.Data());
+      if (bgalcut) {
+	getData(1, k, wsize, bgal.Data());
+	if (fgbgalcopie) putData(3, k, wsize, bgal.Data());
+      }
+
+      fg_last_meansig = false;
+
+      if (kb == 0) {
+	last_meanok = vin.Sum()/wsize;
+	last_sigok = vin.SumX2()/wsize - last_meanok*last_meanok;
+      }
 
       // Calcul moyenne et sigma du bloc 
       nok = 0.;  meanx = 0.;
       mean = 0.;  sig = 0.;
-      
+      meanx_forlast = 0.;
       for(j=0; j<wsize; j++) {
-	if ( vfg(j) ) continue;
+	meanx_forlast += k+j;
+	if ( vfg(j) ) { ns_flgcut++;  continue; }
+	if (bgalcut && (bgal(j) > bmincut) && (bgal(j) < bmaxcut) ) 
+	  { ns_bgalcut++;  continue; }
 	mean += vin(j);
 	sig += vin(j)*vin(j);
@@ -172,5 +213,6 @@
 	meanx /= nok;
 	sig = sig/nok-mean*mean;
-	if (sig < 1.e-6*mean) sig = 1.e-6*mean;
+	if (sig < 1.e-10*mean) sig = 1.e-10*mean;
+	if (sig < 1.e-39) sig = 1.e-39;
 	int kk = nbblkok%nPtFit;
 	nbblkok++;
@@ -178,11 +220,21 @@
 	YErr(kk) = sig;
 	X(kk) = meanx;
-      }
-      
-      if ((nok>3.5) && (nbblkok >= nPtFit) ) {
+	last_meanok = mean;
+	last_sigok = sig;
+      }
+      else {
+	int kk = nbblkok%nPtFit;
+	Y(kk) = last_meanok; 
+	YErr(kk) = last_sigok*10.;
+	X(kk) = meanx_forlast/wsize;	
+	fg_last_meansig = true;
+      }
+      
+      if (((nok>3.5) || fg_last_meansig) && (nbblkok >= nPtFit) ) {
 	// On force le fit a garder une certaine continuite
 	Y(nPtFit) = y_last;
 	double smin, smax;
 	YErr(Range(0,0,nPtFit)).MinMax(smin, smax);
+	if (smax < 1.e-39)   smax = 1.e-39;
 	YErr(nPtFit) = smax/10.;
 	X(nPtFit) = sn_last;
@@ -190,5 +242,15 @@
 	X0 = X;
 	X0 -= sn0;
-	poly.Fit(X0,Y,YErr,degPol,errCoef);
+	Vector polsave(degPol);
+	for(int jj=0; jj<=poly.Degre(); jj++) polsave(jj) = poly[jj];
+	try {
+	  poly.Fit(X0,Y,YErr,degPol,errCoef);
+	}
+	catch (CaughtSignalExc& excsig) {
+	  cout << " -- simoffset.cc/ catched CaughtSignalExc - Msg= " 
+	       << excsig.Msg() << " kb=" << kb << " k=" << k << endl;
+	  cout << " X0= " << X0 << " Y=" << Y << " YErr=" << YErr << endl;
+	  for(int jj=0; jj<=poly.Degre(); jj++) poly[jj] = polsave(jj);
+	}
       }
       else {
@@ -208,5 +270,5 @@
       }
       if (ntpoly) {    // Remplissage du XNTuple de controle 
-	float xnt[10];
+	float xnt[12];
 	xnt[0] = k;
 	xnt[1] = sn0;
@@ -218,4 +280,6 @@
 	xnt[7] = poly[2];
 	xnt[8] = poly(k-sn0);
+	xnt[9] = nok;
+	xnt[10] = nbblkok;
 	xntp.Fill(NULL, xnt, NULL, NULL);
       }
@@ -244,4 +308,5 @@
 	vinc.ReSize(wszt); 
 	vfgc.ReSize(wszt);
+	bgal.ReSize(wszt); 
 	int jb = 0;
 	for(int kbsc=kbs; kbsc<nblk; kbsc++) {
@@ -270,34 +335,19 @@
 	putData(1, k, wsize, vinc.Data(), vfgc.Data());
       }
-      
-      if (fga0) {
-	vout = poly[0];
-	putData(3, k, wsize, vout.Data());
-      }
-      if (fga1) {
-	vout = poly[1];
-	putData(4, k, wsize, vout.Data());
-      }
-      if (fga2) {
-	vout = poly[2];
-	putData(5, k, wsize, vout.Data());
-      }
-      if (fgsn0) {
-	vout = sn0;
-	putData(6, k, wsize, vout.Data());
-      }
+      /*
       if (fgmeany) {
 	vout = mean;
-	putData(7, k, wsize, vout.Data());
+	putData(4, k, wsize, vout.Data());
       }
       if (fgsigy) {
 	vout = sig;
-	putData(8, k, wsize, vout.Data());
+	putData(5, k, wsize, vout.Data());
       }
 
       if (fgmeanx) {
 	vout = meanx;
-	putData(9, k, wsize, vout.Data());
-      }
+	putData(6, k, wsize, vout.Data());
+      }
+      */
 
       klast+=wsize;
@@ -316,4 +366,10 @@
       k = klast+1;
       getData(0, k, wsize, vin.Data(), vfg.Data());
+
+      if (bgalcut) {
+	getData(1, k, wsize, bgal.Data());
+	if (fgbgalcopie) putData(3, k, wsize, bgal.Data());
+      }
+
       for(j=0; j<wsize; j++) 
 	voff(j) = poly(k+j-sn0);
@@ -324,34 +380,20 @@
 	putData(1, k, wsize, vin.Data(), vfg.Data());
       }
-      
-      if (fga0) {
-	vout = poly[0];
-	putData(3, k, wsize, vout.Data());
-      }
-      if (fga1) {
-	vout = poly[1];
-	putData(4, k, wsize, vout.Data());
-      }
-      if (fga2) {
-	vout = poly[2];
-	putData(5, k, wsize, vout.Data());
-      }
-      if (fgsn0) {
-	vout = sn0;
-	putData(6, k, wsize, vout.Data());
-      }
+
+      /*      
       if (fgmeany) {
 	vout = mean;
-	putData(7, k, wsize, vout.Data());
+	putData(4, k, wsize, vout.Data());
       }
       if (fgsigy) {
 	vout = sig;
-	putData(8, k, wsize, vout.Data());
+	putData(5, k, wsize, vout.Data());
       }
       
       if (fgmeanx) {
 	vout = meanx;
-	putData(9, k, wsize, vout.Data());
-      }
+	putData(6, k, wsize, vout.Data());
+      }
+      */
       
       klast+=wsize;
