Commit e37109b7 authored by Reza  ANSARI's avatar Reza ANSARI
Browse files

Modifs cosmetiques ds gacfit.h et ajout embryon fichier gcxfit.h pour fit...

Modifs cosmetiques ds gacfit.h et ajout embryon fichier gcxfit.h pour fit cross-cor, mais  , Reza 20/12/2018
parent 269f56a6
/* PAON4 analysis software
Class representing an autocorrelation beam (adapted from JSkyMap code)
Class for fitting auto-correlations (antenna geometry)
R. Ansari, December 2018 */
#ifndef GACFIT_H_SEEN
......@@ -11,7 +11,7 @@
#include "acbeam.h"
static int debug_galfit=0;
static int debug_gacfit=0;
//-----------------------------------------------------------------
class MyACSignal {
......@@ -51,7 +51,7 @@ public:
virtual int getExpectedSignal(size_t j, vector<double> & sig, double D_dish, double theta_b, double phi_b, double A, double B)
{
//DBG if (debug_galfit > 0) cout<<"*DBG*A*getExpectedSignal() j="<<j<<" size="<<v_time_data[j].size()<<endl;
//DBG if (debug_gacfit > 0) cout<<"*DBG*A*getExpectedSignal() j="<<j<<" size="<<v_time_data[j].size()<<endl;
sig.resize(v_time_data[j].size());
vector<double> & vtime = v_time_data[j];
vector<double> & vdata = vv_data[j][nac_];
......@@ -97,7 +97,7 @@ public:
{
vector<double> sig;
//DBG cout << "*DBG*getExpectedSignal() : calling getExpectedSignal(j="<<j<<" ,sig,...)"<<endl;
//DBG debug_galfit=1;
//DBG debug_gacfit=1;
getExpectedSignal(j, sig, D_dish, theta_b, phi_b, A, B);
TVector<double> rvs(sig.size());
for(size_t k=0; k<sig.size(); k++) rvs(k)=sig[k];
......
/* PAON4 analysis software
Class for fitting auto-correlations (antenna geometry)
R. Ansari, December 2018 */
#ifndef GCXFIT_H_SEEN
#define GCXFIT_H_SEEN
#include "generalfit.h"
#include "simplex.h"
#include "sunitpcst.h"
#include "cxbeam.h"
class MyCxSignal {
public:
MyCxSignal(vector< vector<double> > & v_time_data_, vector< vector< vector< complex<double> > > > & vv_cxdata_,
vector< vector< vector<double> > > & vv_cxerr_, vector<double> & v_freqs_,
vector< SLinInterp1D > & v_interp_elev_, vector< SLinInterp1D > & v_interp_sazim_,
size_t ncx=0, bool fggauss_beam=true)
: v_time_data(v_time_data_), vv_data(vv_data_), vv_err(vv_err_), v_freqs(v_freqs_),
v_interp_elev(v_interp_elev_), v_interp_sazim(v_interp_sazim_),
nac_(nac),
{
cout << " MyACSignal() Antenna=nac+1= "<<nac_+1<<" NPar="<< 3+2*v_time_data_.size() << endl;
double clight = PhysQty::c().SIValue();
v_lambdas.resize(v_freqs.size());
for(size_t k=0; k<v_freqs.size(); k++)
v_lambdas[k] = clight/(v_freqs[k]*1.e6);
}
virtual int getDataSignal(size_t j, vector<double> & sig)
{
sig.resize(v_time_data[j].size());
vector<double> & vdata = vv_data[j][nac_];
for(size_t k=0; k<sig.size(); k++) sig[k]=vdata[k];
return (int)sig.size();
}
virtual TVector<double> getDataSignal(size_t j)
{
vector<double> sig;
getDataSignal(j, sig);
TVector<double> rvs(sig.size());
for(size_t k=0; k<sig.size(); k++) rvs(k)=sig[k];
return rvs;
}
virtual int getExpectedSignal(size_t j, vector<double> & sig, double D_dish, double theta_b, double phi_b, double A, double B)
{
//DBG if (debug_gacfit > 0) cout<<"*DBG*A*getExpectedSignal() j="<<j<<" size="<<v_time_data[j].size()<<endl;
sig.resize(v_time_data[j].size());
vector<double> & vtime = v_time_data[j];
vector<double> & vdata = vv_data[j][nac_];
vector<double> & verr = vv_err[j][nac_];
for(size_t k=0; k<vtime.size(); k++) {
double elev=v_interp_elev[j].YInterp(vtime[k]);
double thetasrc=Angle(90.-elev,Angle::Degree).ToRadian();
double sazim=v_interp_sazim[j].YInterp(vtime[k]);
while (sazim > 360.) sazim -= 360.;
double phisrcdeg = 90.-sazim;
if (phisrcdeg < 0.) phisrcdeg+=360.;
// double phisrc=acos(v_interp_cphi[j].YInterp(vtime[k]));
// if (phisrc<0.) phisrc+=Angle::TwoPiCst();
ACBeam beam(D_dish, theta_b, phi_b, v_lambdas[j]);
beam.setGaussianLobe(fggauss_beam_);
double acb=beam.Value(thetasrc, Angle(phisrcdeg,Angle::Degree).ToRadian());
sig[k] = A*acb+B;
//DBG cout << " *DBG* "<<k<<" Theta,Phi-src = " << Angle(thetasrc).ToDegree() << " , " << Angle(phisrc).ToDegree()
//DBG << " elev="<<elev<<" phisrc="<<Angle(phisrc).ToDegree()<<" --> acb= "<<acb<< " sig="<<sig[k]<<endl;
}
return (int)sig.size();
}
virtual int getExpectedSignal(size_t j, vector<double> & sig, const double* parm)
{
/*
param[0] = D (Dish diameter)
parm[1] = delta_elev = theta
parm[2] = azim = phi
parm[3+2*j] = Aj (Amplitude) for track j
parm[4+2*j] = Bj (Baseline / Offset) for track j
*/
double D_dish = parm[0];
double theta_b = parm[1];
double phi_b = parm[2];
double A = parm[3+2*j];
double B = parm[4+2*j];
return getExpectedSignal(j, sig, D_dish, theta_b, phi_b, A, B);
}
virtual TVector<double> getExpectedSignal(size_t j, double D_dish, double theta_b, double phi_b, double A, double B)
{
vector<double> sig;
//DBG cout << "*DBG*getExpectedSignal() : calling getExpectedSignal(j="<<j<<" ,sig,...)"<<endl;
//DBG debug_gacfit=1;
getExpectedSignal(j, sig, D_dish, theta_b, phi_b, A, B);
TVector<double> rvs(sig.size());
for(size_t k=0; k<sig.size(); k++) rvs(k)=sig[k];
return rvs;
}
virtual double getXi2(const double* parm, int & npts)
{
double xi2=0.;
npts=0;
for(size_t j=0; j<v_time_data.size(); j++) {
vector<double> sig;
npts += getExpectedSignal(j, sig, parm);
vector<double> & vdata = vv_data[j][nac_];
vector<double> & verr = vv_err[j][nac_];
for(size_t k=0; k<vdata.size(); k++) {
double xx=(sig[k]-vdata[k])/verr[k];
xi2 += (xx*xx);
}
}
return xi2;
}
vector< vector<double> > & v_time_data;
vector< vector< vector<double> > > & vv_data;
vector< vector< vector<double> > > & vv_err;
vector<double> & v_freqs;
vector<double> v_lambdas;
vector< SLinInterp1D > & v_interp_elev;
vector< SLinInterp1D > & v_interp_sazim;
size_t nac_;
bool fggauss_beam_;
};
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment