Commit 4d013c7c authored by Reza  ANSARI's avatar Reza ANSARI
Browse files

Suite debug fit simultane des 6 cx-cor, Reza 21/02/2019

parent f030f383
......@@ -15,13 +15,25 @@
class My6CxSignalsB {
public:
My6CxSignalsB(vector< AcxDataSet> & acxd, vector< TrackSet > & trks)
: n_cxcor(6), acxd_(acxd), trk_(trks)
: acxd_(acxd), trk_(trks), prtlevel_(0)
{
if (acxd_.size() != trk_.size())
throw ParmError("My6CxSignalsB::My6CxSignalsB(acxd,trks)/ERROR acxd_.size() != trk_.size()");
if (acxd_.size() < 1)
throw ParmError("My6CxSignalsB::My6CxSignalsB(acxd,trks)/ERROR acxd_.size() < 1");
n_cxcor=acxd_[0].getNbCrossCor();
v_xi2.resize(acxd_[0].getNbCrossCor());
v_npts_xi2.resize(acxd_[0].getNbCrossCor());
for(size_t kcx=0; kcx<n_cxcor; kcx++) {
v_npts_xi2[kcx]=0;
for(size_t i=0; i<acxd_.size(); i++)
v_npts_xi2[kcx]+=acxd_[i].v_time_data.size();
}
//DBG cout << " My6CxSignalsB() "<< endl;
}
inline void SetPrintLevel(int lev=0) { prtlevel_=lev; }
// get the the TimeVec in dataset i , corresponding to the j-th track
virtual TVector<double> getTimeVec(size_t i, size_t j)
{
......@@ -64,12 +76,15 @@ public:
sig.resize(vtime.size());
// vector< complex<double> > & vcxdata = vv_cxdata[j][nac_];
// vector<double> & vcxerr = vv_cxerr[j][nac_];
double A = acxd_[i].v_Acx[j][kcx];
complex<double> B = acxd_[i].v_Bcx[j][kcx];
double A = acxd_[i].v_Acx[kcx][j];
complex<double> B = acxd_[i].v_Bcx[kcx][j];
CxBeam beam = getBeam(i, kcx);
beam.ShiftBaseline(baseline_shift);
if (prtlevel_>1) {
cout << "My6CxSignalsB::getExpectedSignal() A="<<A<<" beam.baseline="<<beam.baseline_
<< " a1.D="<<beam.a1_.getDdish()<<endl;
}
for(size_t l=0; l<vtime.size(); l++) {
double elev=trk_[i].v_interp_elev[j].YInterp(vtime[l]);
double thetasrc=Angle(90.-elev,Angle::Degree).ToRadian();
......@@ -110,15 +125,23 @@ public:
size_t na0[6]={0,0,0,0,0,1};
size_t nb0[6]={0,1,2,1,2,2};
size_t ja=0, jb=0;
if (kcx<3) {
phase=parm[kcx];
baseline_shift=Vector3d(parm[3*kcx+3], parm[3*kcx+4], parm[3*kcx+5]);
}
else {
size_t ja=na0[kcx] , jb=nb0[kcx];
ja=na0[kcx] , jb=nb0[kcx];
phase=parm[jb]-parm[ja];
baseline_shift=Vector3d(parm[3*jb+3]-parm[3*ja+3], parm[3*jb+4]-parm[3*ja+4], parm[3*jb+5]-parm[3*ja+5]);
}
if (prtlevel_>1) {
cout << "My6CxSignalsB::getFromFitParam(kcx="<<kcx<<") ja="<<ja<<" jb="<<jb
<<" ->phase="<<phase<<" BaselineShift="<<baseline_shift<<endl<<" parm[0..11]= ";
for(size_t k=0; k<12; k++) cout << parm[k]<<" ; ";
cout<<endl;
}
return;
}
// get the expected signal for trackset i , corresponding to the j-th track for cross-correlation kcx (0<=kcx<=5)
......@@ -156,21 +179,34 @@ public:
{
double xi2=0.;
npts=0;
for(size_t i=0; i<acxd_.size(); i++) {
for(size_t j=0; j<acxd_[i].v_time_data.size(); j++) { // loop over tracks in the data set
for(size_t kcx=0; kcx<n_cxcor; kcx++) { // loop over the six cross-correlations
double myphase=0.;
Vector3d myblshift;
for(size_t kcx=0; kcx<n_cxcor; kcx++) { // loop over the six cross-correlations
v_xi2[kcx]=0.;
double myxi2=0.;
getFromFitParam(kcx, parm, myphase, myblshift);
for(size_t i=0; i<acxd_.size(); i++) {
for(size_t j=0; j<acxd_[i].v_time_data.size(); j++) { // loop over tracks in the data set
vector< complex<double> > sig;
npts += getExpectedSignal(i, j, kcx, sig, parm);
npts += getExpectedSignal(i, j, kcx, sig, myphase, myblshift);
vector< complex<double> > & vcxdata = acxd_[i].vv_cxdata[j][kcx];
vector<double> & verr = acxd_[i].vv_cxerr[j][kcx];
for(size_t l=0; l<vcxdata.size(); l++) {
double xx=(sig[l].real()-vcxdata[l].real())/verr[l];
xi2 += (xx*xx);
myxi2 += (xx*xx);
xx=(sig[l].imag()-vcxdata[l].imag())/verr[l];
xi2 += (xx*xx);
myxi2 += (xx*xx);
}
}
}
v_xi2[kcx] = myxi2/(double)(2*v_npts_xi2[kcx]);
xi2 += myxi2;
}
if (prtlevel_>1) {
cout << "My6CxSignalsB::getXi2() : xi2[1...6]= ";
for(size_t kcx=0; kcx<n_cxcor; kcx++) cout << v_xi2[kcx]<<" ; ";
cout<<endl;
}
npts *= 2;
return xi2;
......@@ -179,6 +215,9 @@ public:
size_t n_cxcor; // Number of Cx correlations ( 6 for PAON4 )
vector< AcxDataSet> & acxd_;
vector< TrackSet > & trk_;
vector< double > v_xi2;
vector< size_t > v_npts_xi2;
int prtlevel_;
};
//-----------------------------------------------------------------
......@@ -211,5 +250,20 @@ public:
}
};
//-----------------------------------------------------------------
class My6CxMinZFunc : public MinZFunction, public My6CxSignalsB {
public:
My6CxMinZFunc(vector< AcxDataSet> & acxd, vector< TrackSet > & trks)
: MinZFunction(12) , // 12 = 3 phases + 3*3 position shifts
My6CxSignalsB(acxd, trks)
{
}
virtual double Value(double const parm[])
{
int ndataused;
double xi2=getXi2(parm, ndataused);
return (xi2/(double)ndataused);
}
};
#endif
......@@ -142,6 +142,9 @@ int main (int narg, char* arg[])
CxBaselineFitter cxbfit(v_acxd, v_trk);
string cxbofile = "cxb6_"+outfilename;
string cxbckfile = "cxb6_"+checkfilename;
//cxbfit.doSimplexMinimize();
cxbfit.initFitParams();
// cxbfit.doCheck();
cxbfit.dofit(cxbofile,fg_fixbaseline);
cxbfit.saveExpectedCx(cxbckfile);
}
......
......@@ -540,14 +540,9 @@ int ACxSetFitter::doCxfit(string outfilenamecx, bool useAac)
cout << "======================================================================================"<<endl;
cout << "---------- ACxSetFitter::doCxfit() ; Performing cross-cor phase fit ..."<<endl;
if (useAac) cout << " ... Using Amplitude from auto-correlations fit for initial fit parameter value..."<<endl;
/* DEL POutPersist * pox = NULL ;
if (outcheckfilenamecx.length()>0) {
cout << "... expected cross-cor for fitted params (and cadata) will be saved to file "<<outcheckfilenamecx<<endl;
pox = new POutPersist(outcheckfilenamecx);
} */
ofstream ofr(outfilenamecx.c_str());
ofr << "#### cross-cor phase fit (ACxSetFitter::doCxfit() ) "<<endl
<< "## NumCxCor Xi2red Phase err_Phase (deg) A0 err_A0 A1 err_A1 ..."<<endl;
<< "## NumCxCor RcFit Xi2red Phase err_Phase (deg) A0 err_A0 A1 err_A1 ..."<<endl;
int tot_npoints_fit = 0;
for(size_t j=0; j<NTRK; j++) tot_npoints_fit += 2*(acxd_.v_time_data[j].size());
size_t Anum1[6]={0,0,0,1,1,2};
......@@ -598,21 +593,12 @@ int ACxSetFitter::doCxfit(string outfilenamecx, bool useAac)
Vector asig = SOPHYA::abs(signal);
double mins, maxs;
asig.MinMax(mins, maxs);
/* if (pox) {
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;
/*DBG if (pox) {
expsignal *= complex<double>(A,0.);
sprintf(oname,"ecx_%d_%d",(int)ii+1,(int)j+1);
(*pox) << PPFNameTag(oname)<<expsignal;
} */
}
double fparm[50]; fparm[0]=0.;
......@@ -669,42 +655,36 @@ int ACxSetFitter::doCxfit(string outfilenamecx, bool useAac)
if(rcfit>0) {
v_xi2red_cx[ii]=mFit.GetChi2Red();
// cout<< "-------------------------- Result for Cross No " << ii << endl;
cout<< "------ Fit result for Cross No "<<ii<<" Reduce_Chisquare = " << mFit.GetChi2Red()
cout<< "------ Fit result for Cross No "<<ii+1<<" Reduce_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);
if (phase<0.) phase += 2.*M_PI;
if (phase>2.*M_PI) phase -= 2.*M_PI;
cout <<"Phase= "<<setw(10)<<Angle(phase).ToDegree()<<" +/- "<<setw(10)<<Angle(err_phase).ToDegree()<<" deg."<<endl;
ofr <<setw(8)<<Angle(phase).ToDegree()<<" "<<setw(8)<<Angle(err_phase).ToDegree()<<" ";
v_phase[ii]=phase;
v_err_phase[ii]=err_phase;
for(size_t j=0; j<NTRK; j++) {
double Aac=sqrt(v_A[Anum1[ii]][j] * v_A[Anum2[ii]][j]);
double Ai=(useAac?Aac:v_amp[j]);
double A=mFit.GetParm(1+3*j); double err_A=mFit.GetParmErr(1+3*j);
cout << " Trk["<<j<<"] A= "<<A<<" +/- "<<err_A<<" (A/Ai="<<A/Ai<<")"<<endl;
v_Acx[ii][j]=A;
v_Bcx[ii][j]=complex<double>(0.,0.);
v_err_Acx[ii][j]=err_A;
ofr <<setw(8)<<A<<" "<<setw(8)<<err_A<<" ";
}
ofr << endl;
}
else {
cout << "---Fit failed for "<<ii<<" Fit_Error, rc = " << rcfit << " nstep="<<mFit.GetNStep()<<endl;
ofr <<setw(4)<<ii+1<<" ERROR FIT RC="<<rcfit<<" nstep="<<mFit.GetNStep()<<endl;
if (_prtlevel_>0) mFit.PrintFitErr(rcfit);
}
/* DEL if (pox) {
for(size_t j=0; j<NTRK; j++) {
TVector< complex<double> > expsignal = gxi2.getExpectedSignal(j, v_phase[ii], v_Acx[ii][j]);
sprintf(oname,"simcx_%d_%d",(int)ii+1,(int)j+1);
(*pox) << PPFNameTag(oname)<<expsignal;
}
} */
ofr <<setw(4)<<ii+1<<" "<<setw(5)<<rcfit<<setw(8)<<mFit.GetChi2Red()<<" ";
double phase=mFit.GetParm(0); double err_phase=mFit.GetParmErr(0);
if (phase<0.) phase += 2.*M_PI;
if (phase>2.*M_PI) phase -= 2.*M_PI;
cout <<"Phase= "<<setw(10)<<Angle(phase).ToDegree()<<" +/- "<<setw(10)<<Angle(err_phase).ToDegree()<<" deg."<<endl;
ofr <<setw(8)<<Angle(phase).ToDegree()<<" "<<setw(8)<<Angle(err_phase).ToDegree()<<" ";
v_phase[ii]=phase;
v_err_phase[ii]=err_phase;
for(size_t j=0; j<NTRK; j++) {
double Aac=sqrt(v_A[Anum1[ii]][j] * v_A[Anum2[ii]][j]);
double Ai=(useAac?Aac:v_amp[j]);
double A=mFit.GetParm(1+3*j); double err_A=mFit.GetParmErr(1+3*j);
cout << " Trk["<<j<<"] A= "<<A<<" +/- "<<err_A<<" (A/Ai="<<A/Ai<<")"<<endl;
v_Acx[ii][j]=A;
v_Bcx[ii][j]=complex<double>(0.,0.);
v_err_Acx[ii][j]=err_A;
ofr <<setw(8)<<A<<" "<<setw(8)<<err_A<<" ";
}
ofr << endl;
}
// DEL if (pox) delete pox;
fit_cx_done=true;
acxd_.v_cxbeams=v_cxbeams;
acxd_.v_phase=v_phase;
......@@ -713,6 +693,7 @@ int ACxSetFitter::doCxfit(string outfilenamecx, bool useAac)
return 0;
}
int ACxSetFitter::saveExpectedCx(string outcheckfilename)
{
cout << "ACxSetFitter::saveExpectedCx() saving expected cross-cor (and visi-data) to file "<<outcheckfilename<<endl;
......@@ -745,8 +726,8 @@ int ACxSetFitter::saveExpectedCx(string outcheckfilename)
//------------------------ CxBaselineFitter -------------------------------------
CxBaselineFitter::CxBaselineFitter(vector<AcxDataSet> & v_data, vector<TrackSet> & v_tks)
: v_acxd(v_data), v_trks(v_tks), tot_ntrks(0), fit_done(false), xi2red(-9.e9),
bestfitparam(NULL), err_bestfitparam(NULL)
: v_acxd(v_data), v_trks(v_tks), tot_ntrks(0), fit_done(false), simplex_done(false),
xi2red(-9.e9), bestfitparam(NULL), err_bestfitparam(NULL)
{
if (v_acxd.size() != v_trks.size())
throw ParmError("CxBaselineFitter::CxBaselineFitter(v_data, v_tks) NOT same size v_data,v_tks ");
......@@ -762,6 +743,8 @@ CxBaselineFitter::CxBaselineFitter(vector<AcxDataSet> & v_data, vector<TrackSet>
size_t nparam = 4*(v_acxd[0].getNbAutoCor()-1);
bestfitparam = new double[nparam];
err_bestfitparam = new double[nparam];
initFitParams();
}
CxBaselineFitter::~CxBaselineFitter()
......@@ -770,6 +753,28 @@ CxBaselineFitter::~CxBaselineFitter()
if (err_bestfitparam) delete[] err_bestfitparam;
}
void CxBaselineFitter::initFitParams()
{
size_t NB_ANTENNES=v_acxd[0].getNbAutoCor(); // nombre d'antennes
size_t NB_CXCORS=v_acxd[0].getNbCrossCor();
if (NB_ANTENNES != 4)
throw PError("CxBaselineFitter::initFitParams() NB_ANTENNES != 4 Current version works only for 4 antenna");
v_phases.resize(NB_ANTENNES-1);
v_err_phases.resize(NB_ANTENNES-1);
v_baselineshits.resize(NB_ANTENNES-1);
v_err_baselineshits.resize(NB_ANTENNES-1);
for(size_t i=0; i<(NB_ANTENNES-1); i++) {
v_phases[i]=v_acxd[0].v_phase[i]; v_err_phases[i]=0.;
v_baselineshits[i]=Vector3d(0.,0.,0.);
v_err_baselineshits[i]=Vector3d(0.,0.,0.);
bestfitparam[i]=v_phases[i];
err_bestfitparam[i]=0.;
for(size_t j=0; j<3; j++) {
bestfitparam[3*(i+1)+j]=err_bestfitparam[3*(i+1)+j]=0.;
}
}
}
int CxBaselineFitter::dofit(string outfilename, bool fgfixbaseline)
{
size_t NB_ANTENNES=v_acxd[0].getNbAutoCor(); // nombre d'antennes
......@@ -808,25 +813,18 @@ int CxBaselineFitter::dofit(string outfilename, bool fgfixbaseline)
mFit.SetData(&gdata); // connect data to the fitter , here the data is unused - gxi2 includes its data
mFit.SetMaxStep(1000);
if (NB_ANTENNES != 4)
throw PError("CxBaselineFitter::dofit() NB_ANTENNES != 4 Current version works only for 4 antenna");
v_phases.resize(NB_ANTENNES-1);
v_err_phases.resize(NB_ANTENNES-1);
v_baselineshits.resize(NB_ANTENNES-1);
v_err_baselineshits.resize(NB_ANTENNES-1);
// SetParam(int n,double value, double step,double min=1., double max=-1.);
for(size_t i=0; i<(NB_ANTENNES-1); i++) {
char pname[32];
sprintf(pname,"Phase_%d",(int)(i+2));
mFit.SetParam(i,pname,v_acxd[0].v_phase[i],M_PI/180.,0.,2.5*M_PI);
v_phases[i]=v_acxd[0].v_phase[i]; v_err_phases[i]=0.;
mFit.SetParam(i,pname,v_phases[i],M_PI/180.,0.,2.5*M_PI);
v_err_phases[i]=0.;
sprintf(pname,"BaselineShift_X_%d",(int)(i+2));
mFit.SetParam(3+3*i,pname,0.,0.02,-0.25,0.25);
mFit.SetParam(3+3*i,pname,v_baselineshits[i].X(),0.02,-0.25,0.25);
sprintf(pname,"BaselineShift_Y_%d",(int)(i+2));
mFit.SetParam(4+3*i,pname,0.,0.02,-0.25,0.25);
mFit.SetParam(4+3*i,pname,v_baselineshits[i].Y(),0.02,-0.25,0.25);
sprintf(pname,"BaselineShift_Z_%d",(int)(i+2));
mFit.SetParam(5+3*i,pname,0.,0.02,-0.25,0.25);
mFit.SetParam(5+3*i,pname,v_baselineshits[i].Z(),0.02,-0.25,0.25);
if (fgfixbaseline) {
cout << " ... fitting phases only, fixed baselines "<<endl;
mFit.SetFix(3+3*i); mFit.SetFix(4+3*i); mFit.SetFix(5+3*i);
......@@ -860,6 +858,122 @@ int CxBaselineFitter::dofit(string outfilename, bool fgfixbaseline)
return 0;
}
int CxBaselineFitter::doSimplexMinimize()
{
size_t NB_ANTENNES=v_acxd[0].getNbAutoCor(); // nombre d'antennes
size_t NB_CXCORS=v_acxd[0].getNbCrossCor();
cout << "======================================================================================"<<endl;
cout << "------- CxBaselineFitter::doSimplexMinimize() Performing baseline/phase determination using the 6 cross-cors "<<" TotNbTracks="<<tot_ntrks<<endl;
if (NB_ANTENNES != 4)
throw PError("CxBaselineFitter::doSimplexMinimize() NB_ANTENNES != 4 Current version works only for 4 antenna");
My6CxMinZFunc mzfunc(v_acxd, v_trks);
MinZSimplex simplex(&mzfunc);
// Guess the center and step for constructing the initial simplex
size_t nparam = 4*(NB_ANTENNES-1);
Vector P0(nparam);
Vector step(nparam);
for(size_t i=0; i<(NB_ANTENNES-1); i++) {
P0(i)=v_acxd[0].v_phase[i];
step(i)=M_PI/6.;
for(size_t j=0;j<3;j++) {
P0((i+1)*3+j)=0.;
step((i+1)*3+j)=0.05;
}
}
cout << " Initial Point: "<<P0.Transpose()<<endl;
cout << " Initial Step: "<<step.Transpose()<<endl;
cout << " Initial Xi2= " << mzfunc.Value(P0.Data())<<endl;
simplex.SetInitialPoint(P0);
simplex.SetInitialStep(step);
simplex.SetPrtLevel(_prtlevel_);
Vector oparm(nparam);
int rc = simplex.Minimize(oparm);
if (rc != 0) {
string srt;
int sr = simplex.StopReason(srt);
cout << " Convergence Pb, StopReason= " << sr << " : " << srt << endl;
}
else {
cout << " Converged: NStep= " << simplex.NbIter() << " Best Xi2="<< mzfunc.Value(oparm.Data()) << endl;
simplex_done=true;
for(size_t i=0; i<(NB_ANTENNES-1); i++) {
v_phases[i]=oparm(i);
double xs=oparm(i*3+3);
double ys=oparm(i*3+4);
double zs=oparm(i*3+5);
v_baselineshits[i]=Vector3d(xs,ys,zs);
cout << " ANTENNE["<<i+2<<"] : Phase="<<v_phases[i]<<" BaseLineShift="<<v_baselineshits[i]<<endl;
}
}
return 0;
}
int CxBaselineFitter::doCheck(bool fgminimize)
{
size_t NB_ANTENNES=v_acxd[0].getNbAutoCor(); // nombre d'antennes
size_t NB_CXCORS=v_acxd[0].getNbCrossCor();
cout << "======================================================================================"<<endl;
cout << "------- CxBaselineFitter::doCheck() Performing baseline/phase determination using the 6 cross-cors "<<" TotNbTracks="<<tot_ntrks<<endl;
if (NB_ANTENNES != 4)
throw PError("CxBaselineFitter::doCheck() NB_ANTENNES != 4 Current version works only for 4 antenna");
My6CxMinZFunc mzfunc(v_acxd, v_trks);
mzfunc.SetPrintLevel(_prtlevel_);
// Guess the center and step for constructing the initial simplex
size_t nparam = 4*(NB_ANTENNES-1);
Vector P0(nparam), PC(nparam);
Vector step(nparam);
for(size_t i=0; i<(NB_ANTENNES-1); i++) {
P0(i)=v_acxd[0].v_phase[i];
step(i)=M_PI;
for(size_t j=0;j<3;j++) {
P0((i+1)*3+j)=0.;
step((i+1)*3+j)=0.20;
}
}
cout << " ---- Initial Point: "<<P0.Transpose()<<endl;
cout << " Initial Xi2= " << mzfunc.Value(P0.Data())<<endl;
double pstep=M_PI/20.;
double zstep=0.05;
double bestxi2=9.e19;
Vector oparm(nparam);
size_t cnt=0;
for(int i1=-1; i1<=1; i1++) {
PC(0)=P0(0)+(double)i1*pstep;
for(int i2=-1; i2<=1; i2++) {
PC(1)=P0(1)+(double)i2*pstep;
for(int i3=-1; i3<=1; i3++) {
PC(2)=P0(2)+(double)i3*pstep;
for(int j1=-1; j1<=1; j1++) {
PC(5)=P0(5)+(double)j1*zstep;
for(int j2=-1; j2<=1; j2++) {
PC(8)=P0(8)+(double)j2*zstep;
for(int j3=-1; j3<=0; j3++) {
PC(11)=P0(11)+(double)j3*zstep;
double xi2=mzfunc.Value(PC.Data());
if (xi2<bestxi2) { bestxi2=xi2; oparm=PC; }
cnt++;
}
}
}
}
}
}
cout << "End of Check-Loop Count= " << cnt << " Best Xi2="<< mzfunc.Value(oparm.Data()) << " for :"<<endl;
cout << oparm.Transpose();
return 0;
}
int CxBaselineFitter::saveExpectedCx(string outcheckfilename)
{
......@@ -869,17 +983,19 @@ int CxBaselineFitter::saveExpectedCx(string outcheckfilename)
char oname[48];
My6CxSignalsB cxsigb(v_acxd, v_trks);
cxsigb.SetPrintLevel(_prtlevel_);
for(size_t i=0; i<v_acxd.size(); i++) {
for(size_t j=0; j<v_acxd[i].NbTrk(); j++) {
Vector tmvec = cxsigb.getTimeVec(i,j);
sprintf(oname,"tim_%d_%d",(int)i+1,(int)j+1);
sprintf(oname,"tim_%d_%d",(int)j+1,(int)i+1);
pox << PPFNameTag(oname)<<tmvec;
for(size_t kcx=0; kcx<NB_CXCORS; kcx++) {
TVector< complex<double> > signal = cxsigb.getDataSignal(i,j,kcx);
sprintf(oname,"cx_%d_%d_%d",(int)kcx+1,(int)i+1,(int)j+1);
sprintf(oname,"cx_%d_%d_%d",(int)kcx+1,(int)j+1,(int)i+1);
pox << PPFNameTag(oname)<<signal;
TVector< complex<double> > expsignal = cxsigb.getExpectedSignal(i,j,kcx,bestfitparam);
sprintf(oname,"simcx_%d_%d_%d",(int)kcx+1,(int)i+1,(int)j+1);
sprintf(oname,"simcx_%d_%d_%d",(int)kcx+1,(int)j+1,(int)i+1);
pox << PPFNameTag(oname)<<expsignal;
}
}
......
......@@ -173,8 +173,13 @@ public:
CxBaselineFitter(vector<AcxDataSet> & v_data, vector<TrackSet> & v_tks);
~CxBaselineFitter();
void initFitParams();
int dofit(string outfilename, bool fgfixbaseline=false); // if fgfixbaseline = true, fit phases only
int doCheck(bool fgminimize=false);
int doSimplexMinimize();
int saveExpectedCx(string outcheckfilename); // outfilename PPF (or fits ?) file with the expected signals
vector<AcxDataSet> & v_acxd;
......@@ -182,6 +187,7 @@ public:
size_t tot_ntrks;
bool fit_done;
bool simplex_done;
int rcfit;
double xi2red;
vector<double> v_phases;
......
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