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

Suite implementation / debug du fit des phases sur les cross-corr, Reza 21/12/2018

parent 0b361dd2
......@@ -158,7 +158,7 @@ public:
};
//-----------------------------------------------------------------
class MyACMinZ : public MinZFunction, MyACSignal {
class MyACMinZ : public MinZFunction, public MyACSignal {
public:
MyACMinZ(vector< vector<double> > & v_time_data_, vector< vector< vector<double> > > & vv_data_,
vector< vector< vector<double> > > & vv_err_, vector<double> & v_freqs_,
......
......@@ -32,7 +32,7 @@ public:
return (int)sig.size();
}
virtual TVector<double> getDataSignal(size_t j)
virtual TVector< complex<double> > getDataSignal(size_t j)
{
vector< complex<double> >sig;
getDataSignal(j, sig);
......@@ -42,7 +42,7 @@ public:
}
virtual int getExpectedSignal(size_t j, vector< complex<double> > & sig, double phase, double A)
virtual int getExpectedSignal(size_t j, vector< complex<double> > & sig, double phase, double A, complex<double> B=complex<double>(0.,0.))
{
//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());
......@@ -68,11 +68,12 @@ public:
virtual int getExpectedSignal(size_t j, vector< complex<double> > & sig, const double* parm)
{
double phase = parm[0];
double A = parm[1+j];
double A = parm[1+3*j];
complex<double> B = complex<double>(parm[2+3*j], parm[3+3*j]);
return getExpectedSignal(j, sig, phase, A);
}
virtual TVector< complex<double> > getExpectedSignal(size_t j, double phase, double A)
virtual TVector< complex<double> > getExpectedSignal(size_t j, double phase, double A, complex<double> B=complex<double>(0.,0.))
{
vector< complex<double> > sig;
//DBG cout << "*DBG*getExpectedSignal() : calling getExpectedSignal(j="<<j<<" ,sig,...)"<<endl;
......@@ -113,13 +114,13 @@ public:
};
//-----------------------------------------------------------------
class MyCxGenXi2 : public GeneralXi2, MyCxSignal {
class MyCxGenXi2 : public GeneralXi2, public MyCxSignal {
public:
MyCxGenXi2(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_,
CxBeam & beam, size_t ncx=0)
: GeneralXi2(1+v_time_data_.size()) ,
: GeneralXi2(1+3*v_time_data_.size()) ,
MyCxSignal(v_time_data_, vv_cxdata_, vv_cxerr_, v_interp_elev_, v_interp_sazim_, beam, ncx)
{
}
......@@ -127,7 +128,9 @@ public:
virtual double Value(GeneralFitData& data, double* parm, int& ndataused)
/*
param[0] = phase
parm[1+j] = Aj (Amplitude) for track j
parm[1+3*j] = Aj (Amplitude) for track j
parm[2+3*j] = Bj-real Offset, real part
parm[3+3*j] = Bj-imag Offset, imag part
*/
{
return getXi2(parm, ndataused);
......
......@@ -629,6 +629,7 @@ double P4Coords::getLat(){return(paon_lat);}
Vector3d P4Coords::getAntennaPosition(sa_size_t num)
{
if ((num<1)||(num>4)) throw ParmError("P4Coords::getAntennaPosition() out of range antenna position <1 or>4 !");
num--;
return Vector3d(paon4_posx[num],paon4_posy[num],paon4_posz[num]);
}
......@@ -636,6 +637,7 @@ Vector3d P4Coords::getBaseline(sa_size_t num1, sa_size_t num2)
{
if ((num1<1)||(num1>4)) throw ParmError("P4Coords::getBaseline() out of range antenna position num1 <1 or>4 !");
if ((num2<1)||(num2>4)) throw ParmError("P4Coords::getBaseline() out of range antenna position num2 <1 or>4 !");
num1--; num2--;
return Vector3d(paon4_posx[num2]-paon4_posx[num1],paon4_posy[num2]-paon4_posy[num1],paon4_posz[num2]-paon4_posz[num1]);
}
......
......@@ -261,9 +261,9 @@ public :
static double getLon();
//! static method returning PAON4 (central antenna) latitude in decimal degrees (positive -> north)
static double getLat();
//! return the antenna position for antenna \b num (1..4) in the PAON4 coordinate system
//! return the antenna position for antenna \b num (1 <= num <= 4) in the PAON4 coordinate system
static Vector3d getAntennaPosition(sa_size_t num);
//! return the antenna pair baseline from num1 -> num2
//! return the antenna pair baseline from num1 -> num2 (antenna identification numbers 1 <= num1, num2 <= 4)
static Vector3d getBaseline(sa_size_t num1, sa_size_t num2);
//! static method Converting a TU time into RA (Right Ascension sky coordinates) in decimal degrees (uses old EROS code in SOPHYA, datime.h)
static double RAFromTimeTU(TimeStamp const & tu);
......
......@@ -74,6 +74,7 @@ static vector<double> tstart, tend; // time interval in seconds (start,end tim
static vector<double> v_freqs; // Central frequency associated to each visibilty = f(time) dataset
static vector<string> trkflnm; // sources (satellites, sky sources ...) track datatable file names
static double zenang=0.; // zenith angle = shift with declination, used for initial value of fitted angles
static double theta_0=0., phi_0=0.; // Theta, phi angles corresponding to the theoretical pointing defined by zenang
//--- End of list of global parameters
......@@ -82,7 +83,8 @@ int decode_args(int narg, char* arg[]);
int read_p4_autocxcors(vector< vector<double> > & v_time_data, vector< vector< vector<double> > > & vv_data,
vector< vector< vector<double> > > & vv_err,
vector< vector<double> > & v_min_data, vector< vector<double> > & v_max_data,
vector< vector< vector< complex<double> > > > & vv_cxdata, vector< vector< vector<double> > > & vv_cxerr);
vector< vector< vector< complex<double> > > > & vv_cxdata, vector< vector< vector<double> > > & vv_cxerr,
vector< vector<double> > & v_min_cxdata, vector< vector<double> > & v_max_cxdata);
int read_srctracks(vector< vector<double> > & v_time_sat, vector< vector<double> > & v_sat_elev,
vector< vector<double> > & v_sat_azim, vector< SLinInterp1D > & v_interp_elev,
......@@ -96,7 +98,8 @@ int do_p4_trkfit(vector< vector<double> > & v_time_data, vector< vector< vector<
vector< vector<double> > & v_A, vector< vector<double> > & v_B);
int do_p4_cxtrkfit(vector< vector<double> > & v_time_data, vector< vector< vector< complex<double> > > > & vv_cxdata,
vector< vector< vector<double> > > & vv_cxerr, vector<double> & v_freqs,
vector< vector< vector<double> > > & vv_cxerr, vector<double> & v_freqs,
vector< vector<double> > & v_min_cxdata, vector< vector<double> > & v_max_cxdata,
vector< SLinInterp1D > & v_interp_elev, vector< SLinInterp1D > & v_interp_sazim,
vector<double> & v_Ddish, vector<double> & v_thetaant, vector<double> & v_phiant,
vector<double> & v_phase, vector< vector<double> > & v_A);
......@@ -123,8 +126,10 @@ int main (int narg, char* arg[])
vector< vector<double> > v_max_data(NTRK);
vector< vector< vector< complex<double> > > > vv_cxdata(NTRK);
vector< vector< vector<double> > > vv_cxerr(NTRK);
vector< vector<double> > v_min_cxdata(NTRK);
vector< vector<double> > v_max_cxdata(NTRK);
int tot_npoints = read_p4_autocxcors(v_time_data, vv_data, vv_err, v_min_data, v_max_data, vv_cxdata, vv_cxerr);
int tot_npoints = read_p4_autocxcors(v_time_data, vv_data, vv_err, v_min_data, v_max_data, vv_cxdata, vv_cxerr, v_min_cxdata, v_max_cxdata);
vector< vector<double> > v_time_sat(NTRK);
vector< vector<double> > v_sat_elev(NTRK);
......@@ -167,9 +172,10 @@ int main (int narg, char* arg[])
}
} //End-of-check , saving data & expected signal to outcheckfilename
vector<double> v_phase;
int rcf2 = do_p4_cxtrkfit(v_time_data, vv_cxdata, vv_cxerr, v_freqs,
v_interp_elev, v_interp_sazim, v_Ddish, v_thetaant, v_phiant, v_phase, v_A);
vector<double> v_phase(NB_CXCORS);
vector< vector<double> > v_Acx(NB_CXCORS);
int rcf2 = do_p4_cxtrkfit(v_time_data, vv_cxdata, vv_cxerr, v_freqs, v_min_cxdata, v_max_cxdata,
v_interp_elev, v_interp_sazim, v_Ddish, v_thetaant, v_phiant, v_phase, v_Acx);
} // End of try bloc
catch (PThrowable & exc) {
......@@ -198,7 +204,9 @@ int main (int narg, char* arg[])
int read_p4_autocxcors(vector< vector<double> > & v_time_data, vector< vector< vector<double> > > & vv_data,
vector< vector< vector<double> > > & vv_err,
vector< vector<double> > & v_min_data, vector< vector<double> > & v_max_data,
vector< vector< vector< complex<double> > > > & vv_cxdata, vector< vector< vector<double> > > & vv_cxerr)
vector< vector< vector< complex<double> > > > & vv_cxdata, vector< vector< vector<double> > > & vv_cxerr,
vector< vector<double> > & v_min_cxdata, vector< vector<double> > & v_max_cxdata)
{
int tot_npoints = 0; // total number of points for fit
cout << "---- trkacfit/read_p4_autocxcors() ; reading 4 PAON4 auto-correlation & 6 Cross-cor signals/DataTables ..."<<endl;
......@@ -234,6 +242,8 @@ int read_p4_autocxcors(vector< vector<double> > & v_time_data, vector< vector< v
vector<double> vetmp;
vv_cxdata[j].push_back(vtmp);
vv_cxerr[j].push_back(vetmp);
v_min_cxdata[j].push_back(9.e19);
v_max_cxdata[j].push_back(-9.e19);
}
vector< vector<double> > & v_data = vv_data[j];
......@@ -254,16 +264,22 @@ int read_p4_autocxcors(vector< vector<double> > & v_time_data, vector< vector< v
for(size_t ii=0; ii<NB_CXCORS; ii++) { // 6 cross-correlations
vector< complex<double> > & vcx = v_vcx[ii];
v_cxdata[ii].push_back(vcx[k]);
v_cxerr[ii].push_back(0.1*std::abs(vcx[k]));
double acx=std::abs(vcx[k]);
v_cxerr[ii].push_back(0.1*acx);
if (acx<v_min_cxdata[j][ii]) v_min_cxdata[j][ii]=acx;
if (acx>v_max_cxdata[j][ii]) v_max_cxdata[j][ii]=acx;
}
}
tot_npoints += v_time_data[j].size(); // total number of points for fit
cout << " ... Done for " << j+1 << " data size="<<v_time_data[j].size()<<endl;
cout << " Data-AutoCor Min[A1...A4]=";
for(size_t ii=0; ii<4; ii++) cout<<setw(10)<<v_min_data[j][ii]<< " , "; cout << endl;
cout << " Data-AutoCor Max[A1...A4]=";
for(size_t ii=0; ii<4; ii++) cout<<setw(10)<<v_max_data[j][ii]<< " , "; cout << endl;
cout << " Data-AutoCor Min,Max[A1...A4]=";
for(size_t ii=0; ii<NB_ANTENNES; ii++)
cout<<setw(10)<<v_min_data[j][ii]<<","<<setw(10)<<v_max_data[j][ii]<<" ; "; cout << endl;
cout << " Data-CxCorr (abs) Min,Max[Cx1...Cx6]=";
for(size_t ii=0; ii<NB_ANTENNES; ii++)
cout<<setw(10)<<v_min_cxdata[j][ii]<<","<<setw(10)<<v_max_cxdata[j][ii]<<" ; "; cout << endl;
}
return tot_npoints;
}
......@@ -330,6 +346,8 @@ int do_p4_trkfit(vector< vector<double> > & v_time_data, vector< vector< vector<
vector<double> & v_Ddish, vector<double> & v_thetaant, vector<double> & v_phiant,
vector< vector<double> > & v_A, vector< vector<double> > & v_B)
{
cout << "--------------------------------------------------------------------------------------"<<endl;
cout << "--------------------------------------------------------------------------------------"<<endl;
cout << "---- trkacfit/do_p4_trkfit() ; Performing antenna pointing fit ..."<<endl;
ofstream ofr(outfilename.c_str());
ofr << "#### Pointing/dish diameter fit on autocorrelation (programme trkacfit.cc ) "<<endl
......@@ -340,7 +358,7 @@ int do_p4_trkfit(vector< vector<double> > & v_time_data, vector< vector< vector<
int tot_npoints = 0;
for(size_t j=0; j<NTRK; j++) tot_npoints += v_time_data[j].size();
for(size_t ii=0; ii<NB_ANTENNES; ii++) {
cout << "3."<<ii+1<<" Creating General Fit for AutoCor Antenna= " << ii+1 << endl;
cout << "do_p4_trkfit 1."<<ii+1<<" Creating General Fit for AutoCor Antenna= " << ii+1 << endl;
GeneralFitData gdata(1, tot_npoints);
for(size_t j=0; j<NTRK; j++) {
vector< vector<double> > & v_data = vv_data[j];
......@@ -385,7 +403,7 @@ int do_p4_trkfit(vector< vector<double> > & v_time_data, vector< vector< vector<
}
//DBG mFit.PrintFit();
cout << "4."<<ii+1<<" Performing the fit for AutoCor Antenna= " << ii+1 << endl;
cout << "do_p4_trkfit 2."<<ii+1<<" Performing the fit for AutoCor Antenna= " << ii+1 << endl;
int rcfit = mFit.Fit();
if (prtlevel>1) mFit.PrintFit();
if(rcfit>0) {
......@@ -434,12 +452,18 @@ int do_p4_trkfit(vector< vector<double> > & v_time_data, vector< vector< vector<
/* --Fonction-- */
int do_p4_cxtrkfit(vector< vector<double> > & v_time_data, vector< vector< vector< complex<double> > > > & vv_cxdata,
vector< vector< vector<double> > > & vv_cxerr, vector<double> & v_freqs,
vector< vector< vector<double> > > & vv_cxerr, vector<double> & v_freqs,
vector< vector<double> > & v_min_cxdata, vector< vector<double> > & v_max_cxdata,
vector< SLinInterp1D > & v_interp_elev, vector< SLinInterp1D > & v_interp_sazim,
vector<double> & v_Ddish, vector<double> & v_thetaant, vector<double> & v_phiant,
vector<double> & v_phase, vector< vector<double> > & v_A)
{
cout << "--------------------------------------------------------------------------------------"<<endl;
cout << "--------------------------------------------------------------------------------------"<<endl;
cout << "---- trkacfit/do_p4_cxtrkfit() ; Performing cross-cor phase fit ..."<<endl;
POutPersist pox("cxdbg.ppf");
ofstream ofr("cxcorfit.txt");
ofr << "#### cross-cor phase fit (programme trkacfit.cc ) "<<endl
<< "## NumCxCor Xi2red Phase err_Phase A0 err_A0 A1 err_A1 ..."<<endl;
......@@ -448,37 +472,98 @@ int do_p4_cxtrkfit(vector< vector<double> > & v_time_data, vector< vector< vecto
size_t Anum1[6]={0,0,0,1,1,2};
size_t Anum2[6]={1,2,3,2,3,3};
for(size_t ii=0; ii<NB_CXCORS; ii++) {
cout << "3."<<ii+1<<" Creating General Fit for AutoCor Antenna= " << ii+1 << endl;
Vector3d baseline=P4Coords::getBaseline(Anum1[ii]+1,Anum2[ii]+1);
cout << "do_p4_cxtrkfit() 1."<<ii+1<<" Creating General Fit for CrossCor= " << ii << " FxF= " << Anum1[ii]+1<<"x"<<Anum2[ii]+1
<<" Baseline="<<baseline<<endl;
GeneralFitData gdata(1, tot_npoints);
/* for(size_t j=0; j<NTRK; j++) {
vector< vector<double> > & v_data = vv_data[j];
vector< vector<double> > & v_err = vv_err[j];
for(size_t j=0; j<NTRK; j++) {
vector< vector< complex<double> > > & v_cxdata = vv_cxdata[j];
vector< vector<double> > & v_cxerr = vv_cxerr[j];
for(size_t k=0; k<v_time_data[j].size(); k++) {
gdata.AddData1(v_time_data[j][k],v_data[ii][k],v_err[ii][k]); // Fill x, y and error on y
gdata.AddData1(v_time_data[j][k],v_cxdata[ii][k].real(),v_cxerr[ii][k]); // Fill x, y and error on y
gdata.AddData1(v_time_data[j][k],v_cxdata[ii][k].imag(),v_cxerr[ii][k]); // Fill x, y and error on y
}
} */
}
double clight = PhysQty::c().SIValue();
double lambda = clight/(v_freqs[0]*1.e6);
ACBeam acb1(v_Ddish[Anum1[ii]], v_thetaant[Anum1[ii]], v_phiant[Anum1[ii]], lambda);
acb1.setGaussianLobe(fggaussbeam);
ACBeam acb2(v_Ddish[Anum2[ii]], v_thetaant[Anum2[ii]], v_phiant[Anum2[ii]], lambda);
Vector3d baseline=P4Coords::getBaseline(Anum1[ii]+1,Anum2[ii]+1);
acb2.setGaussianLobe(fggaussbeam);
CxBeam cxbeam(acb1, acb2, baseline);
/*
cout << "do_p4_cxtrkfit()*DBG*A0* Theta_0= "<<theta_0<<" Phi_0= "<<phi_0<<endl;
cout << "do_p4_cxtrkfit()*DBG*A1* Theta1="<<v_thetaant[Anum1[ii]]<<" Phi1="<<v_phiant[Anum1[ii]]<< " Beam1="<<acb1.Value(theta_0, phi_0)<<endl;
cout << "do_p4_cxtrkfit()*DBG*A2* Theta2="<<v_thetaant[Anum2[ii]]<<" Phi2="<<v_phiant[Anum2[ii]]<< " Beam2="<<acb2.Value(theta_0, phi_0)<<endl;
*/
MyCxGenXi2 gxi2( v_time_data, vv_cxdata, vv_cxerr, v_interp_elev, v_interp_sazim, cxbeam, ii);
GeneralFit mFit(&gxi2);
mFit.SetData(&gdata); // connect data to the fitter , here the data is unused - gxi2 includes its data
mFit.SetMaxStep(1000);
// SetParam(int n,double value, double step,double min=1., double max=-1.);
mFit.SetParam(0,"Phase",M_PI,M_PI/360.,0.,2.2*M_PI);
mFit.SetParam(0,"Phase",0.,M_PI/360.,0.,2.2*M_PI);
char oname[32];
vector<double> v_amp(NTRK);
for(size_t j=0; j<NTRK; j++) {
double A=1.; // v_max_cxdata[j][ii];
TVector< complex<double> > signal = gxi2.getDataSignal(j);
Vector asig = SOPHYA::abs(signal);
double mins, maxs;
asig.MinMax(mins, maxs);
sprintf(oname,"cx_%d_%d",(int)ii+1,(int)j+1);
pox << PPFNameTag(oname)<<signal;
TVector< complex<double> > expsignal = gxi2.getExpectedSignal(j, 0., A);
Vector aexpsig = SOPHYA::abs(expsignal);
double mine, maxe;
aexpsig.MinMax(mine, maxe);
A=maxs/maxe;
v_amp[j]=A;
expsignal *= complex<double>(A, 0.);
//DBG cout << "do_p4_cxtrkfit()*DBG*B* j="<<j<<" A="<<A<<" max-sig= "<<maxs<<" max-exp="<<maxe<<endl;
sprintf(oname,"simcx_%d_%d",(int)ii+1,(int)j+1);
pox << PPFNameTag(oname)<<expsignal;
}
double fparm[50]; fparm[0]=0.;
double bestxi2 = 9.e19;
double bestphase=0.;
int bestnpts,npts;
int bestafact;
double afact[12]={0.5,0.6,0.7,0.8,0.9,1.0,1.1,1.2,1.35,1.5,1.75,2.0};
for(int ia=0; ia<12; ia++) {
for(size_t j=0; j<NTRK; j++) {
fparm[1+3*j]=afact[ia]*v_amp[j]; fparm[2+3*j]=fparm[3+3*j]=0.;
}
for(double ph=0.; ph<360.; ph += 1) {
fparm[0]=Angle(ph, Angle::Degree).ToRadian();
double xi2 = gxi2.getXi2(fparm, npts);
if (xi2 < bestxi2) {
bestxi2 = xi2; bestphase=fparm[0]; bestnpts=npts; bestafact=afact[ia];
}
}
}
mFit.SetParam(0,"Phase",bestphase,M_PI/360.,0.,2.2*M_PI);
cout << "do_p4_cxtrkfit() 2."<<ii+1<<" bestxi2="<<bestxi2<<" bestphase="<<Angle(bestphase).ToDegree()<<" bestnpts="<<bestnpts
<<" bestafact="<<bestafact<< " A= "; for(size_t j=0; j<NTRK; j++) cout << v_amp[j] << " , "; cout << endl;
for(size_t j=0; j<NTRK; j++) {
char pname[32];
sprintf(pname,"A%d",(int)(j+1));
double A = sqrt(v_A[j][Anum1[ii]] * v_A[j][Anum2[ii]]);;
mFit.SetParam(1+j,pname,A,A/10.,A/2,A*2);
sprintf(pname,"A%d",(int)(j+1));
double A=v_amp[j]; // v_max_cxdata[j][ii];
// double A = sqrt(v_A[j][Anum1[ii]] * v_A[j][Anum2[ii]]);
mFit.SetParam(1+3*j,pname,A,A/10.,A/4,A*3);
sprintf(pname,"Bre%d",(int)(j+1));
mFit.SetParam(2+3*j,pname,0.,A/25.,-A/5,A/5.);
sprintf(pname,"Bim%d",(int)(j+1));
mFit.SetParam(3+3*j,pname,0.,A/25.,-A/5,A/5.);
mFit.SetFix(2+3*j,0.);
mFit.SetFix(3+3*j,0.);
}
//DBG mFit.PrintFit();
cout << "4."<<ii+1<<" Performing the fit for CrossCor " << ii << endl;
cout << "do_p4_cxtrkfit() 3."<<ii+1<<" Performing the fit for CrossCor " << ii << " FxF= " << Anum1[ii]+1<<"x"<<Anum2[ii]+1<<endl;
int rcfit = mFit.Fit();
if (prtlevel>1) mFit.PrintFit();
if(rcfit>0) {
......@@ -486,11 +571,11 @@ int do_p4_cxtrkfit(vector< vector<double> > & v_time_data, vector< vector< vecto
cout<< "---Fit resultReduce_Chisquare = " << mFit.GetChi2Red()<< " nstep="<<mFit.GetNStep() << " rc="<<rcfit<<endl;
ofr <<setw(4)<<ii+1<<" "<<setw(8)<<mFit.GetChi2Red()<<" ";
double phase=mFit.GetParm(0); double err_phase=mFit.GetParmErr(0);
cout <<setw(16)<<"Phase= "<<setw(10)<<phase<<" +/- "<<setw(10)<<err_phase<<" m."<<endl;
cout <<setw(16)<<"Phase= "<<setw(10)<<Angle(phase).ToDegree()<<" +/- "<<setw(10)<<Angle(err_phase).ToDegree()<<" deg."<<endl;
ofr <<setw(8)<<phase<<" "<<setw(8)<<err_phase<<" ";
// v_Ddish[ii]=Dfit;
for(size_t j=0; j<NTRK; j++) {
double A=mFit.GetParm(1+j); double err_A=mFit.GetParmErr(1+j);
double A=mFit.GetParm(1+3*j); double err_A=mFit.GetParmErr(1+3*j);
cout << " Trk/Sat["<<j<<"] -> A= "<<A<<" +/- "<<err_A<<endl;
// v_A[ii][j]=A; v_B[ii][j]=B;
// ofr <<setw(8)<<A<<" "<<setw(8)<<err_A<<" "<<setw(8)<<B<<" "<<setw(8)<<err_B<<" ";
......@@ -592,7 +677,13 @@ int decode_args(int narg, char* arg[])
cout << " Beam: D_dish (initial guess)="<<D_dish<<" GaussianBeam ? "<<(fggaussbeam?"Yes":"No")
<< " dec-shift(zenithAngle)= "<<zenang<<endl;
cout << " Input Base Path= "<< input_base_path << " NTrack="<<NTRK<<" OutFileName= "<<outfilename<<endl;
if (zenang<0.) {
theta_0 = Angle(-zenang, Angle::Degree).ToRadian(); phi_0 = Angle::PioTwoCst()+Angle::OnePiCst();
}
else {
theta_0 = Angle(+zenang, Angle::Degree).ToRadian(); phi_0 = Angle::PioTwoCst();
}
NTRK=lastargs.size();
dataflnm.resize(NTRK);
tstart.resize(NTRK); tend.resize(NTRK);
......
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