Changeset 2045 in Sophya for trunk/ArchTOIPipe/ProcWSophya/simoffset.cc

Timestamp:

Jun 7, 2002, 4:21:22 PM (23 years ago)

Author:

ansari

Message:

Remplacement de l'objet poly (pb plantage ds lors du Fit) par un vecteur - Reza 7/6/2002

File:

• trunk/ArchTOIPipe/ProcWSophya/simoffset.cc (modified) (14 diffs)

trunk/ArchTOIPipe/ProcWSophya/simoffset.cc

@@ -12 +12 @@
 #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 +26 @@
   bgalcut = false;
   ns_flgcut = ns_bgalcut = 0;
+  npb_fitpoly = 0;
   SavePolyNTuple();
 }
@@ -51 +54 @@
   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 +60 @@
   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 +122 @@
   NTuple xntp(11, nomsnt);
 
+  char ans[20];
+
   try {
 
@@ -168 +174 @@
     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 +245 @@
         fg_last_meansig = true;
       }
-      
+
+            
       if (((nok>3.5) || fg_last_meansig) && (nbblkok >= nPtFit) ) {
         // On force le fit a garder une certaine continuite
@@ -248 +258 @@
         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 +283 @@
         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 +292 @@
       }
 
-      /*
-      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 +324 @@
       // 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 +332 @@
         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 +355 @@
 
       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 +401 @@
   }
 }
+
+
+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;
+}
+