Skip to content

GitLab

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

Codage de la gestion du flag NOAC ds le fit des auto-correlations, Reza 26/02/2019

parent f797d4b5
Hide whitespace changes
Inline Side-by-side
Showing with 77 additions and 30 deletions
Scripts/p4sat.cmd
View file @ adc5825d
......@@ -137,6 +137,15 @@ source /pbs/throng/baoradio/Library/Sophya/env.csh
cd /sps/baoradio/Reza/AnaP4/SatLD14jan19
set PRGBASE = /sps/baoradio/Reza/AnaP4/AnaPAON4
### Pour compiler le code au CC
cd $PRGBASE
git pull
make
## Pour compiler avec MINUIT
setenv MINUITDIR /pbs/throng/baoradio/Library/ExtLibs/pro/
make mn_trkacxfit
########################
cd /sps/baoradio/Reza/AnaP4/SatLD14jan19
## Extraction du signal pour les satellites entre 1275-1282 MHz
${PRGBASE}/Objs/tfm2dt ChezOlivier/LDOutputs/cor2_LD14jan19.ppf 1275,1282 dt_LD14jan19.ppf
......@@ -169,8 +178,12 @@ ${PRGBASE}/Satellites/Objs/trk2dt trk_40890_20190122.txt trk_40890_20190122.ppf
${PRGBASE}/Satellites/Objs/predictsatsgp4 -T "2019/01/22 14:35:00" -H 180.,55 -K 41859 TLE_20190125/gal*
${PRGBASE}/Satellites/Objs/trk2dt trk_41859_20190122.txt trk_41859_20190122.ppf 2019/01/14
## Extraction du signal pour les sources du ciel entre 1300-1310 MHz
${PRGBASE}/Objs/tfm2dt ChezOlivier/LDOutputs/cor2_LD14jan19.ppf 1300,1312 dtsky_LD14jan19.ppf
## Extraction du signal pour les sources du ciel entre 1342-1350 MHz -> FreqCentral=1346 MHz
${PRGBASE}/Objs/tfm2dt ../ChezOlivier/LDOutputs/cor2_LD14jan19.ppf 1342,1350 dtsky1346_LD14jan19.ppf
## Extraction du signal pour les sources du ciel entre 1392-1400 MHz -> FreqCentral=1396 MHz
${PRGBASE}/Objs/tfm2dt ../ChezOlivier/LDOutputs/cor2_LD14jan19.ppf 1392,1400 dtsky1396_LD14jan19.ppf
## Extraction du signal pour les sources du ciel entre 1410-1418 MHz -> FreqCentral=1416 MHz
${PRGBASE}/Objs/tfm2dt ../ChezOlivier/LDOutputs/cor2_LD14jan19.ppf 1410,1418 dtsky1416_LD14jan19.ppf
#### Fabrication des DataTables track a partir des fichiers texte fourni par JEC
foreach f ( CasA_survey_2019-1-14 CasA_survey_2019-1-15 CasA_survey_2019-1-16 CygA_survey_2019-1-17 CygA_survey_2019-1-18 )
......@@ -197,13 +210,18 @@ Sources avec observation vers CasA (z=+11.42)
@trk dt_LD14jan19 1115,1160 1278.5 trk_40889_20190114
@trk dt_LD14jan19 1460,1500 1278.5 trk_43566_20190115
@trk dt_LD14jan19 2455,2495 1278.5 trk_43055_20190115
@trk dtsky_LD14jan19 920,960 1305 CasA_survey_2019-1-14
@trk dtsky_LD14jan19 2355,2395 1305 CasA_survey_2019-1-15
@trk dtsky_LD14jan19 3790,3830 1305 CasA_survey_2019-1-16
trk dt_LD14jan19 920,960 1278.5 CasA_survey_2019-1-14
trk dt_LD14jan19 2355,2395 1278.5 CasA_survey_2019-1-15
trk dt_LD14jan19 3790,3830 1278.5 CasA_survey_2019-1-16
# Sky at 1346 MHz
@trk dtsky1346_LD14jan19 920,960 1346 CasA_survey_2019-1-14
@trk dtsky1346_LD14jan19 2355,2395 1346 CasA_survey_2019-1-15
@trk dtsky1346_LD14jan19 3790,3830 1346 CasA_survey_2019-1-16
# Sky at 1396 MHz
@trk dtsky1396_LD14jan19 920,960 1396 CasA_survey_2019-1-14
@trk dtsky1396_LD14jan19 2355,2395 1396 CasA_survey_2019-1-15
@trk dtsky1396_LD14jan19 3790,3830 1396 CasA_survey_2019-1-16
# Sky at 1416 MHz
@trk dtsky1416_LD14jan19 920,960 1416 CasA_survey_2019-1-14
@trk dtsky1416_LD14jan19 2355,2395 1416 CasA_survey_2019-1-15
@trk dtsky1416_LD14jan19 3790,3830 1416 CasA_survey_2019-1-16
#### ldtrk_cygA.d
--- Observation Longue duree Janvier 2019
......@@ -211,10 +229,12 @@ Sources avec observation vers CygA (z=-6.5)
@zenang -6.
@trk dt_LD14jan19 4130,4175 1278.5 trk_43057_20190116
@trk dt_LD14jan19 6120,6160 1278.5 trk_41175_20190118
@trk dtsky_LD14jan19 5025,5065 1305 CygA_survey_2019-1-17
@trk dtsky_LD14jan19 6460,6500 1305 CygA_survey_2019-1-18
trk dt_LD14jan19 5025,5065 1278.5 CygA_survey_2019-1-17
trk dt_LD14jan19 6460,6500 1278.5 CygA_survey_2019-1-18
@trk dtsky1346_LD14jan19 5025,5065 1346 CygA_survey_2019-1-17
@trk dtsky1346_LD14jan19 6460,6500 1346 CygA_survey_2019-1-18
@trk dtsky1396_LD14jan19 5025,5065 1396 CygA_survey_2019-1-17
@trk dtsky1396_LD14jan19 6460,6500 1396 CygA_survey_2019-1-18
@trk dtsky1416_LD14jan19 5025,5065 1416 CygA_survey_2019-1-17
@trk dtsky1416_LD14jan19 6460,6500 1416 CygA_survey_2019-1-18
#### ldtrk_virA.d
--- Observation Longue duree Janvier 2019
......@@ -230,6 +250,8 @@ ${PRGBASE}/Objs/trkacxfit -docx6s -docx6f -D 4.5 -out ldj.txt -ckf ldj.ppf ldtrk
# Ou
${PRGBASE}/Objs/trkacxfit -docx6s -docx6f -D 4.5 -out casA.txt -ckf casA.ppf ldtrk_casA.d
${PRGBASE}/Objs/trkacxfit -nocxf -D 4.5 -out casA.txt -ckf casA.ppf ldtrk_casA.d
#### Utiliser le script piapp : ckfitres.pic et les scripts definis dedans pour verifier les resultats
piapp> exec ckfitres.pic cxb6_casA.ppf
piapp> ck6 1 1
......
gacfit.h
View file @ adc5825d
......@@ -33,6 +33,15 @@ public:
v_lambdas[k] = clight/(v_freqs[k]*1.e6);
ndatapoints_=0;
for(size_t j=0; j<v_time_data.size(); j++) ndatapoints_+=v_time_data[j].size();
v_idxparam.resize(v_freqs.size());
int idx=0;
for(size_t j=0; j<v_time_data.size(); j++) {
if (v_noAC[j]) {
v_idxparam[j]=-1;
continue;
}
v_idxparam[j]=idx; idx++;
}
}
virtual TVector<double> getTimeVec(size_t j) const
......@@ -101,8 +110,13 @@ public:
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];
double A = 1.;
double B = 1.;
if (v_idxparam[j]>=0) {
int jj = v_idxparam[j];
double A = parm[3+2*jj];
double B = parm[4+2*jj];
}
return getExpectedSignal(j, sig, D_dish, theta_b, phi_b, A, B);
}
......@@ -122,7 +136,7 @@ public:
double xi2=0.;
npts=0;
for(size_t j=0; j<v_time_data.size(); j++) {
// if (v_noAC[j]) continue; On ne peut pas faire ca , il faut modifier sinon la liste des parametres
if (v_noAC[j]) continue; // On ne peut pas faire ca , il faut modifier sinon la liste des parametres
vector<double> sig;
npts += getExpectedSignal(j, sig, parm);
vector<double> & vdata = vv_data[j][nac_];
......@@ -140,7 +154,8 @@ public:
vector< vector< vector<double> > > & vv_err;
vector<double> & v_freqs;
vector<double> v_lambdas;
vector<bool> & v_noAC;
vector<bool> & v_noAC;
vector<int> v_idxparam; // vecteur de numero pour les parametres, prenant en compte les data/track non active ( v_noAC[j]=true )
vector< SLinInterp1D > & v_interp_elev;
vector< SLinInterp1D > & v_interp_sazim;
size_t nac_;
......
trkfit.cc
View file @ adc5825d
......@@ -466,21 +466,22 @@ int ACxSetFitter::doACfit(string outfilename)
phiAntenne=Angle(90.,Angle::Degree).ToRadian();
}
}
mFit.SetParam(1,"ThetaAntenne",thetaAntenne,M_PI/1440,0.,thetaAntenne+M_PI/36.);
mFit.SetParam(1,"ThetaAntenne",thetaAntenne,M_PI/1440,0.,M_PI/3.);
mFit.SetParam(2,"PhiAntenne",phiAntenne,M_PI/180.,0.,2.*M_PI);
// mFit.SetFix(1, thetaAntenne);
// mFit.SetFix(2, phiAntenne);
size_t jj=0;
for(size_t j=0; j<NTRK; j++) {
if (acxd_.v_noAC[j]) continue;
char pname[32];
sprintf(pname,"A%d",(int)(j+1));
double A = acxd_.v_max_data[j][ii];
mFit.SetParam(2*j+3,pname,A,A/10.,A/20,A*5);
mFit.SetParam(2*jj+3,pname,A,A/10.,A/20,A*5);
sprintf(pname,"B%d",(int)(j+1));
double B = acxd_.v_min_data[j][ii];
mFit.SetParam(2*j+4,pname,B,B/10.,B/20,B*5);
mFit.SetFix(2*j+4, B);
mFit.SetParam(2*jj+4,pname,B,B/10.,B/20,B*5);
// mFit.SetFix(2*jj+4, B);
}
//DBG mFit.PrintFit();
// cout << "do_p4_trkfit 2."<<ii+1<<" Performing the fit for AutoCor Antenna= " << ii+1 << endl;
......@@ -520,14 +521,23 @@ int ACxSetFitter::doACfit(string outfilename)
cout <<setw(16)<<"PhiAntenne= "<<setw(12)<<Angle(phiant).ToDegree()<< " +/- "
<<setw(12)<<Angle(err_phiant).ToDegree()<<" (azimuth ="
<<setw(8)<<azimdeg<<" +/- "<<setw(8)<<err_azimdeg<<" ) deg."<<endl;
ofr <<setw(8)<<azimdeg<<" "<<setw(8)<<err_azimdeg<<" ";
ofr <<setw(8)<<azimdeg<<" "<<setw(8)<<err_azimdeg<<" ";
jj=0;
for(size_t j=0; j<NTRK; j++) {
double A=mFit.GetParm(3+2*j); double err_A=mFit.GetParmErr(3+2*j);
double B=mFit.GetParm(4+2*j); double err_B=mFit.GetParmErr(4+2*j);
cout << " Trk/Sat["<<j<<"] -> A= "<<A<<" +/- "<<err_A<<" B= "<<B<<" +/- "<<err_B<<endl;
v_A[ii][j]=A; v_B[ii][j]=B;
v_err_A[ii][j]=err_A; v_err_B[ii][j]=err_B;
ofr <<setw(8)<<A<<" "<<setw(8)<<err_A<<" "<<setw(8)<<B<<" "<<setw(8)<<err_B<<" ";
double A = acxd_.v_max_data[j][ii];
double B = acxd_.v_min_data[j][ii];
double err_A = 0.;
double err_B = 0.;
if (!acxd_.v_noAC[j]) {
A=mFit.GetParm(3+2*jj); err_A=mFit.GetParmErr(3+2*jj);
B=mFit.GetParm(4+2*jj); err_B=mFit.GetParmErr(4+2*jj);
jj++;
}
cout << " Trk/Sat["<<j<<"] -> A= "<<A<<" +/- "<<err_A<<" B= "<<B<<" +/- "<<err_B<<(acxd_.v_noAC[j]?" NOT Fitted":"")<<endl;
if (acxd_.v_noAC[j])
ofr <<setw(8)<<A<<" "<<setw(8)<<" NOFIT "<<" "<<setw(8)<<B<<" "<<setw(8)<<" NOFIT "<<" ";
else
ofr <<setw(8)<<A<<" "<<setw(8)<<err_A<<" "<<setw(8)<<B<<" "<<setw(8)<<err_B<<" ";
}
ofr << endl;
double clight = PhysQty::c().SIValue();
......
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