Commit d53df5e6 authored by Pilar Cossio's avatar Pilar Cossio
Browse files

Fourier Space CTF Maximizing Param

parent 70793f6f
......@@ -639,29 +639,43 @@ int bioem::run()
ofstream outputProbFile;
outputProbFile.open ("Output_Probabilities");
outputProbFile <<"************************* HEADER:: NOTATION *******************************************\n";
outputProbFile << " RefMap: MapNumber ; LogProb: natural logarithm of posterior Probability ; Constant: Numerical Const. for adding Probabilities \n";
outputProbFile << " RefMap: MapNumber ; Maximizing Param: MaxLogProb - alpha - beta - gamma - PSF amp - PSF phase - PSF envelope - center x - center y - normalization - offsett \n";
if(param.writeCTF) outputProbFile << " RefMap: MapNumber ; CTFMaxParm: defocus - b-factor (ref. Penzeck 2010)\n";
outputProbFile <<"************************* HEADER:: NOTATION *******************************************\n\n";
for (int iRefMap = 0; iRefMap < RefMap.ntotRefMap; iRefMap ++)
{
// **** Total Probability ***
bioem_Probability_map& pProbMap = pProb.getProbMap(iRefMap);
outputProbFile << "RefMap " << iRefMap << " Probability " << log(pProbMap.Total) + pProbMap.Constoadd + 0.5 * log(M_PI) + (1 - param.param_device.Ntotpi * 0.5)*(log(2 * M_PI) + 1) + log(param.param_device.volu) << " Constant " << pProbMap.Constoadd << "\n";
outputProbFile << "RefMap: " << iRefMap << " LogProb: " << log(pProbMap.Total) + pProbMap.Constoadd + 0.5 * log(M_PI) + (1 - param.param_device.Ntotpi * 0.5)*(log(2 * M_PI) + 1) + log(param.param_device.volu) << " Constant: " << pProbMap.Constoadd << "\n";
outputProbFile << "RefMap " << iRefMap << " Maximizing Param: ";
outputProbFile << "RefMap: " << iRefMap << " Maximizing Param: ";
// *** Param that maximize probability****
outputProbFile << (pProbMap.Constoadd + 0.5 * log(M_PI) + (1 - param.param_device.Ntotpi * 0.5) * (log(2 * M_PI) + 1) + log(param.param_device.volu)) << " ";
outputProbFile << param.angles[pProbMap.max.max_prob_orient].pos[0] << " ";
outputProbFile << param.angles[pProbMap.max.max_prob_orient].pos[1] << " ";
outputProbFile << param.angles[pProbMap.max.max_prob_orient].pos[2] << " ";
outputProbFile << param.CtfParam[pProbMap.max.max_prob_conv].pos[0] << " ";
outputProbFile << param.CtfParam[pProbMap.max.max_prob_conv].pos[1] << " ";
outputProbFile << param.CtfParam[pProbMap.max.max_prob_conv].pos[2] << " ";
outputProbFile << pProbMap.max.max_prob_cent_x << " ";
outputProbFile << pProbMap.max.max_prob_cent_y << " " ;
outputProbFile << pProbMap.max.max_prob_norm << " " ; // param.param_device.Ntotpi / param.param_device.Ntotpi << " " ;
outputProbFile << pProbMap.max.max_prob_mu ;
outputProbFile << param.angles[pProbMap.max.max_prob_orient].pos[0] << " [rad] ";
outputProbFile << param.angles[pProbMap.max.max_prob_orient].pos[1] << " [rad] ";
outputProbFile << param.angles[pProbMap.max.max_prob_orient].pos[2] << " [rad] ";
outputProbFile << param.CtfParam[pProbMap.max.max_prob_conv].pos[0] << " [ ] ";
outputProbFile << param.CtfParam[pProbMap.max.max_prob_conv].pos[1] << " [1./A²] ";
outputProbFile << param.CtfParam[pProbMap.max.max_prob_conv].pos[2] << " [1./A²] ";
outputProbFile << pProbMap.max.max_prob_cent_x << " [pix] ";
outputProbFile << pProbMap.max.max_prob_cent_y << " [pix] " ;
outputProbFile << pProbMap.max.max_prob_norm << " [ ] " ;
outputProbFile << pProbMap.max.max_prob_mu << " [ ] ";
outputProbFile << "\n";
if(param.writeCTF){
myfloat_t denomi;
denomi = param.CtfParam[pProbMap.max.max_prob_conv].pos[1] * param.CtfParam[pProbMap.max.max_prob_conv].pos[1] + param.CtfParam[pProbMap.max.max_prob_conv].pos[2] * param.CtfParam[pProbMap.max.max_prob_conv].pos[2];
outputProbFile << "RefMap: " << iRefMap << " CTFMaxParam: ";
outputProbFile << 2*M_PI*param.CtfParam[pProbMap.max.max_prob_conv].pos[1]/denomi/param.elecwavel << " [micro-m]; ";
outputProbFile << "2*(" << sqrt(4*M_PI*param.CtfParam[pProbMap.max.max_prob_conv].pos[2]/denomi) << ")² [1./A²] \n";
}
// *** For individual files*** //angProbfile.open ("ANG_PROB_"iRefMap);
if(param.param_device.writeAngles)
......
......@@ -42,6 +42,8 @@ public:
int forprintBest(const char* fileinput);
void PrepareFFTs();
bool doaaradius;
bool writeCTF;
myfloat_t elecwavel;
bioem_param_device param_device;
......
......@@ -76,12 +76,13 @@ int bioem_param::readParameters(const char* fileinput,const char* fileangles)
//Default not flipped micrographs
param_device.flipped=false;
param_device.debugterm=false;
param_device.writeCC=false;
writeCTF=false;
elecwavel=0.0;
//Assigning name of angles
//Storing angle file name if existing.
inanglef=std::string(fileangles);
printf("Reading Input Data %s %s\n",inanglef.c_str(),fileangles);
NotUn_angles=0;
ifstream input(fileinput);
......@@ -261,6 +262,18 @@ int bioem_param::readParameters(const char* fileinput,const char* fileangles)
NotUn_angles=int(atoi(token));
cout << "Not uniform total angel: " << NotUn_angles<< "\n";
}
else if (strcmp(token, "WRITE_CTF_PARAM") == 0)//Number of Euler angle tripplets in non uniform Euler angle sampling
{
writeCTF=true;
token = strtok(NULL," ");
elecwavel=atof(token);
cout << "Writing CTF parameters from PSF param that maximize the posterior \n";
if(elecwavel < 200 ){
cout << "Wrong electron wave length " << elecwavel << "\n";
cout << "Has to be in Angstrom (A)\n";
exit(1);}
cout << "Electron wave length in (A) is: " << elecwavel << "\n";
}
}
......
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