/*  ------------------------ Projet BAORadio -------------------- 
  Programme de verification comportement en loi de puissance 
    R. Ansari , Sep 2011

  Usage: ckpowl In3DMap  [InRenormFactor] 
---------------------------------------------------------------  */

#include "machdefs.h"
#include "sopnamsp.h"
#include <iostream>
#include <string>
#include <math.h>

#include <typeinfo>

#include "specpk.h"
#include "histats.h"

#include "qhist.h" 
#include "lobe.h" 
#include "cubedef.h" 
#include "fgndsub.h" 

#include "histinit.h"
#include "fftwserver.h"
#include "randr48.h"      

#include "ctimer.h"

typedef ThSDR48RandGen RandomGenerator ;

//-------------------------------------------------------------------------
//      ------------------ MAIN PROGRAM ------------------------------
//-------------------------------------------------------------------------
int main(int narg, const char* arg[])
{
  if (narg<3) {
    cout << " Usage: calcpk In3DMap_PPFName OutName [InRenormFactor] [PixNoiseLevel] " << endl;
    return 1;
  }
  Timer tm("calcpk");
  int rc = 0;
  try {
    string inppfname = arg[1];
    string outname = arg[2];
    double renfact=1.;
    bool fgrenfact=false;
    if (narg>3) {
      renfact=atof(arg[3]);
      fgrenfact=true;
    }
    double pixsignoise = 0.;
    bool fgaddnoise=false;
    if (narg>4) {
      pixsignoise=atof(arg[4]);
      if (pixsignoise>1.e-9)  fgaddnoise=true;
    }


    cout << "ckpowl[1] : reading 3D map from file " << inppfname << endl;
    TArray<r_4> inmap;
    {
      PInPersist pin(inppfname);
      pin >> inmap;
    }
    if (fgrenfact) { 
      cout << "  Applying RenormFactor inmap = inmap*rfact, rfact=" << renfact << endl;
      inmap *= renfact;
    }
    double mean, sigma;
    MeanSigma(inmap, mean, sigma);
    cout << " InMap sizes " << inmap.InfoString() << endl;
    inmap.Show(); 
    cout << "  Mean=" << mean << " Sigma=" << sigma << endl;
    if (fgaddnoise) { 
      BeamEffect::AddNoise(inmap, pixsignoise);
    }

    tm.Split(" After read ");
    cout << "ckpowl[2] : Checking power law behaviour " << endl;
    PowerLawChecker plck(inmap);
    plck.CheckXYMean();
    tm.Split(" Done PowerLawChecker::CheckXYMean() ");    
  }  // End of try bloc 
  catch (PThrowable & exc) {  // catching SOPHYA exceptions
    cerr << " ckpowl.cc: Catched Exception (PThrowable)" << (string)typeid(exc).name() 
         << "\n...exc.Msg= " << exc.Msg() << endl;
    rc = 99;
  }
  catch (std::exception & e) {  // catching standard C++ exceptions
    cerr << " ckpowl.cc: Catched std::exception "  << " - what()= " << e.what() << endl;
    rc = 98;
  }
  catch (...) {  // catching other exceptions
    cerr << " ckpowlx.cc: some other exception (...) was caught ! " << endl;
    rc = 97;
  }
  cout << " ==== End of ckpowl.cc program  Rc= " << rc << endl;
  return rc;    
}