Index: trunk/ArchTOIPipe/ProcWSophya/simoffset.cc
===================================================================
--- trunk/ArchTOIPipe/ProcWSophya/simoffset.cc	(revision 2004)
+++ trunk/ArchTOIPipe/ProcWSophya/simoffset.cc	(revision 2008)
@@ -20,4 +20,5 @@
   mWSz = (mwsz > 8) ? mwsz : 8;
   nPtFit = (nptfit > degpol+2) ? nptfit : degpol+2;
+  if (nPtFit%2 == 0) nPtFit++;
   degPol = (degpol > 1)?degpol:1;
   totnscount = 0;
@@ -95,20 +96,23 @@
     int wsize = mWSz;
 
-    Vector vin(wsize);
+    int hisblk = nPtFit/2+1;
+    Vector vinhist(wsize*hisblk);
+    TVector<uint_8> vfghist(wsize*hisblk);
     Vector voff(wsize);
     Vector vout(wsize);
-    TVector<uint_8> vfg(wsize);
+    Vector vinc(wsize);
+    TVector<uint_8> vfgc(wsize);
+
     
     // Pour le fit 
-    Vector errCoef(3);
-    
-    Vector X(nPtFit);
-    Vector X0(nPtFit);
-    Vector Y(nPtFit);
-    Vector YErr(nPtFit);
+    Vector errCoef(degPol+1);
+    
+    Vector X(nPtFit+1);
+    Vector X0(nPtFit+1);
+    Vector Y(nPtFit+1);
+    Vector YErr(nPtFit+1);
 
     // Variables diverses 
-    int k,i,j,klast;
-    int nks = 0;
+    int k,kb,j,klast;
     klast = snb-1;
     totnbblock = 0;
@@ -124,10 +128,17 @@
     double sig = 0.;
     double meanx = 0.;
-    
+    double sn_last = 0.;
+    double y_last = 0.;
+
     // Boucle sur les sampleNum
-    // 1er partie, on traite par paquets de wsize
-
-    for(k=snb;k<=sne-wsize+1;k+=wsize) {
+    // 1ere partie, on traite par paquets de wsize
+    int nblk = (sne-snb+1)/wsize; 
+    for(kb=0; kb<nblk; kb++) {
+      k = kb*wsize+snb;
+      //    for(k=snb;k<=sne-wsize+1;k+=wsize) {
       // Lecture d'un bloc de donnees
+      Vector vin( vinhist(Range((kb%hisblk)*wsize, 0, wsize) ) );
+      TVector<uint_8> vfg( vfghist(Range((kb%hisblk)*wsize, 0, wsize) ) );
+
       getData(0, k, wsize, vin.Data(), vfg.Data());
 
@@ -144,5 +155,4 @@
       }
       
-      sn0 = (double)(k+wsize/2);
       
       if (!fginiXYdone) {
@@ -152,5 +162,5 @@
 	  sig = sig/nok-mean*mean;	
 	}
-	X = RegularSequence(k+wsize*0.5, (double)wsize);
+	X = RegularSequence((kb+0.5)*wsize, (double)wsize);
 	Y = mean;
 	YErr = (nok > 0.5) ? sqrt(mean) : 1.;
@@ -158,8 +168,9 @@
       }
 
