Index: trunk/ArchTOIPipe/ProcWSophya/simoffset.cc
===================================================================
--- trunk/ArchTOIPipe/ProcWSophya/simoffset.cc	(revision 2040)
+++ trunk/ArchTOIPipe/ProcWSophya/simoffset.cc	(revision 2045)
@@ -12,9 +12,11 @@
 #include "ctimer.h"
 #include "ntuple.h"
+#include "sopemtx.h"
 
 #include "flagtoidef.h"
 
+
 SimpleOffsetEstimator::SimpleOffsetEstimator(int mwsz, int nptfit, int degpol)
-  : poly((degpol > 1)?degpol:1)
+  : poly((degpol > 1)?degpol+1:2)
 {
   SetParams(mwsz, nptfit, degpol);
@@ -24,4 +26,5 @@
   bgalcut = false;
   ns_flgcut = ns_bgalcut = 0;
+  npb_fitpoly = 0;
   SavePolyNTuple();
 }
@@ -51,5 +54,5 @@
   os << "\n ------------------------------------------------------ \n" 
      << " SimpleOffsetEstimator::PrintStatus() - MeanWSize= " << mWSz << " NPtFit=" 
-     << nPtFit << " DegPoly=" << degPol << " poly.Degre()=" << poly.Degre() << endl;
+     << nPtFit << " DegPoly=" << degPol << endl;
   if (bgalcut) 
     os << " bGalCut = Yes , bGalMin = " << bmincut << " bGalMax= " << bmaxcut << endl;
@@ -57,5 +60,6 @@
   TOIProcessor::PrintStatus(os);
   os << " ProcessedSampleCount=" << ProcessedSampleCount() 
-     << " NS_FlgCut= " << ns_flgcut << " NS_BGalCut= " << ns_bgalcut << endl;
+     << " NS_FlgCut= " << ns_flgcut << " NS_BGalCut= " << ns_bgalcut 
+     << " NPb_FitPoly= " << npb_fitpoly << endl;
   os << " ------------------------------------------------------ " << endl;  
 }
@@ -118,4 +122,6 @@
   NTuple xntp(11, nomsnt);
 
+  char ans[20];
+
   try {
 
@@ -168,4 +174,7 @@
     int nblk = (sne-snb+1)/wsize; 
     for(kb=0; kb<nblk; kb++) {
+      //      cout << " SimpleOffsetEstimator::run - Loop " << kb << " NbBlkOK= "  
+      //	   << nbblkok << " <CR> to continue / q --> QUIT ... \n" ; 
+      //     gets(ans); if (ans[0] == 'q') break;
       k = kb*wsize+snb;
       //    for(k=snb;k<=sne-wsize+1;k+=wsize) {
@@ -236,5 +245,6 @@
 	fg_last_meansig = true;
       }
-      
+
+            
       if (((nok>3.5) || fg_last_meansig) && (nbblkok >= nPtFit) ) {
 	// On force le fit a garder une certaine continuite
@@ -248,31 +258,22 @@
 	X0 = X;
 	X0 -= sn0;
-	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);
-	}
+	FitPoly(X0, Y, YErr);
       }
       else {
 	if (nbblkok < 2) {
 	  sn0 = k+1;
-	  poly[0] = mean;
-	  for(int jj=1; jj<=poly.Degre(); jj++) poly[jj] = 0.;
+	  poly(0) = mean;
+	  for(int jj=1; jj<=degPol; jj++) poly(jj) = 0.;
 	}
 	else if (nbblkok < nPtFit) {
-	  poly[0] = 0.5*(Y(nbblkok-1)+Y(0));
-	  poly[1] = (Y(nbblkok-1) - Y(0))/(X(nbblkok-1)-X(0));
+	  poly(0) = 0.5*(Y(nbblkok-1)+Y(0));
+	  poly(1) = (Y(nbblkok-1) - Y(0))/(X(nbblkok-1)-X(0));
 	  sn0 =  0.5*(X(nbblkok-1)+X(0));
-	  for(int jj=2; jj<=poly.Degre(); jj++) poly[jj] = 0.;
+	  for(int jj=2; jj<=degPol; jj++) poly(jj) = 0.;
 	} 
 	sn_last = sn0;
 	y_last = poly(sn_last-sn0);
       }
+      
       if (ntpoly) {    // Remplissage du XNTuple de controle 
 	float xnt[12];
@@ -282,8 +283,8 @@
 	xnt[3] = mean;
 	xnt[4] = sig;
-	xnt[5] = poly[0];
-	xnt[6] = poly[1];
-	xnt[7] = poly[2];
-	xnt[8] = poly(k-sn0);
+	xnt[5] = poly(0);
+	xnt[6] = poly(1);
+	xnt[7] = poly(2);
+	xnt[8] = ApplyPoly(k-sn0);
 	xnt[9] = nok;
 	xnt[10] = nbblkok;
@@ -291,18 +292,9 @@
       }
 
-      /*
-      if (nbblkok < 8) {
-	cout << "------ DBG-X " << nbblkok << "," << nok << " degre=" << poly.Degre() << endl;
-	cout << "DBG-A X0=" << X0 << endl;
- 	cout << "DBG-A Y=" << Y << endl;
- 	cout << "DBG-A YErr=" << YErr << endl;
-	cout << "DBG-A poly= " << poly << endl;
-      }
-      */
-      
+     
       if (fgincopie) putData(2, k, wsize, vin.Data(), vfg.Data());
 
       if (kb < nPtFit/2) continue;
