#include "machdefs.h"

//---- System et stdc++ include files
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <ctype.h>
#include <string.h>
#include <iostream>
#include <fstream>
#include <complex>

#include <typeinfo>
#include <string>
#include <vector>
#include <map>
#include <functional>
#include <list>
#include <algorithm> //JEC
#include <numeric>   //JEC

//---- Sophya include files
#include "sopnamsp.h"
#include "basetools.h"
#include "systools.h"
#include "sutils.h"
#include "ntools.h"
#include "array.h"
#include "histats.h"

//---- Spiapp include files
#include "nobjmgr.h"
#include "servnobjm.h"

//---- Include files from additionnal modules

//---- function to compare bits on double 
int_8 BitCmp64(double v,int_8 flg) 
{return ((int_8)((v<0.) ? v-0.1 : v+0.1))&flg;} 

//---- function for Adding and displaying Objects  
void Keep_Object(AnyDataObj & obj, string const & nom) 
{ 
  string name = nom; 
  NamedObjMgr om; 
  if (om.GetObj(name)) 
    cerr << "KeepObj()/Warning Already kept object " << endl; 
  else om.AddObj(obj, name); 
} 

void Display_Object(AnyDataObj & obj, string const & opt, string const & nom) 
 { 
  string name = nom; 
  NamedObjMgr om; 
  if (!om.GetObj(name)) 
    om.AddObj(obj, name); 
  om.DisplayObj(name, opt); 
} 

//---- function for getting and setting ObjectManager variables  
void Set_ObjMgrVar(MuTyV v, string const & nom) 
{ 
  NamedObjMgr om; 
  om.SetVar(nom, (string)v); 
} 

MuTyV Get_ObjMgrVar(const char * nom) 
{ 
  string name = nom; NamedObjMgr om; 
  MuTyV v = om.GetVar(name); 
  return v; 
} 

//---- Macro for Objects and variables saving
#define KeepObj(obj) Keep_Object(obj, #obj)
#define GetOMVar(var) Get_ObjMgrVar( #var )
#define SetOMVar(var) Set_ObjMgrVar(var, #var )

//---- Macro Displaying objects and command execution
#define DispObj(obj, att) Display_Object(obj, att, #obj); 

#define ExecCmd(cmd) srvo.ExecuteCommand(cmd); 



//-------------------------------------------------//
//----------------- User Functions ----------------//
//-------------------------------------------------//

extern "C" {
  int fcalib( vector<string>& args );
   int func( vector<string>& args );
}