-      if (nok > 3.) {
+      if (nok > 3.5) {
 	mean /= nok;
 	meanx /= nok;
 	sig = sig/nok-mean*mean;
+	if (sig < 1.e-6*mean) sig = 1.e-6*mean;
 	int kk = nbblkok%nPtFit;
 	nbblkok++;
@@ -169,65 +180,32 @@
       }
       
-      X0 = X;
-      X0 -= sn0;
-      if (nbblkok > poly.Degre()+1) 
+      if ((nok>3.5) && (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);
+	YErr(nPtFit) = smax/10.;
+	X(nPtFit) = sn_last;
+	sn0 = (double)(k-nPtFit*wsize*0.5);
+	X0 = X;
+	X0 -= sn0;
 	poly.Fit(X0,Y,YErr,degPol,errCoef);
+      }
       else {
-	poly[0] = mean;
-	for(int jj=1; jj<=poly.Degre(); jj++) poly[jj] = 0.;
-      }
-      /*
-      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;
-      }
-      */
-      
-      // Calcul des valeurs d'offset en sortie
-      for(j=0; j<wsize; j++) 
-	voff(j) = poly(k+j-sn0);
-      
-      if (fgoffset) putData(0, k, wsize, voff.Data());
-      if (fgincopie) putData(2, k, wsize, vin.Data(), vfg.Data());
-      if (fgout) { 
-	vin -= voff;
-	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());
-      }
-      if (fgsigy) {
-	vout = sig;
-	putData(8, k, wsize, vout.Data());
-      }
-
-      if (fgmeanx) {
-	vout = meanx;
-	putData(9, k, wsize, vout.Data());
-      }
-
+	if (nbblkok < 2) {
+	  sn0 = k+1;
+	  poly[0] = mean;
+	  for(int jj=1; jj<=poly.Degre(); 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));
+	  sn0 =  0.5*(X(nbblkok-1)+X(0));
+	  for(int jj=2; jj<=poly.Degre(); jj++) poly[jj] = 0.;
+	} 
+	sn_last = sn0;
+	y_last = poly(sn_last-sn0);
+      }
       if (ntpoly) {    // Remplissage du XNTuple de controle 
-	char * nomsnt[] = {"sncur", "sn0", "meanx", "meany", "sigy", "a0", "a1", "a2", "ycur"};
 	float xnt[10];
 	xnt[0] = k;
@@ -243,4 +221,84 @@
       }
 
+      /*
+      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;
+      TVector<uint_8> vfgc;
+      int kbs = kb-nPtFit/2;
+      k = kbs*wsize+snb;
+      if (kb == nblk-1) {
+	int wszt = wsize*hisblk;
+	voff.ReSize(wszt);
+	vout.ReSize(wszt);	
+	vinc.ReSize(wszt); 
+	vfgc.ReSize(wszt);
+	int jb = 0;
+	for(int kbsc=kbs; kbsc<nblk; kbsc++) {
+	  vinc(Range(jb*wsize, 0, wsize)) = 
+	       vinhist(Range((kbsc%hisblk)*wsize, 0, wsize) ) ;
+	  vfgc(Range(jb*wsize, 0, wsize)) = 
+	       vfghist(Range((kbsc%hisblk)*wsize, 0, wsize) ) ;
+	  jb++;     
+	}
+	wsize = wszt;
+      }
+      else {
+	vinc = vinhist(Range((kbs%hisblk)*wsize, 0, wsize) ) ;
+	vfgc = vfghist(Range((kbs%hisblk)*wsize, 0, wsize) ) ;
+      }
+
+      // Calcul des valeurs d'offset en sortie
+      for(j=0; j<wsize; j++) 
+	voff(j) = poly(k+j-sn0);
+      sn_last = k+wsize;
+      y_last = poly(sn_last-sn0);
+
+      if (fgoffset) putData(0, k, wsize, voff.Data());
+      if (fgout) { 
+	vinc -= voff;
+	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());
+      }
+      if (fgsigy) {
+	vout = sig;
+	putData(8, k, wsize, vout.Data());
+      }
+
+      if (fgmeanx) {
+	vout = meanx;
+	putData(9, k, wsize, vout.Data());
+      }
+
       klast+=wsize;
       totnscount+=wsize;
@@ -249,11 +307,12 @@
     } // Fin boucle sur les samples, par pas de wsize
 
-    // 2eme partie, on traite la fin du bloc d'echantillons si necessaire
+    // 3eme partie, on traite la fin du bloc d'echantillons si necessaire
     if (klast < sne) {
       wsize = sne-klast;
-      vin.ReSize(wsize);
+      Vector vin(wsize);
       voff.ReSize(wsize);
       vout.ReSize(wsize);
-      vfg.ReSize(wsize);
+      TVector<uint_8> vfg(wsize);
+      k = klast+1;
       getData(0, k, wsize, vin.Data(), vfg.Data());
       for(j=0; j<wsize; j++) 
