#ifndef _define_h_
#define _define_h_
#include<iostream.h>
#include<fstream.h>
#include<math.h>
#include <vector>
#include <algorithm>
#include<string>
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>
#include "../../randomc.h"
#ifndef _NO_NAMESPACE
using namespace std;
#define STD std
#else
#define STD
#endif
#endif
void
Hist (vector < double >&y, vector < double >&out, int bins)
{
for (int i = 0; i < bins; i++)
out.push_back (0.);
sort (y.begin (), y.end ());
double width = y[y.size () - 1] / bins;
for (int i = 0; i < bins; i++)
for (unsigned int j = 0; j < y.size (); j++)
{
if (y[j] >= i * width && y[j] < (i + 1) * width)
out[i] += 1;
}
}
void
FirstOrderStat (vector < double >&y, double &mean, double &std)
{
double sum = 0., sumsq = 0.;
sort (y.begin (), y.end ());
for (unsigned int i = 0; i < y.size (); i++)
sum += y[i];
mean = sum / y.size ();
for (unsigned int i = 0; i < y.size (); i++)
sumsq += pow ((y[i] - mean), 2.0);
std = sqrt (sumsq / (y.size () - 1));
}
void
Frate (double *t, vector < double >&y, double &r, double resolution,int filter)
{
double sum = 0;
for (unsigned int i = 0; i < y.size (); i++)
{
if (filter == 1)
{
if (*t - y[i] > -resolution/2 && *t - y[i] <= resolution/2)
sum += 1. / resolution;
}
if (filter == 2)
{
if (*t - y[i] > 0)
sum += pow((1.0/resolution),2.0) * (*t - y[i]) * exp (-(1.0/resolution) * (*t - y[i]));
}
if (filter == 3)
{
sum +=
(1. / (sqrt (2 * 3.14) * resolution)) *
exp (-pow ((*t - y[i]) / (resolution), 2.));
}
}
r = sum;
}
void Correlate(vector<double>&x,vector<double>&y,vector <double>&C)
{
double mx,my,sigx,sigy;
double temp;
FirstOrderStat(x,mx,sigx);
FirstOrderStat(y,my,sigy);
for (unsigned int i=0;i<x.size();i++)
{
temp=0;
for (unsigned int j=0;j<x.size();j++)
{
if (i+j<y.size())
temp+=(x[j]-mx)*(y[j+i]-my)/(sigx*sigy);
}
C.push_back(temp/(x.size()-1));
}
}