int fcalib( vector<string>& args )
{
// Some definitions to help using spiapp;
NamedObjMgr omg;
Services2NObjMgr& srvo = *omg.GetServiceObj();

//-------------- Object List --------------
//Number of objects = 3
string ___nomobj;

___nomobj = "datav0";
NTuple * ___datav0 = dynamic_cast< NTuple  * >(omg.GetObj(___nomobj));
if(___datav0==NULL) throw NullPtrError("CxxExecutor::PutObject: Non existing object datav0... please update file");
NTuple & datav0 = (*___datav0);

___nomobj = "datav1";
NTuple * ___datav1 = dynamic_cast< NTuple  * >(omg.GetObj(___nomobj));
if(___datav1==NULL) throw NullPtrError("CxxExecutor::PutObject: Non existing object datav1... please update file");
NTuple & datav1 = (*___datav1);

___nomobj = "ntcal";
NTuple * ___ntcal = dynamic_cast< NTuple  * >(omg.GetObj(___nomobj));
if(___ntcal==NULL) throw NullPtrError("CxxExecutor::PutObject: Non existing object ntcal... please update file");
NTuple & ntcal = (*___ntcal);




//--------------------------------------------//
//----------------- User Code ----------------//
//--------------------------------------------//

 r_4 tRef     = atof(args[0].data());    // start of SCA calibration
 r_4 cycle = atof(args[1].data()); // numero du cycle
 int_4 nentry = datav0.NEntry(); //suppose datav1.NEntry= datav0.NEntry
 r_4 minInt0 = atof(args[2].data()); //min intensity Ch 0
 r_4 maxInt0 = atof(args[3].data()); //max intensity Ch 0
 r_4 minInt1 = atof(args[4].data()); //min intensity Ch 1
 r_4 maxInt1 = atof(args[5].data()); //max intensity Ch 1

 // cout << "min/max 0/1: " << minInt0 << " " << maxInt0 
 //     << " " <<  minInt1 << " " << maxInt1 << endl;
 r_4 deltathres0 = 0.80 * (maxInt0-minInt0); // seuil sur l'increment de l'intensite du a la DAB Ch0
 r_4 deltathres1 = 0.80 * (maxInt1-minInt1); // seuil sur l'increment de l'intensite du a la DAB Ch1
 // select time E [twin0,twin1] U [twin2,twin3] for Tsys
 // time unit = second 
 r_4 twin0 = -3.;
 r_4 twin1 = -1.;
 r_4 twin2 =  6.;
 r_4 twin3 =  8;

 std::vector<r_4> lowerIntv0;
 std::vector<r_4> lowerIntv1;
 std::vector<r_4> upperIntv0;
 std::vector<r_4> upperIntv1;

 // select values to compute background Tsys
 for (int_4 ie=0;ie<nentry;ie++){
   r_4 difftime0 = datav0.GetVal(ie,"time")-tRef;
   r_4 x0 = datav0.GetVal(ie,"x0");

   r_4 difftime1 = datav1.GetVal(ie,"time")-tRef;
   r_4 x1 = datav1.GetVal(ie,"x1");

   if ( (difftime0 >= twin0 && difftime0 <= twin1) ||
	(difftime0 >= twin2 && difftime0 <= twin3)
	) lowerIntv0.push_back(x0);
   if ( (difftime1 >= twin0 && difftime1 <= twin1) ||
	(difftime1 >= twin2 && difftime1 <= twin3)
	) lowerIntv1.push_back(x1);
 } //end for

 //compute the Tsys for both channels
 std::nth_element(lowerIntv0.begin(),
		  lowerIntv0.begin()+lowerIntv0.size()/2,
		  lowerIntv0.end());
 r_4 tsys0 = *(lowerIntv0.begin()+lowerIntv0.size()/2);

 std::nth_element(lowerIntv1.begin(),
		  lowerIntv1.begin()+lowerIntv1.size()/2,
		  lowerIntv1.end());
 r_4 tsys1 = *(lowerIntv1.begin()+lowerIntv1.size()/2);
 
 //set the threshold for DAB detection
 r_4 thres0 =  tsys0 + deltathres0;
 r_4 thres1 =  tsys1 + deltathres1;
 
 cout << "thres 0: " << thres0 << " 1: " << thres1 <<endl;

 for (int_4 ie=0;ie<nentry;ie++){
   r_4 x0 = datav0.GetVal(ie,"x0");
   r_4 x1 = datav1.GetVal(ie,"x1");
   
   if(x0 >= thres0) upperIntv0.push_back(x0);
   if(x1 >= thres1) upperIntv1.push_back(x1);

 }//end for

 cout << "size upper 0: " << upperIntv0.size() << " 1: " << upperIntv1.size() << endl;

 std::vector<r_4> upperIntAdjDiffv0(upperIntv0.size());
 std::adjacent_difference(upperIntv0.begin(),
			  upperIntv0.end(),
			  upperIntAdjDiffv0.begin());

 //collect data for at most 2 DAB states
 std::vector<r_4> upIntStatev0[2];
 int state=-1;
 for (int i=1;i<upperIntv0.size();i++) {//skip first value
   if (fabs(upperIntAdjDiffv0[i])<upperIntv0[i]*0.010) {
     if(state == -1) state=0;
     if(state == -2) state=1;
     upIntStatev0[state].push_back(upperIntv0[i]);
   } else {
     if (state == 0) state=-2;
   }
 }

 r_4 mean = 0;
 r_4 nval = 0;
 for (int i=0;i<2;i++) {
   if (!upIntStatev0[i].empty()) {
     std::nth_element(upIntStatev0[i].begin(),
		  upIntStatev0[i].begin()+upIntStatev0[i].size()/2,
		  upIntStatev0[i].end());
     mean += *(upIntStatev0[i].begin()+upIntStatev0[i].size()/2);
     nval++;
   } 
 }
 mean /= nval;

 r_4 deltaInt0 = mean - tsys0;


 ////
 std::vector<r_4> upperIntAdjDiffv1(upperIntv1.size());
 std::adjacent_difference(upperIntv1.begin(),
			  upperIntv1.end(),
			  upperIntAdjDiffv1.begin());

 //collect data for at most 2 DAB states
 std::vector<r_4> upIntStatev1[2];
 state=-1;
 for (int i=1;i<upperIntv1.size();i++) {//skip first value
   if (fabs(upperIntAdjDiffv1[i])<upperIntv1[i]*0.010) {
     if(state == -1) state=0;
     if(state == -2) state=1;
     upIntStatev1[state].push_back(upperIntv1[i]);
   } else {
     if (state == 0) state=-2;
   }
 }

 mean = 0;
 nval = 0;
 for (int i=0;i<2;i++) {
   if (!upIntStatev1[i].empty()) {
     std::nth_element(upIntStatev1[i].begin(),
		  upIntStatev1[i].begin()+upIntStatev1[i].size()/2,
		  upIntStatev1[i].end());
     mean += *(upIntStatev1[i].begin()+upIntStatev1[i].size()/2);
     nval++;
   } 
 }
 mean /= nval;

 r_4 deltaInt1 = mean - tsys1;


 //save & return
 r_4 out[5];
 out[0] = cycle;
 out[1] = tsys0; out[2] = tsys1;
 out[3] = deltaInt0; out[4] = deltaInt1;

//  for (int i=0;i<5;i++)
//    cout << "out["<<i<<"]="<<out[i]<<endl;

 ntcal.Fill(out);

 return 0;
}