-      Vector vinc;
+      Vector vinc; 
       TVector<uint_8> vfgc;
       int kbs = kb-nPtFit/2;
@@ -332,8 +324,7 @@
       // Calcul des valeurs d'offset en sortie
       for(j=0; j<wsize; j++) 
-	voff(j) = poly(k+j-sn0);
+	voff(j) = ApplyPoly(k+j-sn0);
       sn_last = k+wsize;
-      y_last = poly(sn_last-sn0);
-
+      y_last = ApplyPoly(sn_last-sn0);
       if (fgoffset) putData(0, k, wsize, voff.Data());
       if (fgout) { 
@@ -341,24 +332,9 @@
 	putData(1, k, wsize, vinc.Data(), vfgc.Data());
       }
-      /*
-      if (fgmeany) {
-	vout = mean;
-	putData(4, k, wsize, vout.Data());
-      }
-      if (fgsigy) {
-	vout = sig;
-	putData(5, k, wsize, vout.Data());
-      }
-
-      if (fgmeanx) {
-	vout = meanx;
-	putData(6, k, wsize, vout.Data());
-      }
-      */
 
       klast+=wsize;
       totnscount+=wsize;
       totnbblock++;
-      
+      //      cout << " SimpleOffset::run " << totnbblock << " totnscount " << totnscount << endl;
     } // Fin boucle sur les samples, par pas de wsize
 
@@ -379,5 +355,5 @@
 
       for(j=0; j<wsize; j++) 
-	voff(j) = poly(k+j-sn0);
+	voff(j) = ApplyPoly(k+j-sn0);
       if (fgoffset) putData(0, k, wsize, voff.Data());
       if (fgincopie) putData(2, k, wsize, vin.Data(), vfg.Data());
@@ -425,2 +401,52 @@
   }
 }
+
+
+double SimpleOffsetEstimator::ApplyPoly(double x)
+{
+  if (degPol == 0) return(poly(0));
+  else if (degPol == 1)  return(poly(0)+poly(1)*x);
+  else if (degPol == 2)  return(poly(0)+poly(1)*x+poly(2)*x*x);
+  else if (degPol == 3)  return(poly(0)+poly(1)*x+poly(2)*x*x+poly(3)*x*x*x);
+  else {
+    double ret = poly(0);
+    double xx = x;
+    for(int k=1; k<degPol+1; k++)  { ret += xx*poly(k);  xx *= x; }
+    return(ret);
+  }
+}
+
+// Fonction pour faire LinFit d'un polynome
+r_8 rzpoly_f(uint_4 k, r_8 x )
+{
+  if (k == 0) return(1.);
+  else if (k == 1) return(x);
+  else if (k == 2) return(x*x);
+  else if (k == 3) return(x*x*x);
+  else {
+    r_8 ret = x;
+    for(int i=1; i<k; i++) ret *= x;
+    return(ret);
+  }
+}
+
+void SimpleOffsetEstimator::FitPoly(Vector& X0, Vector& Y, Vector& YErr)
+{
+  Vector cpol(degPol+1);
+  Vector errcoef(degPol+1);
+  LinFitter<r_8> lfit;
+  try {
+    lfit.LinFit(X0, Y, YErr, degPol+1, rzpoly_f, cpol, errcoef);
+    poly = cpol;
+  }
+  catch (PException & exc) {
+    if (npb_fitpoly < 50) 
+      cout << " -- SimpleOffsetEstimator::FitPoly/ Catched Exception " 
+	   << (string)typeid(exc).name() << "\n .... Msg= " << exc.Msg() << endl;
+    if (npb_fitpoly < 10) 
+      cout << " X0= " << X0 << " Y=" << Y << " YErr=" << YErr << endl;
+    npb_fitpoly++;
+  }
+  return;
+}
+
Index: trunk/ArchTOIPipe/ProcWSophya/simoffset.h
===================================================================
--- trunk/ArchTOIPipe/ProcWSophya/simoffset.h	(revision 2040)
+++ trunk/ArchTOIPipe/ProcWSophya/simoffset.h	(revision 2045)
@@ -11,5 +11,4 @@
 #include "toiprocessor.h"
 #include "tvector.h"
-#include "poly.h"
 
 // ---- Calcul de ligne de base 
@@ -42,4 +41,7 @@
   inline  void  SavePolyNTuple(bool fg=false, string name="") { ntpoly = fg; ntpolyname=name; }
 protected:
+  virtual double ApplyPoly(double x);
+  virtual void   FitPoly(Vector& X0, Vector& Y, Vector& YErr);
+
   int_8 totnscount;   // Nombre total d'echantillon processe
   int_8 totnbblock;   // Nombre total de blocs 
@@ -47,5 +49,7 @@
   int nPtFit;
   int degPol;
-  Poly poly;
+  //  Poly poly;   Pb SEGV lors de Poly.Fit
+  Vector poly;
+  int_4 npb_fitpoly;
   bool ntpoly;
   string ntpolyname;
