1st tentative to unify phase description with polynomial fit. Tweak the...

1st tentative to unify phase description with polynomial fit. Tweak the TransformFitData function to trigger on correct transition point (tricky)

AnaFringes/hisold.pic

-exec fitphi.pic cor_CasA5mai18_0
+exec tmp2H.pic

AnaFringes/history.pic

-exec fitphi.pic cor_CasA6mai18
+exec tmp2H.pic

AnaFringes/postFringeV3.cc

@@ -67,6 +67,8 @@ struct PARAM {
 // Utilities
 //------------------------
 
+bool abscomp(r_8 a, r_8 b) {return fabs(a)<fabs(b);}
+
 void split(const string& str, const string& delimiters , vector<string>& tokens) {
     // Skip delimiters at beginning.
     string::size_type lastPos = str.find_first_not_of(delimiters, 0);
@@ -202,11 +204,15 @@ void TransformFitData(GeneralFitData& data, std::string tag=""){
     throw SzMismatchError("TransformFitData: FATAL: x empty");
   
   //adjacent difference to obtain the transition -Pi -> +Pi or vice versa
+
   vector<r_8>dy(y.size());
   std::adjacent_difference(y.begin(),y.end(),dy.begin());
   
-  r_8 minval = *min_element(dy.begin(),dy.end());
-  r_8 maxval = *max_element(dy.begin(),dy.end());
+  
+
+  //  r_8 minval = *min_element(dy.begin(),dy.end(),abscomp);
+  r_8 minval = median(dy);
+  r_8 maxval = *max_element(dy.begin(),dy.end(),abscomp);
   maxval = std::max(fabs(minval),fabs(maxval));
   minval = std::min(fabs(minval),fabs(maxval));
 
@@ -214,27 +220,23 @@ void TransformFitData(GeneralFitData& data, std::string tag=""){
   cout << "min/max dy(1): " << minval << " , " << maxval << endl;
 
   //Clean pathological points
-  
+
   vector<bool>nonvalidPts(x.size(),false);
   for(int i=1;i<dy.size()-1;i++){
-    if(fabs(dy[i])>minval*1.01 && fabs(dy[i])<0.99*maxval){
-      cout << "Patho: i,x,y,dy[i-1],dy[i],dy[i+1]: " 
-	   << i << ", " 
-	   << x[i] << ", "  
-	   << y[i] << ", "  
-	   << dy[i-1] << ", "  
-	   << dy[i] << ", "  
-	   << dy[i+1] << "\n";  
+    if(fabs(dy[i])>minval*1.5 && fabs(dy[i])<0.99*maxval){
+//       cout << "Patho: i,x,y,dy[i-1],dy[i],dy[i+1]: " 
+// 	   << i << ", " 
+// 	   << x[i] << ", "  
+// 	   << y[i] << ", "  
+// 	   << dy[i-1] << ", "  
+// 	   << dy[i] << ", "  
+// 	   << dy[i+1] << "\n";  
       nonvalidPts[i]=true;
       i++;
     }
   }
 
   //clean x,y,ey vectors: proceed backward to preserve iterator
-//   cout << "Dump BEFORE erase: ndata= " << x.size() <<  endl;
-//   for(int i=0; i<x.size(); i++){
-//     cout << "i,x,y: " << i << " " << x[i] << " " << y[i] << endl;
-//   }
   for(int i=nonvalidPts.size()-1;i>=0;i--){
     if(nonvalidPts[i]){
       cout << "erasing pt["<<i<<"]"<<endl;
@@ -243,10 +245,6 @@ void TransformFitData(GeneralFitData& data, std::string tag=""){
       ey.erase(ey.begin()+i);
     }
   }
-//   cout << "Dump AFTER erase: ndata= " << x.size() <<  endl;
-//   for(int i=0; i<x.size(); i++){
-//     cout << "i,x,y: " << i << " " << x[i] << " " << y[i] << endl;
-//   }
 
 
   {
@@ -705,53 +703,53 @@ void estimePhase(){
 
     GeneralFitData dataRes(1,freq.Size());
 
-    if(freqToTest.empty() || npos < 10 || nneg < 10){
+//     if(freqToTest.empty() || npos < 10 || nneg < 10){
 
-      cout << "Linear model choosen" << endl;
+//       cout << "Linear model choosen" << endl;
 
       
-      LinearFunc fLin;      
-      r_8 Par[3];
+//       LinearFunc fLin;      
+//       r_8 Par[3];
       
-      //reference point x:freq, y:phi
-      Par[0] = median(x);
-      Par[1] = median(y);
+//       //reference point x:freq, y:phi
+//       Par[0] = median(x);
+//       Par[1] = median(y);
 
-      //slope dy/dx
-      Par[2] = medslope;
+//       //slope dy/dx
+//       Par[2] = medslope;
       
-      cout << " ------------------------- " << endl;
-      cout << "The model selected is: linear with\n"
-	   << setprecision(20) 
-	   << " xref = " << Par[0] << " "
-	   << " yref = " << Par[1] << " "
-	   << " slope= " << Par[2] 
-	   << endl;
-      cout << " ------------------------- " << endl;
-
-      ofs << setprecision(20)
-	  << vPhaseName[ip] 
-	  << " Linear " 
-	  << Par[0] << " "
-	  << Par[1] << " "
-	  << Par[2] 
-	  << endl;
+//       cout << " ------------------------- " << endl;
+//       cout << "The model selected is: linear with\n"
+// 	   << setprecision(20) 
+// 	   << " xref = " << Par[0] << " "
+// 	   << " yref = " << Par[1] << " "
+// 	   << " slope= " << Par[2] 
+// 	   << endl;
+//       cout << " ------------------------- " << endl;
+
+//       ofs << setprecision(20)
+// 	  << vPhaseName[ip] 
+// 	  << " Linear " 
+// 	  << Par[0] << " "
+// 	  << Par[1] << " "
+// 	  << Par[2] 
+// 	  << endl;
 
 	
       
-      //compute the residuals
+//       //compute the residuals
 
-      for(int i=0; i<freq.Size(); i++){
-	r_8 xp[1]; xp[0]=freq(i); 
-	r_8 residual = phase(i) - fLin.Value(xp,Par);
-	dataRes.AddData1(xp[0],residual);
-      }
+//       for(int i=0; i<freq.Size(); i++){
+// 	r_8 xp[1]; xp[0]=freq(i); 
+// 	r_8 residual = phase(i) - fLin.Value(xp,Par);
+// 	dataRes.AddData1(xp[0],residual);
+//       }
 
       
-    } else {
+//     } else {
 
 
-      cout << "Sawtooth model choosen" << endl;
+//       cout << "Sawtooth model choosen" << endl;
       //On transforme en modele continu
       TransformFitData(dataFitReduced,vPhaseName[ip]);
       
@@ -772,6 +770,7 @@ void estimePhase(){
       }
       cout << "PolFit Success: dump: \n";
       pol.Print(cout);
+      cout << endl;
 
       ofs << setprecision(20)
  	  <<  vPhaseName[ip] 
@@ -793,7 +792,7 @@ void estimePhase(){
 
 
 
-    }//good model
+//     }//good model
     
     //save residuals      
       string nameTag = "Res"+vPhaseName[ip];
@@ -1136,6 +1135,7 @@ void rotateVisi() {
 //----------------------------------------------------------
 
 int main(int narg, char* arg[]) {
+
   SophyaInit();
 
   //message used in Exceptions