int func( vector<string>& args )
{
// Some definitions to help using spiapp;
NamedObjMgr omg;
Services2NObjMgr& srvo = *omg.GetServiceObj();

//-------------- Object List --------------
//Number of objects = 1
string ___nomobj;

___nomobj = "inarr";
TVector< r_4 > * ___inarr = dynamic_cast< TVector< r_4 >  * >(omg.GetObj(___nomobj));
if(___inarr==NULL) throw NullPtrError("CxxExecutor::PutObject: Non existing object inarr... please update file");
TVector< r_4 > & inarr = (*___inarr);



//--------------------------------------------//
//----------------- User Code ----------------//
//--------------------------------------------//
//#include "/home/bao/tmp/PIATmp_iAIt8R/cxx_spiapp.h"

//Input: inarr vector: X=> freq
//Output: outarr vector: X=> freq
//Center a window @ 1410 MHz; the DAQ takes 250MHz starting at 1250MHz
//RT filter width for calibration : 6.25MHz  
 r_4 f0 = atof(args[0].data());
 r_4 RTFilterWidth = atof(args[1].data());
 sa_size_t c0 = (sa_size_t)inarr.NElts()*(f0-1250.)/250.;
 sa_size_t dc = (sa_size_t)inarr.NElts()*RTFilterWidth/250.;
 sa_size_t fMin = c0-dc/2; //lower bound of the freq. filter
 sa_size_t fMax = c0+dc/2; //upper bound of the freq. filter

 TVector<r_4> outarr(inarr.SubVector(Range(fMin,fMax)),false);//freq. bin 0 is a DC level
 
 KeepObj(outarr);


 return 0;
}
