28 #include <boost/format.hpp>
30 #include "cctbx/sgtbx/space_group.h"
32 #include "ObjCryst/ObjCryst/PowderPattern.h"
33 #include "ObjCryst/ObjCryst/Molecule.h"
34 #include "ObjCryst/ObjCryst/PowderPatternBackgroundBayesianMinimiser.h"
35 #include "ObjCryst/RefinableObj/Simplex.h"
36 #include "ObjCryst/Quirks/VFNDebug.h"
37 #include "ObjCryst/Quirks/VFNStreamFormat.h"
38 #include "ObjCryst/ObjCryst/CIF.h"
39 #include "ObjCryst/Quirks/Chronometer.h"
41 #include "ObjCryst/wxCryst/wxPowderPattern.h"
55 #define POSSIBLY_UNUSED(expr) (void)(expr)
68 CylinderAbsCorr::~CylinderAbsCorr()
74 const static string mName=
"CylinderAbsCorr";
80 const static string className=
"CylinderAbsCorr";
89 TAU_PROFILE(
"CylinderAbsCorr::CalcCorr()",
"void ()",TAU_DEFAULT);
102 const REAL t0 = 16.0/(3.*M_PI);
103 const REAL t1 = (25.99978-0.01911*pow(s2,0.25))*exp(-0.024551*s2)+ 0.109561*sqrt(s2)-26.04556;
104 const REAL t2 = -0.02489-0.39499*s2+1.219077*pow(s2,1.5)- 1.31268*pow(s2,2)+0.871081*pow(s2,2.5)-0.2327*pow(s2,3);
105 const REAL t3 = 0.003045+0.018167*s2-0.03305*pow(s2,2);
106 const REAL t = -t0*muR-t1*pow(muR,2)-t2*pow(muR,3)-t3*pow(muR,4);
111 const REAL t1 = 1.433902+11.07504*s2-8.77629*s2*s2+ 10.02088*s2*s2*s2-3.36778*s2*s2*s2*s2;
112 const REAL t2 = (0.013869-0.01249*s2)*exp(3.27094*s2)+ (0.337894+13.77317*s2)/pow(1.0+11.53544*s2, 1.555039);
113 const REAL t3 = 1.933433/pow(1.0+23.12967*s2, 1.686715) -0.13576*sqrt(s2)+1.163198;
114 const REAL t4 = 0.044365-0.04259/pow(1.0+0.41051*s2, 148.4202);
115 const REAL t = (t1-t4)/pow(1+t2*(muR-3),t3)+t4;
131 PowderPatternComponent::PowderPatternComponent():
132 mIsScalable(false),mpParentPowderPattern(0)
135 mClockMaster.AddChild(mClockBraggLimits);
138 PowderPatternComponent::PowderPatternComponent(
const PowderPatternComponent &old):
139 mIsScalable(old.mIsScalable),
140 mpParentPowderPattern(old.mpParentPowderPattern)
142 mClockMaster.AddChild(mClockBraggLimits);
143 if(mpParentPowderPattern!=0)
145 mClockMaster.AddChild(mpParentPowderPattern->GetClockPowderPatternPar());
146 mClockMaster.AddChild(mpParentPowderPattern->GetClockPowderPatternXCorr());
147 mClockMaster.AddChild(mpParentPowderPattern->GetClockPowderPatternRadiation());
151 PowderPatternComponent::~PowderPatternComponent()
158 const static string className=
"PowderPatternComponent";
167 std::map<RefinablePar*,CrystVector_REAL>& PowderPatternComponent::GetPowderPattern_FullDeriv(std::set<RefinablePar *> &vPar)
169 this->CalcPowderPattern_FullDeriv(vPar);
170 return mPowderPattern_FullDeriv;
173 std::map<RefinablePar*,CrystVector_REAL>& PowderPatternComponent::GetPowderPatternIntegrated_FullDeriv(std::set<RefinablePar *> &vPar)
176 return mPowderPatternIntegrated_FullDeriv;
200 void PowderPatternComponent::CalcPowderPattern_FullDeriv(std::set<RefinablePar *> &vPar)
202 TAU_PROFILE(
"PowderPatternComponent::CalcPowderPattern_FullDeriv()",
"void ()",TAU_DEFAULT);
203 mPowderPattern_FullDeriv.clear();
204 for(std::set<RefinablePar *>::iterator par=vPar.begin();par!=vPar.end();par++)
211 const REAL step=(*par)->GetDerivStep();
212 (*par)->Mutate(step);
214 (*par)->Mutate(-2*step);
216 (*par)->Mutate(step);
217 mPowderPattern_FullDeriv[*par]/= step*2;
218 if(MaxAbs(mPowderPattern_FullDeriv[*par])==0)
219 mPowderPattern_FullDeriv[*par].resize(0);
225 TAU_PROFILE(
"PowderPatternComponent::CalcPowderPatternIntegrated_FullDeriv()",
"void ()",TAU_DEFAULT);
226 mPowderPatternIntegrated_FullDeriv.clear();
227 for(std::set<RefinablePar *>::iterator par=vPar.begin();par!=vPar.end();par++)
234 const REAL step=(*par)->GetDerivStep();
235 (*par)->Mutate(step);
237 (*par)->Mutate(-2*step);
239 (*par)->Mutate(step);
240 mPowderPatternIntegrated_FullDeriv[*par]/= step*2;
241 if(MaxAbs(mPowderPatternIntegrated_FullDeriv[*par])==0)
242 mPowderPatternIntegrated_FullDeriv[*par].resize(0);
251 PowderPatternBackground::PowderPatternBackground():
252 mBackgroundNbPoint(0),
253 mMaxSinThetaOvLambda(10),mModelVariance(0)
259 PowderPatternBackground::PowderPatternBackground(
const PowderPatternBackground &old):
260 mBackgroundNbPoint(old.mBackgroundNbPoint),
261 mBackgroundInterpPointX(old.mBackgroundInterpPointX),
262 mBackgroundInterpPointIntensity(old.mBackgroundInterpPointIntensity),
263 mMaxSinThetaOvLambda(10),mModelVariance(0)
267 mInterpolationModel.SetChoice(old.mInterpolationModel.GetChoice());
270 PowderPatternBackground::~PowderPatternBackground(){}
273 const static string className=
"PowderPatternBackground";
293 pair<const CrystVector_REAL*,const RefinableObjClock*>
296 VFN_DEBUG_MESSAGE(
"PowderPatternBackground::GetPowderPatternIntegratedCalc()",3)
297 this->CalcPowderPatternIntegrated();
303 VFN_DEBUG_MESSAGE(
"PowderPatternBackground::ImportUserBackground():"<<filename,5)
304 ifstream fin (filename.c_str());
308 Error opening file for input:"+filename);
311 CrystVector_REAL bckgd2Theta(max),bckgd(max);
316 fin >> bckgd2Theta(nbPoints);
317 fin >> bckgd(nbPoints);
319 VFN_DEBUG_MESSAGE(
"Background=" << bckgd(nbPoints)\
320 <<
" at 2theta="<<bckgd2Theta(nbPoints),3)
325 bckgd2Theta.resizeAndPreserve(max);
326 bckgd.resizeAndPreserve(max);
328 }
while(fin.eof() ==
false);
330 bckgd2Theta.resizeAndPreserve(nbPoints);
331 bckgd.resizeAndPreserve(nbPoints);
335 bckgd2Theta*= DEG2RAD;
336 }
else bckgd2Theta*= DEG2RAD;
338 this->SetInterpPoints(bckgd2Theta,bckgd);
339 this->InitRefParList();
343 sprintf(buf,
"Imported %d background points",(
int)nbPoints);
344 (*fpObjCrystInformUser)((string)buf);
347 VFN_DEBUG_MESSAGE(
"PowderPatternBackground::ImportUserBackground():finished",5)
349 void PowderPatternBackground::SetInterpPoints(
const CrystVector_REAL tth,
350 const CrystVector_REAL backgd)
352 VFN_DEBUG_ENTRY(
"PowderPatternBackground::SetInterpPoints():",5)
353 if( (tth.numElements()!=backgd.numElements())
354 ||(tth.numElements()<2))
357 number of points differ or less than 2 points !");
371 this->InitRefParList();
373 VFN_DEBUG_EXIT(
"PowderPatternBackground::SetInterpPoints()",5)
376 const pair<const CrystVector_REAL*,const CrystVector_REAL*> PowderPatternBackground::GetInterpPoints()
const
382 CrystVector_uint & groupIndex,
383 unsigned int &first)
const
386 unsigned int index=0;
387 VFN_DEBUG_MESSAGE(
"PowderPatternBackground::GetGeneGroup()",4)
389 for(
long j=0;j<this->
GetNbPar();j++)
392 if(index==0) index=first++;
398 const bool enableRestraints)
409 pair<const CrystVector_REAL*,const RefinableObjClock*>
412 this->CalcPowderPatternIntegrated();
419 #ifdef USE_BACKGROUND_MAXLIKE_ERROR
432 VFN_DEBUG_ENTRY(
"PowderPatternBackground::OptimizeBayesianBackground()",5);
433 TAU_PROFILE(
"PowderPatternBackground::OptimizeBayesianBackground()",
"void ()",TAU_DEFAULT);
441 cout<<
"Initial Chi^2(BayesianBackground)="<<llk<<endl;
447 sprintf(buf,
"Optimizing Background, Cycle %d, Chi^2(Background)=%f",
449 (*fpObjCrystInformUser)((string)buf);
454 cout<<ct<<
", Chi^2(BayesianBackground)="<<llk<<endl;
459 sprintf(buf,
"Done Optimizing Bayesian Background, Chi^2(Background)=%f",(
float)llk);
460 (*fpObjCrystInformUser)((string)buf);
467 try {lsq.
Refine(10,
true,
false);}
471 VFN_DEBUG_EXIT(
"PowderPatternBackground::OptimizeBayesianBackground()",5);
494 TAU_PROFILE(
"PowderPatternBackground::CalcPowderPattern()",
"void ()",TAU_DEFAULT);
495 VFN_DEBUG_MESSAGE(
"PowderPatternBackground::CalcPowderPattern()",3);
502 case POWDER_BACKGROUND_LINEAR:
504 VFN_DEBUG_MESSAGE(
"PowderPatternBackground::CalcPowderPattern()..Linear",2)
512 VFN_DEBUG_MESSAGE(
"PowderPatternBackground::CalcPowderPattern()"<<nb,2)
520 for(
unsigned long i=0;i<nb;i++)
533 *b = (b1*(p2-i)+b2*(i-p1))/(p2-p1) ;
538 case POWDER_BACKGROUND_CUBIC_SPLINE:
549 VFN_DEBUG_MESSAGE(
"PowderPatternBackground::CalcPowderPattern()",3);
550 #ifdef USE_BACKGROUND_MAXLIKE_ERROR
556 for(
long i=0;i<nb;i++) {*p++ = var;var +=step;}
561 VFN_DEBUG_MESSAGE(
"PowderPatternBackground::CalcPowderPattern():End",3);
564 void PowderPatternBackground::CalcPowderPattern_FullDeriv(std::set<RefinablePar*> &vPar)
566 TAU_PROFILE(
"PowderPatternBackground::CalcPowderPattern_FullDeriv()",
"void ()",TAU_DEFAULT);
568 mPowderPattern_FullDeriv.clear();
570 for(std::set<RefinablePar*>::iterator par=vPar.begin();par!=vPar.end();++par)
577 const REAL step=(*par)->GetDerivStep();
578 (*par)->Mutate(step);
580 (*par)->Mutate(-2*step);
582 (*par)->Mutate(step);
583 mPowderPattern_FullDeriv[*par]/= step*2;
584 if(MaxAbs(mPowderPattern_FullDeriv[*par])==0)
585 mPowderPattern_FullDeriv[*par].resize(0);
590 void PowderPatternBackground::CalcPowderPatternIntegrated()
const
599 VFN_DEBUG_ENTRY(
"PowderPatternBackground::CalcPowderPatternIntegrated()",3)
600 TAU_PROFILE("PowderPatternBackground::CalcPowderPatternIntegrated()","
void ()",TAU_DEFAULT);
604 const
long numInterval=pMin->numElements();
607 for(
int j=0;j<numInterval;j++)
609 const long max=(*pMax)(j);
612 for(
int k=(*pMin)(j);k<=max;k++) *p2 += *p1++;
615 #ifdef USE_BACKGROUND_MAXLIKE_ERROR
618 for(
int j=0;j<numInterval;j++)
620 const long max=(*pMax)(j);
623 for(
int k=(*pMin)(j);k<=max;k++) *p2 += *p1++;
629 VFN_DEBUG_EXIT(
"PowderPatternBackground::CalcPowderPatternIntegrated()",3)
634 TAU_PROFILE(
"PowderPatternBackground::CalcPowderPatternIntegrated_FullDeriv()",
"void ()",TAU_DEFAULT);
637 mPowderPatternIntegrated_FullDeriv.clear();
639 for(std::set<RefinablePar*>::iterator par=vPar.begin();par!=vPar.end();++par)
646 const REAL step=(*par)->GetDerivStep();
647 (*par)->Mutate(step);
649 (*par)->Mutate(-2*step);
651 (*par)->Mutate(step);
652 mPowderPatternIntegrated_FullDeriv[*par]/= step*2;
653 if(MaxAbs(mPowderPatternIntegrated_FullDeriv[*par])==0)
654 mPowderPatternIntegrated_FullDeriv[*par].resize(0);
658 std::map<RefinablePar*, CrystVector_REAL> newDeriv=mPowderPatternIntegrated_FullDeriv;
660 std::vector<const CrystVector_REAL*> v;
662 for(std::map<RefinablePar*, CrystVector_REAL>::reverse_iterator pos=mPowderPatternIntegrated_FullDeriv.rbegin();pos!=mPowderPatternIntegrated_FullDeriv.rend();++pos)
664 cout<<pos->first->GetName()<<
":"<<pos->second.size()<<
","<<newDeriv[pos->first].size()<<endl;
665 if(pos->second.size()==0)
continue;
666 v.push_back(&(pos->second));
667 v.push_back(&(newDeriv[pos->first]));
670 if(v.size()>0) cout<<
"PowderPatternBackground::CalcPowderPatternIntegrated_FullDeriv():"<<endl<<FormatVertVector<REAL>(v,12,1,20)<<endl;
691 void PowderPatternBackground::InitRefParList()
696 string str=
"Background_Point_";
703 false,
true,
true,
false,1.);
704 tmp.SetGlobalOptimStep(10.);
706 tmp.SetDerivStep(1e-3);
709 #ifdef USE_BACKGROUND_MAXLIKE_ERROR
713 true,
true,
true,
false,1.);
715 tmp.SetDerivStep(1e-3);
724 void PowderPatternBackground::InitOptions()
726 VFN_DEBUG_MESSAGE(
"PowderPatternBackground::InitOptions()",5)
727 static
string InterpolationModelName;
728 static
string InterpolationModelChoices[2];
730 static
bool needInitNames=true;
731 if(true==needInitNames)
733 InterpolationModelName=
"Interpolation Model";
734 InterpolationModelChoices[0]=
"Linear";
735 InterpolationModelChoices[1]=
"Spline";
746 void PowderPatternBackground::InitSpline()
const
750 &&(
mClockSpline>this->GetParentPowderPattern().GetRadiation().GetClockWavelength()))
return;
781 WXCrystObjBasic* PowderPatternBackground::WXCreate(wxWindow* parent)
784 mpWXCrystObj=
new WXPowderPatternBackground(parent,
this);
793 PowderPatternDiffraction::PowderPatternDiffraction():
794 mpReflectionProfile(0),
795 mCorrLorentz(*this),mCorrPolar(*this),mCorrSlitAperture(*this),
796 mCorrTextureMarchDollase(*this),mCorrTextureEllipsoid(*this),mCorrTOF(*this),mCorrCylAbs(*this),mExtractionMode(false),
797 mpLeBailData(0),mFrozenLatticePar(6),mFreezeLatticePar(false),mFrozenBMatrix(3,3),mGenHKLBMatrix(3,3)
799 VFN_DEBUG_MESSAGE(
"PowderPatternDiffraction::PowderPatternDiffraction()",10)
802 this->SetProfile(new ReflectionProfilePseudoVoigt);
803 this->SetIsIgnoringImagScattFact(true);
809 for(
unsigned int i=0;i<3;++i) mFrozenLatticePar(i)=5;
810 for(
unsigned int i=3;i<6;++i) mFrozenLatticePar(i)=M_PI/2;
814 PowderPatternDiffraction::PowderPatternDiffraction(const PowderPatternDiffraction &old):
815 mpReflectionProfile(0),
816 mCorrLorentz(*this),mCorrPolar(*this),mCorrSlitAperture(*this),
817 mCorrTextureMarchDollase(*this),mCorrTextureEllipsoid(*this),mCorrTOF(*this),mCorrCylAbs(*this),mExtractionMode(false),
818 mpLeBailData(0),mFrozenLatticePar(6),mFreezeLatticePar(old.FreezeLatticePar()),mFrozenBMatrix(3,3),mGenHKLBMatrix(3,3)
822 this->SetIsIgnoringImagScattFact(
true);
823 this->SetProfile(old.mpReflectionProfile->CreateCopy());
825 if(old.mpLeBailData!=0)
827 mpLeBailData=
new DiffractionDataSingleCrystal(
false);
828 *mpLeBailData = *(old.mpLeBailData);
834 for(
unsigned int i=0;i<6;++i) mFrozenLatticePar(i)=old.GetFrozenLatticePar(i);
838 PowderPatternDiffraction::~PowderPatternDiffraction()
840 if(mpReflectionProfile!=0)
843 delete mpReflectionProfile;
848 const static string className=
"PowderPatternDiffraction";
875 pair<const CrystVector_REAL*,const RefinableObjClock*>
878 this->CalcPowderPatternIntegrated();
883 const REAL w,
const REAL u,
const REAL v,
884 const REAL eta0,
const REAL eta1)
886 VFN_DEBUG_MESSAGE(
"PowderPatternDiffraction::SetReflectionProfilePar()",5)
894 if(p==mpReflectionProfile)
return;
895 if(mpReflectionProfile!=0)
898 delete mpReflectionProfile;
900 mpReflectionProfile= p;
907 return *mpReflectionProfile;
912 return *mpReflectionProfile;
916 void PowderPatternDiffraction::GenHKLFullSpace(
917 const REAL maxTheta,
const bool useMultiplicity)
const
923 void PowderPatternDiffraction::GenHKLFullSpace()
const
925 VFN_DEBUG_ENTRY(
"PowderPatternDiffraction::GenHKLFullSpace():",5)
927 if(this->GetRadiation().GetWavelengthType()==WAVELENGTH_TOF)
933 if((mExtractionMode) && (mFhklObsSq.numElements()!=this->GetNbRefl()))
935 VFN_DEBUG_ENTRY(
"PowderPatternDiffraction::GenHKLFullSpace(): need to resize observed intensities",7)
936 const
int n0 = mFhklObsSq.numElements();
937 mFhklObsSq.resizeAndPreserve(this->GetNbRefl());
938 for(
int i=n0;i<this->GetNbRefl();i++) mFhklObsSq(i)=100;
941 mGenHKLBMatrix = this->GetBMatrix();
945 const
bool enableRestraints)
947 if(mUseFastLessPreciseFunc!=allowApproximations)
949 mClockProfileCalc.Reset();
950 mClockGeomStructFact.Reset();
951 mClockStructFactor.Reset();
954 mUseFastLessPreciseFunc=allowApproximations;
955 this->GetNbReflBelowMaxSinThetaOvLambda();
962 if(mUseFastLessPreciseFunc==
true)
964 mClockProfileCalc.Reset();
965 mClockGeomStructFact.Reset();
966 mClockStructFactor.Reset();
969 mUseFastLessPreciseFunc=
false;
970 this->GetNbReflBelowMaxSinThetaOvLambda();
977 if(mUseFastLessPreciseFunc!=allow)
979 mClockProfileCalc.Reset();
980 mClockGeomStructFact.Reset();
981 mClockStructFactor.Reset();
984 mUseFastLessPreciseFunc=allow;
985 this->GetNbReflBelowMaxSinThetaOvLambda();
990 CrystVector_uint & groupIndex,
991 unsigned int &first)
const
994 unsigned int index=0;
995 VFN_DEBUG_MESSAGE(
"PowderPatternDiffraction::GetGeneGroup()",4)
997 for(
long j=0;j<this->
GetNbPar();j++)
1002 if(index==0) index=first++;
1003 groupIndex(i)=index;
1015 pair<const CrystVector_REAL*,const RefinableObjClock*>
1018 this->CalcPowderPatternIntegrated();
1030 bool reprep=(mpCrystal!=0);
1033 if(mpReflectionProfile!=0)
1034 if(mpReflectionProfile->GetClassName()==
"ReflectionProfileDoubleExponentialPseudoVoigt")
1049 VFN_DEBUG_ENTRY(
"PowderPatternDiffraction::SetExtractionMode(),ExtractionMode="<<mExtractionMode<<
", nbrefl="<<this->GetNbRefl(),7)
1050 mExtractionMode=extract;
1051 bool needInit=
false;
1054 this->FreezeLatticePar(
false);
1056 mFhklObsSq.resizeAndPreserve(this->GetNbRefl());
1058 if(extract && (!init) && (mpLeBailData!=0))
1061 const long nbrefl=this->GetNbReflBelowMaxSinThetaOvLambda();
1062 if(nbrefl==mpLeBailData->GetNbRefl())
1064 for(
int i=0;i<nbrefl;++i)
1066 if( (mpLeBailData->GetH()(i)==this->GetH()(i))
1067 &&(mpLeBailData->GetK()(i)==this->GetK()(i))
1068 &&(mpLeBailData->GetL()(i)==this->GetL()(i)))
1070 mFhklObsSq(i)=mpLeBailData->GetFhklObsSq()(i);
1075 VFN_DEBUG_MESSAGE(
"PowderPatternDiffraction::SetExtractionMode():: Forcing initialize, cannot re-use: hkl list differs",10);
1084 VFN_DEBUG_MESSAGE(
"PowderPatternDiffraction::SetExtractionMode():: Forcing initialize, cannot re-use: different number of reflections",10);
1087 if((extract && init) || needInit)
1089 VFN_DEBUG_MESSAGE(
"PowderPatternDiffraction::SetExtractionMode():: Initializing intensities to 100",10);
1092 if((mExtractionMode==
false)&&(mFhklObsSq.numElements()>0))
1094 VFN_DEBUG_ENTRY(
"PowderPatternDiffraction::SetExtractionMode(),LEAVING Le Bail Mode",7)
1097 mpLeBailData->SetWavelength(this->GetRadiation().GetWavelength()(0));
1098 mpLeBailData->SetRadiationType(this->GetRadiation().GetRadiationType());
1101 mpLeBailData->SetName(
string(buf)+this->GetCrystal().
GetName());
1103 const unsigned long nbrefl=this->GetNbReflBelowMaxSinThetaOvLambda();
1104 CrystVector_REAL iobs(nbrefl),sigma(nbrefl);
1105 CrystVector_long h(nbrefl),k(nbrefl),l(nbrefl);
1107 for(
unsigned long i=0;i<nbrefl;++i)
1112 iobs(i)=mFhklObsSq(i);
1114 mpLeBailData->SetHklIobs(h,k,l,iobs,sigma);
1116 mFhklObsSq.resize(0);
1119 mClockFhklObsSq.Click();
1120 VFN_DEBUG_EXIT(
"PowderPatternDiffraction::SetExtractionMode(),ExtractionMode="<<mExtractionMode<<
", nbrefl="<<this->GetNbRefl(),7)
1127 VFN_DEBUG_ENTRY(
"PowderPatternDiffraction::ExtractLeBail()",7)
1128 TAU_PROFILE(
"PowderPatternDiffraction::ExtractLeBail()",
"void (int)",TAU_DEFAULT);
1130 if(mExtractionMode==
false) this->SetExtractionMode(
true,
true);
1131 if(mFhklObsSq.numElements()!=this->GetNbRefl())
1133 VFN_DEBUG_ENTRY(
"PowderPatternDiffraction::ExtractLeBail() mFhklObsSq.size() != NbRefl !!!!!!",7)
1134 mFhklObsSq.resize(this->GetNbRefl());
1138 CrystVector_REAL obs,iextract,calc;
1139 iextract=mFhklObsSq;
1141 mClockFhklObsSq.Click();
1145 mFhklObsSq=iextract;
1146 mClockFhklObsSq.Click();
1152 for(;nbcycle>0;nbcycle--)
1156 for(
unsigned int k0=0;k0<nbrefl;++k0)
1158 if(mvReflProfile[k0].profile.numElements()==0)
continue;
1161 long last=mvReflProfile[k0].last,first;
1163 if(mvReflProfile[k0].first<0)first=0;
1164 else first=(mvReflProfile[k0].first);
1165 const REAL *p1=mvReflProfile[k0].profile.data()+(first-mvReflProfile[k0].first);
1166 const REAL *p2=calc.data()+first;
1167 const REAL *pobs=obs.data()+first;
1168 for(
long i=first;i<=last;++i)
1170 const REAL s2=*p2++;
1171 const REAL tmp=*pobs++ * *p1++;
1180 if((s1>1e-8)&&(!
ISNAN_OR_INF(s1))) iextract(k0)=s1*mFhklObsSq(k0);
1181 else iextract(k0)=1e-8;
1184 mFhklObsSq=iextract;
1185 if(this->GetCrystal().GetScatteringComponentList().GetNbComponent()>0)
1187 const REAL* p1=this->GetFhklCalcSq() .data();
1188 const REAL* p2=mFhklObsSq.data();
1190 for(
long i=nbrefl;i>0;i--)
1192 tmp1 += (*p1) * (*p2++);
1193 tmp2 += (*p1) * (*p1);
1197 mFhklObsSq*=tmp2/tmp1;
1199 mClockFhklObsSq.Click();
1201 mClockIhklCalc.Reset();
1206 VFN_DEBUG_MESSAGE(
"PowderPatternDiffraction::ExtractLeBail(): creating single crystal extracted data",7)
1207 CrystVector_REAL iobs(nbrefl),sigma(nbrefl);
1208 CrystVector_long h(nbrefl),k(nbrefl),l(nbrefl);
1210 for(
unsigned long i=0;i<nbrefl;++i)
1215 iobs(i)=mFhklObsSq(i);
1217 mpLeBailData->SetHklIobs(h,k,l,iobs,sigma);
1219 VFN_DEBUG_EXIT(
"PowderPatternDiffraction::ExtractLeBail()mFhklObsSq.size()=="<<mFhklObsSq.numElements(),7)
1224 VFN_DEBUG_MESSAGE(
"PowderPatternDiffraction::GetNbReflBelowMaxSinThetaOvLambda(): "<<mNbReflUsed<<
"/"<<mNbRefl<<
" [max sin(theta)/lambda="<<
mMaxSinThetaOvLambda<<
"]",4)
1225 this->CalcPowderReflProfile();
1227 if((mNbReflUsed>0)&&(mNbReflUsed<mNbRefl))
1229 if( (mvReflProfile[mNbReflUsed ].first>nbpoint)
1230 &&(mvReflProfile[mNbReflUsed-1].first<=nbpoint))
return mNbReflUsed;
1233 if((mNbReflUsed==mNbRefl) && (mvReflProfile[mNbReflUsed-1].profile.numElements()>0))
1234 if(mvReflProfile[mNbReflUsed-1].first<=nbpoint)
return mNbReflUsed;
1238 for(i=0;i<mNbRefl;i++)
1240 if(mvReflProfile[i].first>nbpoint)
break;
1245 mClockNbReflUsed.Click();
1247 <<
" nb refl="<<mNbReflUsed,4)
1254 const REAL old=mFrozenLatticePar(i);
1255 cout<<
"PowderPatternDiffraction::SetFrozenLatticePar("<<i<<
":"<<v<<
")"<<endl;
1257 mFrozenLatticePar(i)=v;
1258 this->CalcFrozenBMatrix();
1265 VFN_DEBUG_MESSAGE(
"PowderPatternDiffraction::FreezeLatticePar("<<use<<
")", 10)
1266 if(use==mFreezeLatticePar)
return;
1267 mFreezeLatticePar=use;
1268 mFrozenLatticePar=this->GetCrystal().GetLatticePar();
1269 if(use) this->CalcFrozenBMatrix();
1270 mClockTheta.Reset();
1279 unsigned int irefl=0;
1280 unsigned int ilast=0;
1284 if(irefl>=this->GetNbReflBelowMaxSinThetaOvLambda())
break;
1286 while(irefl<this->GetNbReflBelowMaxSinThetaOvLambda())
1288 const REAL stol = mSinThetaLambda(irefl);
1289 while(mSinThetaLambda(irefl)==stol)
1293 if(irefl>=this->GetNbReflBelowMaxSinThetaOvLambda())
break;
1296 if(nbnew>1) nb += nbnew-1;
1306 if(mpLeBailData==NULL)
return false;
1307 return mpLeBailData->GetFhklObsSq().size() > 0;
1312 if(mpLeBailData==NULL)
1313 throw ObjCrystException(
"PowderPatternDiffraction::GetFhklObsSq(): no extracted intensities available");
1314 return mpLeBailData->GetFhklObsSq();
1319 this->GetNbReflBelowMaxSinThetaOvLambda();
1321 TAU_PROFILE(
"PowderPatternDiffraction::CalcPowderPattern()-Apply profiles",
"void (bool)",TAU_DEFAULT);
1323 VFN_DEBUG_ENTRY(
"PowderPatternDiffraction::CalcPowderPattern():",3)
1325 if(this->GetCrystal().GetSpaceGroup().GetClockSpaceGroup()>mClockHKL)
1327 VFN_DEBUG_MESSAGE(
"PowderPatternDiffraction::CalcPowderPattern():"
1328 "spacegroup has changed, re-generating HKL's",5)
1329 this->GenHKLFullSpace();
1331 else if((!this->
IsBeingRefined()) && (this->GetCrystal().GetClockLatticePar()>mClockHKL))
1340 const REAL r = 0.005;
1341 if(abs(this->GetBMatrix()(0)-mGenHKLBMatrix(0))>(r*this->GetBMatrix()(0))) needgen=
true;
1342 else if(abs(this->GetBMatrix()(4)-mGenHKLBMatrix(4))>(r*this->GetBMatrix()(4))) needgen=
true;
1343 else if(abs(this->GetBMatrix()(1)-mGenHKLBMatrix(1))>(r*this->GetBMatrix()(4))) needgen=
true;
1344 else if(abs(this->GetBMatrix()(8)-mGenHKLBMatrix(8))>(r*this->GetBMatrix()(8))) needgen=
true;
1345 else if(abs(this->GetBMatrix()(2)-mGenHKLBMatrix(2))>(r*this->GetBMatrix()(8))) needgen=
true;
1346 else if(abs(this->GetBMatrix()(5)-mGenHKLBMatrix(5))>(r*this->GetBMatrix()(8))) needgen=
true;
1349 VFN_DEBUG_MESSAGE(
"PowderPatternDiffraction::CalcPowderPattern():"
1350 "lattice parameters have changed by more than 0.5%, re-generating HKL's",5)
1351 this->GenHKLFullSpace();
1356 this->CalcPowderReflProfile();
1363 const long nbRefl=this->GetNbRefl();
1364 VFN_DEBUG_MESSAGE(
"PowderPatternDiffraction::CalcPowderPattern\
1365 Applying profiles for "<<nbRefl<<
" reflections",2)
1370 const bool useML= (mIhklCalcVariance.numElements() != 0);
1377 VFN_DEBUG_MESSAGE(
"PowderPatternDiffraction::CalcPowderPattern() Has variance:"<<useML,2)
1379 for(
long i=0;i<mNbRefl;i += step)
1381 if(mvReflProfile[i].profile.numElements()==0)
1384 if(i>=mNbReflUsed)
break;
1388 VFN_DEBUG_MESSAGE(
"PowderPatternDiffraction::CalcPowderPattern()#"<<i,2)
1395 intensity += mIhklCalc(i + step);
1396 if(useML) var += mIhklCalcVariance(i + step);
1398 if(mpReflectionProfile->IsAnisotropic())
break;
1399 if( (i+step) >= nbRefl)
break;
1400 if(mSinThetaLambda(i+step) > (mSinThetaLambda(i)+1e-5) )
break;
1402 VFN_DEBUG_MESSAGE(
"Apply profile(Monochromatic)Refl("<<i<<
")"\
1403 <<mIntH(i)<<
" "<<mIntK(i)<<
" "<<mIntL(i)<<
" "\
1404 <<
" I="<<intensity<<
" stol="<<mSinThetaLambda(i)\
1405 <<
",pixel #"<<mvReflProfile[i].first<<
"->"<<mvReflProfile[i].last,2)
1407 const long first=mvReflProfile[i].first,last=mvReflProfile[i].last;
1408 const REAL *p2 = mvReflProfile[i].profile.data();
1410 for(
long j=first;j<=last;j++) *p3++ += *p2++ * intensity;
1413 const REAL *p2 = mvReflProfile[i].profile.data();
1415 for(
long j=first;j<=last;j++) *p3++ += *p2++ * var;
1424 FAST option not yet implemented !");
1427 VFN_DEBUG_EXIT(
"PowderPatternDiffraction::CalcPowderPattern: End.",3)
1430 void PowderPatternDiffraction::CalcPowderPattern_FullDeriv(std::set<RefinablePar*> &vPar)
1432 TAU_PROFILE(
"PowderPatternDiffraction::CalcPowderPattern_FullDeriv()",
"void ()",TAU_DEFAULT);
1435 bool notYetDerivProfiles=
true;
1436 mIhkl_FullDeriv.clear();
1437 mvReflProfile_FullDeriv.clear();
1438 mPowderPattern_FullDeriv.clear();
1439 for(std::set<RefinablePar*>::iterator par=vPar.begin();par!=vPar.end();++par)
1441 if(*par==0)
continue;
1442 if((*par)->IsFixed())
continue;
1443 if((*par)->IsUsed()==
false)
continue;
1444 if( (*par)->GetType()->IsDescendantFromOrSameAs(gpRefParTypeScatt)
1445 ||(*par)->GetType()->IsDescendantFromOrSameAs(gpRefParTypeScattPow))
1447 this->CalcIhkl_FullDeriv(vPar);
1449 if(notYetDerivProfiles)
1451 if( (*par)->GetType()->IsDescendantFromOrSameAs(gpRefParTypeRadiation)
1452 ||(*par)->GetType()->IsDescendantFromOrSameAs(gpRefParTypeUnitCell)
1456 this->CalcPowderReflProfile_FullDeriv(vPar);
1457 notYetDerivProfiles=
false;
1463 for(std::set<RefinablePar*>::iterator par=vPar.begin();par!=vPar.end();++par)
1468 if((*par)->IsFixed())
continue;
1469 if((*par)->IsUsed()==
false)
continue;
1470 if(mIhkl_FullDeriv[*par].size()!=0)
1472 const long nbRefl=this->GetNbRefl();
1475 mPowderPattern_FullDeriv[*par].resize(specNbPoints);
1476 mPowderPattern_FullDeriv[*par]=0;
1478 for(
long i=0;i<mNbReflUsed;i += step)
1480 if(mvReflProfile[i].profile.numElements()==0)
1483 if(i>=mNbReflUsed)
break;
1492 intensity += mIhkl_FullDeriv[*par](i + step);
1494 if(mpReflectionProfile->IsAnisotropic())
break;
1495 if( (i+step) >= nbRefl)
break;
1496 if(mSinThetaLambda(i+step) > (mSinThetaLambda(i)+1e-5) )
break;
1499 const long first=mvReflProfile[i].first,last=mvReflProfile[i].last;
1500 const REAL *p2 = mvReflProfile[i].profile.data();
1501 REAL *p3 = mPowderPattern_FullDeriv[*par].data()+first;
1502 for(
long j=first;j<=last;j++) *p3++ += *p2++ * intensity;
1506 if(mvReflProfile_FullDeriv[*par].size()!=0)
1508 const long nbRefl=this->GetNbRefl();
1511 mPowderPattern_FullDeriv[*par].resize(specNbPoints);
1512 mPowderPattern_FullDeriv[*par]=0;
1513 cout<<__FILE__<<
":"<<__LINE__<<
":PowderPatternDiffraction::CalcPowderPattern_FullDeriv():par="<<(*par)->GetName()<<endl;
1514 for(
long i=0;i<mNbReflUsed;i += step)
1516 if(mvReflProfile[i].profile.numElements()==0)
1519 if(i>=mNbReflUsed)
break;
1528 intensity += mIhklCalc(i + step);
1530 if(mpReflectionProfile->IsAnisotropic())
break;
1531 if( (i+step) >= nbRefl)
break;
1532 if(mSinThetaLambda(i+step) > (mSinThetaLambda(i)+1e-5) )
break;
1534 if(mvReflProfile_FullDeriv[*par][i].size()>0)
1536 const long first=mvReflProfile[i].first,last=mvReflProfile[i].last;
1537 const REAL *p2 = mvReflProfile_FullDeriv[*par][i].data();
1538 REAL *p3 = mPowderPattern_FullDeriv[*par].data()+first;
1539 for(
long j=first;j<=last;j++) *p3++ += *p2++ * intensity;
1546 std::map<RefinablePar*, CrystVector_REAL> newDeriv=mPowderPattern_FullDeriv;
1547 this->PowderPatternComponent::CalcPowderPattern_FullDeriv(vPar);
1548 std::vector<const CrystVector_REAL*> v;
1550 for(std::map<RefinablePar*, CrystVector_REAL>::reverse_iterator pos=mPowderPattern_FullDeriv.rbegin();pos!=mPowderPattern_FullDeriv.rend();++pos)
1552 if(pos->first==0)
continue;
1553 if(pos->second.size()==0)
continue;
1554 v.push_back(&(newDeriv[pos->first]));
1555 v.push_back(&(pos->second));
1556 cout<<pos->first->GetName()<<
":"<<pos->second.size()<<
","<<newDeriv[pos->first].size()<<endl;
1559 cout<<FormatVertVector<REAL>(v,16,4,1000)<<endl;
1564 void PowderPatternDiffraction::CalcPowderPatternIntegrated()
const
1566 this->GetNbReflBelowMaxSinThetaOvLambda();
1568 TAU_PROFILE(
"PowderPatternDiffraction::CalcPowderPatternIntegrated()",
"void (bool)",TAU_DEFAULT);
1569 TAU_PROFILE_TIMER(timer1,
"PowderPatternDiffraction::CalcPowderPatternIntegrated()1",
"", TAU_FIELD);
1570 TAU_PROFILE_TIMER(timer2,
"PowderPatternDiffraction::CalcPowderPatternIntegrated()2",
"", TAU_FIELD);
1573 TAU_PROFILE_START(timer1);
1574 this->PrepareIntegratedProfile();
1575 TAU_PROFILE_STOP(timer1);
1581 VFN_DEBUG_ENTRY(
"PowderPatternDiffraction::CalcPowderPatternIntegrated()",3)
1582 const
long nbRefl=this->GetNbRefl();
1587 const
bool useML= (mIhklCalcVariance.numElements() != 0);
1594 const REAL * RESTRICT psith=mSinThetaLambda.data();
1595 const REAL * RESTRICT pI=mIhklCalc.data();
1596 const REAL * RESTRICT pIvar=mIhklCalcVariance.data();
1597 vector< pair<unsigned long, CrystVector_REAL> >::const_iterator pos;
1598 pos=mIntegratedProfileFactor.begin();
1599 TAU_PROFILE_START(timer2);
1600 for(
long i=0;i<mNbReflUsed;)
1602 VFN_DEBUG_MESSAGE(
"PowderPatternDiffraction::CalcPowderPatternIntegrated():"<<i,2)
1605 const REAL thmax=*psith+1e-5;
1611 if(useML) var += *pIvar++;
1612 if( ++i >= nbRefl)
break;
1613 if( *(++psith) > thmax )
break;
1614 if(mpReflectionProfile->IsAnisotropic())
break;
1617 VFN_DEBUG_MESSAGE(
"PowderPatternDiffraction::CalcPowderPatternIntegrated():"<<i,2)
1619 const REAL * RESTRICT pFact=pos->second.data();
1620 const
unsigned long nb=pos->second.numElements();
1622 for(
unsigned long j=nb;j>0;j--)
1625 *pData++ += intensity * *pFact++ ;
1631 const REAL * RESTRICT pFact=pos->second.data();
1633 for(
unsigned long j=nb;j>0;j--) *pVar++ += var * *pFact++ ;
1637 TAU_PROFILE_STOP(timer2);
1639 if(gVFNDebugMessageLevel<3)
1642 CrystVector_REAL integr(nb),min(nb),max(nb),diff(nb),index(nb);
1644 for(
long i=0;i<nb;i++)
1651 j<=mpParentPowderPattern->GetIntegratedProfileMax()(i);j++)
1657 cout <<
"Integrated intensities, Component"<<endl
1663 VFN_DEBUG_EXIT(
"PowderPatternDiffraction::CalcPowderPatternIntegrated",3)
1668 TAU_PROFILE(
"PowderPatternDiffraction::CalcPowderPatternIntegrated_FullDeriv()",
"void ()",TAU_DEFAULT);
1673 this->CalcPowderPatternIntegrated();
1674 this->CalcIhkl_FullDeriv(vPar);
1675 const long nbRefl=this->GetNbRefl();
1677 #ifdef PowderPatternDiffraction_CalcPowderPatternIntegrated_FullDerivDEBUG
1679 std::map<RefinablePar*, CrystVector_REAL> oldDeriv=mPowderPatternIntegrated_FullDeriv;
1683 mPowderPatternIntegrated_FullDeriv.clear();
1684 for(std::set<RefinablePar*>::iterator par=vPar.begin();par!=vPar.end();++par)
1689 if(mIhkl_FullDeriv[*par].size()==0)
continue;
1690 if(mPowderPatternIntegrated_FullDeriv[*par].size()==0)
1692 mPowderPatternIntegrated_FullDeriv[*par].resize(nbprof);
1693 mPowderPatternIntegrated_FullDeriv[*par]=0;
1695 const REAL * RESTRICT psith=mSinThetaLambda.data();
1696 const REAL * RESTRICT pI=mIhkl_FullDeriv[*par].data();
1699 vector< pair<unsigned long, CrystVector_REAL> >::const_iterator pos=mIntegratedProfileFactor.begin();
1701 for(
long i=0;i<mNbReflUsed;)
1704 const REAL thmax=*psith+1e-5;
1708 if( ++i >= nbRefl)
break;
1709 if( *(++psith) > thmax )
break;
1710 if(mpReflectionProfile->IsAnisotropic())
break;
1713 REAL * RESTRICT pData=mPowderPatternIntegrated_FullDeriv[*par].data()+pos->first;
1714 const REAL * RESTRICT pFact=pos->second.data();
1715 const unsigned long nb=pos->second.numElements();
1716 #ifdef PowderPatternDiffraction_CalcPowderPatternIntegrated_FullDerivDEBUG
1717 if((i<5)&&(ctpar<8)) cout<<__FILE__<<
":"<<__LINE__<<
":"<<(*par)->GetName()<<
"i="<<setw(16)<<i<<
":I="<<setw(16)<<intensity<<endl;
1719 for(
unsigned long j=nb;j>0;j--)
1721 #ifdef PowderPatternDiffraction_CalcPowderPatternIntegrated_FullDerivDEBUG
1722 *pData += intensity * *pFact ;
1723 if((i<5)&&((*par)->GetName()==
"Cimetidine_C11_x")&&(pos->first==0)&&(nb==j)) cout<<nb-j<<
" SUM1"<<setw(16)<<*pData<<
", dI="<<setw(16)<<intensity<<
", prof="<<setw(16)<<*pFact<<endl;
1726 *pData++ += intensity * *pFact++ ;
1734 #ifdef PowderPatternDiffraction_CalcPowderPatternIntegrated_FullDerivDEBUG
1735 std::vector<const CrystVector_REAL*> v;
1737 cout<<
"PowderPatternDiffraction::CalcPowderPatternIntegrated_FullDeriv():parameters:"<<endl;
1738 for(std::set<RefinablePar*>::iterator par=vPar.begin();par!=vPar.end();++par)
1740 if(mPowderPatternIntegrated_FullDeriv[*par].size()==0)
continue;
1741 v.push_back(&(mPowderPatternIntegrated_FullDeriv[*par]));
1742 v.push_back(&(oldDeriv[*par]));
1743 cout<<(*par)->GetName()<<
":"<<mPowderPatternIntegrated_FullDeriv[*par].size()<<
","<<oldDeriv[*par].size()<<endl;
1746 cout<<
"PowderPatternDiffraction::CalcPowderPatternIntegrated_FullDeriv():"<<endl<<FormatVertVector<REAL>(v,12,1,20)<<endl;
1753 this->CalcSinThetaLambda();
1755 if( (mClockProfileCalc>mClockProfilePar)
1756 &&(mClockProfileCalc>mpReflectionProfile->GetClockMaster())
1757 &&(mClockProfileCalc>mClockTheta)
1758 &&(mClockProfileCalc>this->GetRadiation().GetClockWavelength())
1760 &&(mClockProfileCalc>mClockHKL)
1763 TAU_PROFILE(
"PowderPatternDiffraction::CalcPowderReflProfile()",
"void (bool)",TAU_DEFAULT);
1764 VFN_DEBUG_ENTRY(
"PowderPatternDiffraction::CalcPowderReflProfile()",5)
1770 unsigned int nbLine=1;
1771 CrystVector_REAL spectrumDeltaLambdaOvLambda;
1772 CrystVector_REAL spectrumFactor;
1773 switch(this->GetRadiation().GetWavelengthType())
1775 case WAVELENGTH_MONOCHROMATIC:
1777 spectrumDeltaLambdaOvLambda.resize(1);spectrumDeltaLambdaOvLambda=0.0;
1778 spectrumFactor.resize(1);spectrumFactor=1.0;
1781 case WAVELENGTH_ALPHA12:
1784 spectrumDeltaLambdaOvLambda.resize(2);
1785 spectrumDeltaLambdaOvLambda(0)
1786 =-this->GetRadiation().GetXRayTubeDeltaLambda()
1787 *this->GetRadiation().GetXRayTubeAlpha2Alpha1Ratio()
1788 /(1+this->GetRadiation().GetXRayTubeAlpha2Alpha1Ratio())
1789 /this->GetRadiation().GetWavelength()(0);
1790 spectrumDeltaLambdaOvLambda(1)
1791 = this->GetRadiation().GetXRayTubeDeltaLambda()
1792 /(1+this->GetRadiation().GetXRayTubeAlpha2Alpha1Ratio())
1793 /this->GetRadiation().GetWavelength()(0);
1795 spectrumFactor.resize(2);
1796 spectrumFactor(0)=1./(1.+this->GetRadiation().GetXRayTubeAlpha2Alpha1Ratio());
1797 spectrumFactor(1)=this->GetRadiation().GetXRayTubeAlpha2Alpha1Ratio()
1798 /(1.+this->GetRadiation().GetXRayTubeAlpha2Alpha1Ratio());
1801 case WAVELENGTH_TOF:
1803 spectrumDeltaLambdaOvLambda.resize(1);spectrumDeltaLambdaOvLambda=0.0;
1804 spectrumFactor.resize(1);spectrumFactor=1.0;
1807 default:
throw ObjCrystException(
"PowderPatternDiffraction::PrepareIntegratedProfile():\
1808 Radiation must be either monochromatic, from an X-Ray Tube, or neutron TOF !!");
1812 VFN_DEBUG_MESSAGE(
"PowderPatternDiffraction::CalcPowderReflProfile():\
1813 Computing all Profiles",5)
1817 CrystVector_REAL vx,reflProfile,tmpV;
1818 mvReflProfile.resize(this->GetNbRefl());
1819 for(
unsigned int i=0;i<this->GetNbRefl();i++)
1821 mvReflProfile[i].first=0;
1822 mvReflProfile[i].last=0;
1823 mvReflProfile[i].profile.resize(0);
1825 VFN_DEBUG_MESSAGE(
"PowderPatternDiffraction::CalcPowderReflProfile()",5)
1827 for(
unsigned int line=0;line<nbLine;line++)
1829 for(
long i=0;i<this->GetNbRefl();i++)
1831 VFN_DEBUG_ENTRY(
"PowderPatternDiffraction::CalcPowderReflProfile()#"<<i,5)
1834 VFN_DEBUG_MESSAGE(
"PowderPatternDiffraction::CalcPowderReflProfile()#"<<i,5)
1838 x0+2*tan(x0/2.0)*spectrumDeltaLambdaOvLambda(line));
1842 if(!mUseFastLessPreciseFunc) fact=5.0;
1843 const REAL halfwidth=mpReflectionProfile->GetFullProfileWidth(0.04,center,mH(i),mK(i),mL(i))*fact;
1848 label<<mIntH(i)<<
" "<<mIntK(i)<<
" "<<mIntL(i);
1849 mvLabel.push_back(make_pair(center,label.str()));
1850 REAL spectrumwidth=0.0;
1851 if(this->GetRadiation().GetWavelengthType()==WAVELENGTH_ALPHA12)
1853 spectrumwidth=2*this->GetRadiation().GetXRayTubeDeltaLambda()
1854 /this->GetRadiation().GetWavelength()(0)*tan(x0/2.0);
1858 if(this->GetRadiation().GetWavelengthType()==WAVELENGTH_TOF)
1866 cout<<
"PowderPatternDiffraction::CalcPowderReflProfile(), line"<<__LINE__<<
"first>last !! :"<<first<<
","<<last<<endl;
1867 first=(first+last)/2;
1872 VFN_DEBUG_MESSAGE(
"PowderPatternDiffraction::CalcPowderReflProfile():"<<first<<
","<<last<<
","<<center,3)
1875 cout<<__FILE__<<__LINE__<<endl;
1880 if(first<0) first=0;
1883 vx.resize(last-first+1);
1886 mvReflProfile[i].first=first;
1887 mvReflProfile[i].last=last;
1891 first=mvReflProfile[i].first;
1892 last=mvReflProfile[i].last;
1894 vx.resize(last-first+1);
1896 vx.resize(last-first+1);
1903 for(
long i=first;i<=last;i++) *p1++ = *p0++;
1906 VFN_DEBUG_MESSAGE(
"PowderPatternDiffraction::CalcPowderReflProfile():"<<first<<
","<<last<<
","<<center,3)
1907 reflProfile=mpReflectionProfile->GetProfile(vx,center,mH(i),mK(i),mL(i));
1908 VFN_DEBUG_MESSAGE(
"PowderPatternDiffraction::CalcPowderReflProfile()",2)
1909 if(nbLine>1) reflProfile *=spectrumFactor(line);
1910 if(line==0) mvReflProfile[i].profile = reflProfile;
1911 else mvReflProfile[i].profile += reflProfile;
1915 mvReflProfile[i].profile.resize(0);
1917 VFN_DEBUG_EXIT(
"PowderPatternDiffraction::CalcPowderReflProfile():\
1918 Computing all Profiles: Reflection #"<<i,5)
1922 mClockProfileCalc.Click();
1923 VFN_DEBUG_EXIT(
"PowderPatternDiffraction::CalcPowderReflProfile()",5)
1928 TAU_PROFILE(
"PowderPatternDiffraction::CalcPowderReflProfile_FullDeriv()",
"void (bool)",TAU_DEFAULT);
1929 cout<<__FILE__<<
":"<<__LINE__<<
":PowderPatternDiffraction::CalcPowderReflProfile_FullDeriv()"<<endl;
1930 this->CalcPowderReflProfile();
1931 unsigned int nbLine=1;
1932 CrystVector_REAL spectrumDeltaLambdaOvLambda;
1933 CrystVector_REAL spectrumFactor;
1934 switch(this->GetRadiation().GetWavelengthType())
1936 case WAVELENGTH_MONOCHROMATIC:
1938 spectrumDeltaLambdaOvLambda.resize(1);spectrumDeltaLambdaOvLambda=0.0;
1939 spectrumFactor.resize(1);spectrumFactor=1.0;
1942 case WAVELENGTH_ALPHA12:
1945 spectrumDeltaLambdaOvLambda.resize(2);
1946 spectrumDeltaLambdaOvLambda(0)
1947 =-this->GetRadiation().GetXRayTubeDeltaLambda()
1948 *this->GetRadiation().GetXRayTubeAlpha2Alpha1Ratio()
1949 /(1+this->GetRadiation().GetXRayTubeAlpha2Alpha1Ratio())
1950 /this->GetRadiation().GetWavelength()(0);
1951 spectrumDeltaLambdaOvLambda(1)
1952 = this->GetRadiation().GetXRayTubeDeltaLambda()
1953 /(1+this->GetRadiation().GetXRayTubeAlpha2Alpha1Ratio())
1954 /this->GetRadiation().GetWavelength()(0);
1956 spectrumFactor.resize(2);
1957 spectrumFactor(0)=1./(1.+this->GetRadiation().GetXRayTubeAlpha2Alpha1Ratio());
1958 spectrumFactor(1)=this->GetRadiation().GetXRayTubeAlpha2Alpha1Ratio()
1959 /(1.+this->GetRadiation().GetXRayTubeAlpha2Alpha1Ratio());
1962 case WAVELENGTH_TOF:
1964 spectrumDeltaLambdaOvLambda.resize(1);spectrumDeltaLambdaOvLambda=0.0;
1965 spectrumFactor.resize(1);spectrumFactor=1.0;
1968 default:
throw ObjCrystException(
"PowderPatternDiffraction::CalcPowderReflProfile_FullDeriv():\
1969 Radiation must be either monochromatic, from an X-Ray Tube, or neutron TOF !!");
1974 CrystVector_REAL vx,reflProfile,tmpV;
1977 vector<CrystVector_REAL> vReflProfile_DerivCenter(mNbReflUsed);
1979 mvReflProfile_FullDeriv.clear();
1980 for(std::set<RefinablePar*>::iterator par=vPar.begin();par!=vPar.end();++par)
1982 if(*par==0)
continue;
1983 if( (*par)->GetType()->IsDescendantFromOrSameAs(gpRefParTypeRadiation)
1984 ||(*par)->GetType()->IsDescendantFromOrSameAs(gpRefParTypeUnitCell)
1988 mvReflProfile_FullDeriv[*par].resize(mNbReflUsed);
1990 for(
unsigned int line=0;line<nbLine;line++)
1992 for(
long i=0;i<mNbReflUsed;i++)
1999 x0+2*tan(x0/2.0)*spectrumDeltaLambdaOvLambda(line));
2003 first=mvReflProfile[i].first;
2004 last=mvReflProfile[i].last;
2006 vx.resize(last-first+1);
2008 vx.resize(last-first+1);
2014 for(
long i=first;i<=last;i++) *p1++ = *p0++;
2021 const REAL step=(*par)->GetDerivStep();
2022 (*par)->Mutate(step);
2023 reflProfile=mpReflectionProfile->GetProfile(vx,center,mH(i),mK(i),mL(i));
2024 (*par)->Mutate(-2*step);
2025 reflProfile-=mpReflectionProfile->GetProfile(vx,center,mH(i),mK(i),mL(i));
2026 (*par)->Mutate(step);
2027 reflProfile/=2*step;
2034 const REAL step=(*par)->GetDerivStep();
2035 (*par)->Mutate(step);
2039 (*par)->Mutate(-2*step);
2043 (*par)->Mutate(step);
2050 if(vReflProfile_DerivCenter[i].size()==0)
2052 const REAL step=1e-4;
2053 vReflProfile_DerivCenter[i] =mpReflectionProfile->GetProfile(vx,center+step,mH(i),mK(i),mL(i));
2054 vReflProfile_DerivCenter[i]-=mpReflectionProfile->GetProfile(vx,center-step,mH(i),mK(i),mL(i));
2055 vReflProfile_DerivCenter[i]/=2*step;
2057 reflProfile=vReflProfile_DerivCenter[i];
2058 reflProfile*=dcenter;
2063 reflProfile.resize(0);
2066 if(reflProfile.size()>0)
2068 if(nbLine>1) reflProfile *=spectrumFactor(line);
2069 if(line==0) mvReflProfile_FullDeriv[*par][i] = reflProfile;
2070 else mvReflProfile_FullDeriv[*par][i] += reflProfile;
2081 bool needRecalc=
false;
2083 this->CalcSinThetaLambda();
2084 if((mClockIntensityCorr<mClockTheta)||(mClockIntensityCorr<this->GetClockNbReflBelowMaxSinThetaOvLambda())) needRecalc=
true;
2086 const CrystVector_REAL *mpCorr[6] = {0, 0, 0, 0, 0, 0};
2088 if(this->GetRadiation().GetWavelengthType()==WAVELENGTH_TOF)
2090 mpCorr[0]=&(mCorrTOF.GetCorr());
2091 if(mClockIntensityCorr<mCorrTOF.GetClockCorr()) needRecalc=
true;
2095 mpCorr[0]=&(mCorrLorentz.GetCorr());
2096 if(mClockIntensityCorr<mCorrLorentz.GetClockCorr()) needRecalc=
true;
2098 if(this->GetRadiation().GetRadiationType()==RAD_XRAY)
2100 mpCorr[1]=&(mCorrPolar.GetCorr());
2101 if(mClockIntensityCorr<mCorrPolar.GetClockCorr()) needRecalc=
true;
2104 mpCorr[2]=&(mCorrSlitAperture.GetCorr());
2105 if(mClockIntensityCorr<mCorrSlitAperture.GetClockCorr()) needRecalc=
true;
2109 mpCorr[5] = &(mCorrCylAbs.GetCorr());
2110 if(mClockIntensityCorr<mCorrCylAbs.GetClockCorr()) needRecalc=
true;
2114 if(mCorrTextureMarchDollase.GetNbPhase()>0)
2116 mpCorr[3]=&(mCorrTextureMarchDollase.GetCorr());
2117 if(mClockIntensityCorr<mCorrTextureMarchDollase.GetClockCorr()) needRecalc=
true;
2119 mpCorr[4]=&(mCorrTextureEllipsoid.GetCorr());
2120 if(mClockIntensityCorr<mCorrTextureEllipsoid.GetClockCorr()) needRecalc=
true;
2122 if(needRecalc==
false)
return;
2124 TAU_PROFILE(
"PowderPatternDiffraction::CalcIntensityCorr()",
"void ()",TAU_DEFAULT);
2125 VFN_DEBUG_MESSAGE(
"PowderPatternDiffraction::CalcIntensityCorr()",2)
2126 mIntensityCorr.resize(mNbRefl);
2127 REAL *pCorr=mIntensityCorr.data();
2128 const REAL *p=mpCorr[0]->data();
2129 for(
long i=mNbReflUsed;i>0;i--) *pCorr++ = *p++;
2130 if(this->GetRadiation().GetWavelengthType()!=WAVELENGTH_TOF)
2132 if(this->GetRadiation().GetRadiationType()==RAD_XRAY)
2134 pCorr=mIntensityCorr.data();
2135 p=mpCorr[1]->data();
2136 const REAL* p2=mpCorr[2]->data();
2137 for(
long i=mNbReflUsed;i>0;i--) *pCorr++ *= *p++ * *p2++;
2141 pCorr=mIntensityCorr.data();
2142 p=mpCorr[2]->data();
2143 for(
long i=mNbReflUsed;i>0;i--) *pCorr++ *= *p++;
2146 if(mCorrTextureMarchDollase.GetNbPhase()>0)
2148 pCorr=mIntensityCorr.data();
2149 p=mpCorr[3]->data();
2150 for(
long i=mNbReflUsed;i>0;i--) *pCorr++ *= *p++;
2152 if(mpCorr[4]->numElements()>0)
2154 pCorr=mIntensityCorr.data();
2155 p=mpCorr[4]->data();
2156 for(
long i=mNbReflUsed;i>0;i--) *pCorr++ *= *p++;
2158 if(mpCorr[5] != NULL)
2159 if(mpCorr[5]->numElements()>0)
2161 pCorr=mIntensityCorr.data();
2162 p=mpCorr[5]->data();
2163 for(
long i=mNbReflUsed;i>0;i--) *pCorr++ *= *p++;
2165 mClockIntensityCorr.Click();
2166 VFN_DEBUG_MESSAGE(
"PowderPatternDiffraction::CalcIntensityCorr():finished",2)
2171 this->CalcIntensityCorr();
2172 if(mExtractionMode==
true)
2174 VFN_DEBUG_MESSAGE(
"PowderPatternDiffraction::CalcIhkl():"<<mFhklObsSq.numElements()<<
","<<mIntensityCorr.numElements()<<
","<<mMultiplicity.numElements(),7);
2175 mIhklCalc=mFhklObsSq;
2176 mIhklCalc*=mIntensityCorr;
2177 mIhklCalc*=mMultiplicity;
2178 mClockIhklCalc.Click();
2181 this->CalcStructFactor();
2182 if( (mClockIhklCalc>mClockIntensityCorr)
2183 &&(mClockIhklCalc>mClockStructFactor)
2184 &&(mClockIhklCalc>mClockNbReflUsed))
return;
2186 VFN_DEBUG_MESSAGE(
"PowderPatternDiffraction::CalcIhkl()",3)
2187 TAU_PROFILE(
"PowderPatternDiffraction::CalcIhkl()",
"void ()",TAU_DEFAULT);
2188 const REAL * RESTRICT pr,* RESTRICT pi,* RESTRICT pcorr;
2189 const int * RESTRICT mult;
2192 pr=mFhklCalcReal.data();
2193 pi=mFhklCalcImag.data();
2194 pcorr=mIntensityCorr.data();
2196 mult=mMultiplicity.data();
2197 mIhklCalc.resize(mNbRefl);
2199 if(mFhklCalcVariance.numElements()>0)
2201 const REAL * RESTRICT pv=mFhklCalcVariance.data();
2202 for(
long i=mNbReflUsed;i>0;i--)
2204 *p++ = *mult++ * (*pr * *pr + *pi * *pi + 2 * *pv++) * *pcorr++;
2211 for(
long i=mNbReflUsed;i>0;i--)
2213 *p++ = *mult++ * (*pr * *pr + *pi * *pi) * *pcorr++;
2218 if(mFhklCalcVariance.numElements()==0)
2220 VFN_DEBUG_MESSAGE(
"PowderPatternDiffraction::CalcIhkl(): No Calc Variance",2)
2221 mIhklCalcVariance.resize(0);
2222 VFN_DEBUG_MESSAGE(endl<<
2232 VFN_DEBUG_MESSAGE(
"PowderPatternDiffraction::CalcIhkl(): Calc Variance",2)
2233 mIhklCalcVariance.resize(mNbRefl);
2234 REAL * RESTRICT pVar2=mIhklCalcVariance.data();
2236 const REAL * RESTRICT pInt=mIhklCalc.data();
2237 const REAL * RESTRICT pVar=mFhklCalcVariance.data();
2238 pcorr=mIntensityCorr.data();
2239 mult=mMultiplicity.data();
2241 for(
long j=mNbReflUsed;j>0;j--)
2243 *pVar2++ = (4* *mult) * *pcorr * *pVar *(*pInt++ - (*mult * *pcorr) * *pVar);
2244 pVar++;mult++;pcorr++;
2246 VFN_DEBUG_MESSAGE(endl<<
2250 mExpectedIntensityFactor,
2253 mvLuzzatiFactor[&(mpCrystal->GetScatteringPowerRegistry().GetObj(0))],
2255 mIhklCalcVariance),2);
2256 VFN_DEBUG_MESSAGE(mNbRefl<<
" "<<mNbReflUsed,2)
2261 mClockIhklCalc.Click();
2262 VFN_DEBUG_MESSAGE(
"PowderPatternDiffraction::CalcIhkl():End",3)
2265 void PowderPatternDiffraction::CalcIhkl_FullDeriv(std::set<RefinablePar*> &vPar)
2267 TAU_PROFILE(
"PowderPatternDiffraction::CalcIhkl_FullDeriv()",
"void ()",TAU_DEFAULT);
2269 this->CalcIntensityCorr();
2270 mIhkl_FullDeriv.clear();
2271 if(mExtractionMode==
true)
2276 this->CalcStructFactor_FullDeriv(vPar);
2278 for(std::set<RefinablePar*>::iterator par=vPar.begin();par!=vPar.end();++par)
2280 if(*par==0) mIhkl_FullDeriv[*par]=mIhklCalc;
2283 if(mFhklCalcReal_FullDeriv[*par].size()==0)
continue;
2284 const REAL * RESTRICT pr,* RESTRICT pi,* RESTRICT prd,* RESTRICT pid,* RESTRICT pcorr;
2285 const int * RESTRICT mult;
2288 pr=mFhklCalcReal.data();
2289 pi=mFhklCalcImag.data();
2290 prd=mFhklCalcReal_FullDeriv[*par].data();
2291 pid=mFhklCalcImag_FullDeriv[*par].data();
2292 pcorr=mIntensityCorr.data();
2294 mult=mMultiplicity.data();
2295 mIhkl_FullDeriv[*par].resize(mNbRefl);
2296 p=mIhkl_FullDeriv[*par].data();
2297 if(mFhklCalcImag_FullDeriv[*par].size()==0)
2298 for(
long i=mNbReflUsed;i>0;i--) *p++ = *mult++ * 2 * *pr++ * *prd++ * *pcorr++;
2300 for(
long i=mNbReflUsed;i>0;i--) *p++ = *mult++ * 2 *(*pr++ * *prd++ + *pi++ * *pid++) * *pcorr++;
2304 std::map<RefinablePar*, CrystVector_REAL> oldDeriv;
2305 std::vector<const CrystVector_REAL*> v;
2312 v.push_back(&mIntensityCorr);
2314 for(std::set<RefinablePar*>::iterator par=vPar.begin();par!=vPar.end();++par)
2316 if((*par)==0)
continue;
2317 if(mIhkl_FullDeriv[*par].size()==0)
continue;
2319 const REAL step=(*par)->GetDerivStep();
2320 (*par)->Mutate(step);
2322 oldDeriv[*par]=mIhklCalc;
2323 (*par)->Mutate(-2*step);
2325 oldDeriv[*par]-=mIhklCalc;
2326 oldDeriv[*par]/=2*step;
2327 (*par)->Mutate(step);
2329 v.push_back(&(mIhkl_FullDeriv[*par]));
2330 v.push_back(&(oldDeriv[*par]));
2331 cout<<(*par)->GetName()<<
":"<<mIhkl_FullDeriv[*par].size()<<
","<<oldDeriv[*par].size()<<
", step="<<setw(16)<<step<<endl;
2334 cout<<
"PowderPatternDiffraction::CalcIhkl_FullDeriv():"<<endl<<FormatVertVectorHKLFloats<REAL>(v,12,1,20)<<endl;
2341 if( (this->GetCrystal().GetSpaceGroup().GetClockSpaceGroup()>mClockHKL)
2342 ||(this->GetCrystal().GetClockLatticePar()>mClockHKL)
2343 ||(this->GetRadiation().GetClockWavelength()>mClockHKL)
2345 this->GenHKLFullSpace();
2348 void PowderPatternDiffraction::InitOptions()
2350 VFN_DEBUG_MESSAGE(
"PowderPatternDiffraction::InitOptions()",5)
2352 static string ReflectionProfileTypeName;
2353 static string ReflectionProfileTypeChoices[3];
2355 static bool needInitNames=
true;
2356 if(
true==needInitNames)
2358 ReflectionProfileTypeName=
"Profile Type";
2359 ReflectionProfileTypeChoices[0]=
"Gaussian";
2360 ReflectionProfileTypeChoices[1]=
"Lorentzian";
2361 ReflectionProfileTypeChoices[2]=
"Pseudo-Voigt";
2363 needInitNames=
false;
2365 mReflectionProfileType.Init(3,&ReflectionProfileTypeName,ReflectionProfileTypeChoices);
2366 this->
AddOption(&mReflectionProfileType);
2371 this->CalcPowderReflProfile();
2372 if((mClockProfileCalc>
mClockBraggLimits)&&(this->GetNbReflBelowMaxSinThetaOvLambda()>0))
2374 VFN_DEBUG_ENTRY(
"PowderPatternDiffraction::GetBraggLimits(*min,*max)",3)
2375 TAU_PROFILE(
"PowderPatternDiffraction::GetBraggLimits()",
"void ()",TAU_DEFAULT);
2379 for(;i<(this->GetNbReflBelowMaxSinThetaOvLambda()-1);++i)
2380 mIntegratedReflLimits(i+1)=(mvReflProfile[i].first+mvReflProfile[i].last+mvReflProfile[i+1].first+mvReflProfile[i+1].last)/4;
2383 VFN_DEBUG_EXIT(
"PowderPatternDiffraction::GetBraggLimits(*min,*max)",3)
2393 if(mFreezeLatticePar)
return mFrozenBMatrix;
2399 VFN_DEBUG_MESSAGE(
"PowderPatternDiffraction::CalcFrozenBMatrix()", 10)
2400 REAL a,b,c,alpha,beta,gamma;
2401 REAL aa,bb,cc,alphaa,betaa,gammaa;
2402 POSSIBLY_UNUSED(alphaa);
2404 a=mFrozenLatticePar(0);
2405 b=mFrozenLatticePar(1);
2406 c=mFrozenLatticePar(2);
2407 alpha=mFrozenLatticePar(3);
2408 beta=mFrozenLatticePar(4);
2409 gamma=mFrozenLatticePar(5);
2411 v=sqrt(1-cos(alpha)*cos(alpha)-cos(beta)*cos(beta)-cos(gamma)*cos(gamma)
2412 +2*cos(alpha)*cos(beta)*cos(gamma));
2418 alphaa=acos( (cos(beta )*cos(gamma)-cos(alpha))/sin(beta )/sin(gamma) );
2419 betaa =acos( (cos(alpha)*cos(gamma)-cos(beta ))/sin(alpha)/sin(gamma) );
2420 gammaa=acos( (cos(alpha)*cos(beta )-cos(gamma))/sin(alpha)/sin(beta ) );
2422 mFrozenBMatrix = aa , bb*cos(gammaa) , cc*cos(betaa) ,
2423 0 , bb*sin(gammaa) ,-cc*sin(betaa)*cos(alpha),
2427 void PowderPatternDiffraction::PrepareIntegratedProfile()
const
2429 this->CalcPowderReflProfile();
2430 this->GetNbReflBelowMaxSinThetaOvLambda();
2432 if( (mClockIntegratedProfileFactor>mClockProfileCalc)
2434 &&(mClockIntegratedProfileFactor>mClockNbReflUsed))
2436 VFN_DEBUG_ENTRY(
"PowderPatternDiffraction::PrepareIntegratedProfile()",7)
2441 const
long numInterval=pMin->numElements();
2443 vector< map<
long, REAL> > vIntegratedProfileFactor;
2444 vIntegratedProfileFactor.resize(mNbReflUsed);
2445 vector< map<
long, REAL> >::iterator pos1;
2446 pos1=vIntegratedProfileFactor.begin();
2448 mIntegratedProfileFactor.resize(mNbReflUsed);
2449 vector< pair<
unsigned long, CrystVector_REAL> >::iterator pos2;
2450 pos2=mIntegratedProfileFactor.begin();
2451 for(
long i=0;i<mNbReflUsed;i++)
2454 long firstInterval=numInterval;
2455 for(
long j=0;j<numInterval;j++)
2457 const long first0 = mvReflProfile[i].first;
2458 const long last0 = mvReflProfile[i].last ;
2459 const long first= first0>(*pMin)(j) ? first0:(*pMin)(j);
2460 const long last = last0 <(*pMax)(j) ? last0 :(*pMax)(j);
2461 if((first<=last) && (mvReflProfile[i].profile.size()>0))
2463 if(firstInterval>j) firstInterval=j;
2464 if(pos1->find(j) == pos1->end()) (*pos1)[j]=0.;
2465 REAL *fact = &((*pos1)[j]);
2466 const REAL *p2 = mvReflProfile[i].profile.data()+(first-first0);
2468 for(
long k=first;k<=last;k++) *fact += *p2++;
2471 pos2->first=firstInterval;
2472 pos2->second.resize(pos1->size());
2473 REAL *pFact=pos2->second.data();
2474 for(map<long, REAL>::const_iterator pos=pos1->begin();pos!=pos1->end();++pos)
2475 *pFact++ = pos->second;
2479 mClockIntegratedProfileFactor.Click();
2481 if(gVFNDebugMessageLevel<3)
2484 for(vector< pair<unsigned long, CrystVector_REAL> >::const_iterator
2485 pos=mIntegratedProfileFactor.begin();
2486 pos!=mIntegratedProfileFactor.end();++pos)
2488 cout <<
"Integrated profile factors for reflection #"<<i++<<
" ";
2489 for(
int j=0;j<pos->second.numElements();++j)
2490 cout << j+pos->first<<
"("<<pos->second(j)<<
") ";
2496 VFN_DEBUG_EXIT(
"PowderPatternDiffraction::PrepareIntegratedProfile()",7)
2499 #ifdef __WX__CRYST__
2500 WXCrystObjBasic* PowderPatternDiffraction::WXCreate(wxWindow* parent)
2503 mpWXCrystObj=
new WXPowderPatternDiffraction(parent,
this);
2504 return mpWXCrystObj;
2513 ObjRegistry<PowderPattern>
2514 gPowderPatternRegistry(
"List of all PowderPattern objects");
2516 PowderPattern::PowderPattern():
2517 mIsXAscending(true),mNbPoint(0),
2518 mXZero(0.),m2ThetaDisplacement(0.),m2ThetaTransparency(0.),
2519 mDIFC(48277.14),mDIFA(-6.7),
2520 mScaleFactor(20),mMuR(0), mUseFastLessPreciseFunc(false),
2525 mPowderPatternComponentRegistry.SetName(
"Powder Pattern Components");
2528 gPowderPatternRegistry.Register(*
this);
2538 PowderPattern::PowderPattern(
const PowderPattern &old):
2539 mIsXAscending(old.mIsXAscending),mNbPoint(old.mNbPoint),
2540 mRadiation(old.mRadiation),
2541 mXZero(old.mXZero),m2ThetaDisplacement(old.m2ThetaDisplacement),
2542 m2ThetaTransparency(old.m2ThetaTransparency),
2543 mDIFC(old.mDIFC),mDIFA(old.mDIFA),
2544 mPowderPatternComponentRegistry(old.mPowderPatternComponentRegistry),
2545 mScaleFactor(old.mScaleFactor),mMuR(old.mMuR),
2546 mUseFastLessPreciseFunc(old.mUseFastLessPreciseFunc),
2547 mStatisticsExcludeBackground(old.mStatisticsExcludeBackground),
2548 mMaxSinThetaOvLambda(old.mMaxSinThetaOvLambda),mNbPointUsed(old.mNbPointUsed)
2553 gPowderPatternRegistry.Register(*
this);
2563 PowderPattern::~PowderPattern()
2565 gPowderPatternRegistry.DeRegister(*
this);
2566 for(
int i=0;i<mPowderPatternComponentRegistry.GetNb();i++)
2568 mPowderPatternComponentRegistry.GetObj(i).DeRegisterClient(*
this);
2570 delete &(mPowderPatternComponentRegistry.GetObj(i));
2576 const static string className=
"PowderPattern";
2582 VFN_DEBUG_ENTRY(
"PowderPattern::AddPowderPatternComponent():"<<comp.
GetName(),5)
2587 mClockIntegratedFactorsPrep.Reset();
2588 mPowderPatternComponentRegistry.Register(comp);
2591 mScaleFactor(mPowderPatternComponentRegistry.GetNb()-1)=1.;
2595 RefinablePar tmp(
"Scale_"+comp.
GetName(),mScaleFactor.data()+mPowderPatternComponentRegistry.GetNb()-1,
2597 false,
true,
true,
false,1.);
2604 VFN_DEBUG_EXIT(
"PowderPattern::AddPowderPatternComponent():"<<comp.
GetName(),5)
2609 VFN_DEBUG_ENTRY(
"PowderPattern::RemovePowderPatternComponent():"<<comp.
GetName(),5)
2613 cout<<
"PowderPattern::RemovePowderPatternComponent: removing 1 scale paramater"<<endl;
2615 this->
RemovePar(&this->
GetPar(mScaleFactor.data()+mPowderPatternComponentRegistry.GetNb()-1));
2622 mClockIntegratedFactorsPrep.Reset();
2623 mPowderPatternComponentRegistry.DeRegister(comp);
2633 VFN_DEBUG_EXIT(
"PowderPattern::RemovePowderPatternComponent():"<<comp.
GetName(),5)
2638 VFN_DEBUG_ENTRY(
"PowderPattern::RemovePowderPatternComponent():"<<i,5)
2640 VFN_DEBUG_EXIT(
"PowderPattern::RemovePowderPatternComponent():"<<i,5)
2646 return mPowderPatternComponentRegistry.GetNb();
2650 (
const string &name)
const
2652 return mPowderPatternComponentRegistry.GetObj(name);
2658 return mPowderPatternComponentRegistry.GetObj(i);
2662 (
const string &name)
2664 return mPowderPatternComponentRegistry.GetObj(name);
2670 return mPowderPatternComponentRegistry.GetObj(i);
2678 for(;i<mPowderPatternComponentRegistry.GetNb();++i)
2680 if(&(mPowderPatternComponentRegistry.GetObj(i))==&comp)
break;
2682 if(i==mPowderPatternComponentRegistry.GetNb())
2683 throw ObjCrystException(
"PowderPattern::GetScaleFactor(comp) : no such component");
2684 return mScaleFactor(i);
2691 for(;i<mPowderPatternComponentRegistry.GetNb();++i)
2693 if(&(mPowderPatternComponentRegistry.GetObj(i))==&comp)
break;
2695 if(i==mPowderPatternComponentRegistry.GetNb())
2696 throw ObjCrystException(
"PowderPattern::GetScaleFactor(comp) : no such component");
2702 unsigned long nbPoint)
2704 VFN_DEBUG_MESSAGE(
"PowderPattern::SetPowderPatternPar():"<<min<<
","<<step<<
","<<nbPoint,3)
2707 for(
unsigned long i=0;i<
mNbPoint;i++)
mX(i)=min+step*i;
2723 VFN_DEBUG_MESSAGE(
"PowderPattern::SetPowderPatternX() is ascending="<<
mIsXAscending,5)
2731 return mNbPointUsed;
2753 VFN_DEBUG_MESSAGE(
"PowderPattern::SetWavelength(lambda)",3)
2759 VFN_DEBUG_MESSAGE(
"PowderPattern::SetWavelength(wavelength)",3)
2771 std::map<RefinablePar*,CrystVector_REAL>& PowderPattern::GetPowderPattern_FullDeriv(std::set<RefinablePar *> &vPar)
2773 this->CalcPowderPattern_FullDeriv(vPar);
2774 return mPowderPattern_FullDeriv;
2823 VFN_DEBUG_ENTRY(
"PowderPattern::GetChi2Cumul()",3)
2826 if((mode!=0) && (mode!=1)) mode = mOptProfileIntegration.GetChoice();
2830 if(mNbIntegrationUsed==0)
2834 const REAL *pObs=mIntegratedObs.data();
2836 const REAL *pWeight;
2837 if(mIntegratedWeight.numElements()==0) pWeight=mIntegratedWeightObs.data();
2838 else pWeight=mIntegratedWeight.data();
2841 for(
int i=0;i<mIntegratedPatternMin(0);i++) *pC2Cu++ = 0;
2842 REAL chi2cumul=0,tmp;
2843 for(
unsigned long j=1;j<mNbIntegrationUsed;j++)
2845 tmp=(*pObs++ - *pCalc++) ;
2846 chi2cumul += *pWeight++ * tmp*tmp;
2847 for(
int i=mIntegratedPatternMin(j-1);i<mIntegratedPatternMin(j);i++) *pC2Cu++ =chi2cumul;
2849 if(mIntegratedPatternMin(j)>(int)mNbPointUsed)
2851 for(
unsigned int i=mIntegratedPatternMin(j);i<
mNbPoint;i++) *pC2Cu++ =chi2cumul;
2855 pC2Cu=
mChi2Cumul.data()+mIntegratedPatternMin(mNbIntegrationUsed-1);
2856 for(
unsigned int i=mIntegratedPatternMin(mNbIntegrationUsed-1);i<
mNbPoint;i++) *pC2Cu++ =chi2cumul;
2866 REAL chi2cumul=0,tmp;
2867 for(
unsigned int i=0;i<mNbPointUsed;i++)
2869 tmp = (*pObs++ - *pCalc++) ;
2870 chi2cumul += *pWeight++ * tmp*tmp;
2871 *pC2Cu++ = chi2cumul;
2874 VFN_DEBUG_EXIT(
"PowderPattern::GetChi2Cumul()",3)
2901 m2ThetaDisplacement=displacement;
2907 m2ThetaTransparency=transparency;
2916 x += m2ThetaDisplacement*cos(x/2) +m2ThetaTransparency*sin(x);
2933 VFN_DEBUG_MESSAGE(
"PowderPattern::X2Pixel()",1)
2935 if((pix>0)&&(pix<((long)
mNbPoint-1)))
2938 const REAL localStep=
mX(pix)-
mX(pix+1);
2939 if(localStep>0) pix -= (long)((x-
mX(pix))/localStep);
2941 VFN_DEBUG_MESSAGE(
"PowderPattern::X2Pixel():"<<x<<
","<<pix,1)
2944 VFN_DEBUG_MESSAGE(
"PowderPattern::X2Pixel():"<<x<<
","<<pix<<
","<<
mX(pix),1)
2949 VFN_DEBUG_MESSAGE(
"PowderPattern::X2Pixel():"<<x<<
","<<pix<<
","<<
mX(pix),1)
2950 if(
mX(pix)>=x)
break;
2958 if(
mX(pix)<=x) {pix--;
break;}
2963 VFN_DEBUG_MESSAGE(
"PowderPattern::X2Pixel():"<<x<<
","<<pix<<
","<<
mX(pix),1)
2964 const REAL localStep=
mX(pix)-
mX(pix+1);
2965 pixx = (REAL)pix-(x-
mX(pix))/localStep;
2966 VFN_DEBUG_MESSAGE(
"PowderPattern::X2Pixel():"<<x<<
","<<pix<<
","<<
mX(pix),1)
2972 if((pix>0)&&(pix<((long)
mNbPoint-1)))
2975 const REAL localStep=
mX(pix+1)-
mX(pix);
2976 if(localStep>0) pix += (long)((x-
mX(pix))/localStep);
2978 VFN_DEBUG_MESSAGE(
"PowderPattern::X2Pixel():"<<x<<
","<<pix,1)
2981 VFN_DEBUG_MESSAGE(
"PowderPattern::X2Pixel():"<<x<<
","<<pix<<
","<<
mX(pix),1)
2986 VFN_DEBUG_MESSAGE(
"PowderPattern::X2Pixel():"<<x<<
","<<pix<<
","<<
mX(pix),1)
2987 if(
mX(pix)<=x)
break;
2995 VFN_DEBUG_MESSAGE(
"PowderPattern::X2Pixel():"<<x<<
","<<pix<<
","<<
mX(pix),1)
2996 if(
mX(pix)>=x) {pix-- ;
break;}
3000 VFN_DEBUG_MESSAGE(
"PowderPattern::X2Pixel():"<<x<<
","<<pix<<
","<<
mX(pix),1)
3003 const REAL localStep=
mX(pix+1)-
mX(pix);
3004 VFN_DEBUG_MESSAGE(
"PowderPattern::X2Pixel():"<<x<<
","<<pix<<
","<<
mX(pix)<<
","<<localStep,1)
3005 pixx = (REAL)pix+(x-
mX(pix))/localStep;
3007 VFN_DEBUG_MESSAGE(
"PowderPattern::X2Pixel():"<<x<<
","<<pixx,1)
3016 VFN_DEBUG_MESSAGE(
"PowderPattern::ImportPowderPatternFullprof() : \
3017 from file : "+filename,5)
3018 ifstream fin(filename.c_str());
3022 Error opening file for input:"+filename);
3025 fin >> min >> step >> max;
3030 VFN_DEBUG_MESSAGE(
"PowderPattern::ImportPowderPatternFullprof() :"\
3031 <<
" 2Theta min=" << min*RAD2DEG <<
" 2Theta max=" << max*RAD2DEG \
3037 char tmpComment[200];
3038 fin.getline(tmpComment,100);
3049 sprintf(buf,
"Imported powder pattern: %d points, 2theta=%7.3f -> %7.3f, step=%6.3f",
3050 (
int)
mNbPoint,min*RAD2DEG,max*RAD2DEG,step*RAD2DEG);
3051 (*fpObjCrystInformUser)((string)buf);
3053 VFN_DEBUG_MESSAGE(
"PowderPattern::ImportFullProfPattern():finished:"<<
mNbPoint<<
" points",5)
3058 VFN_DEBUG_MESSAGE(
"PowderPattern::ImportPowderPatternPSI_DMC() : \
3059 from file : "+filename,5)
3060 ifstream fin(filename.c_str());
3064 Error opening file for input:"+filename);
3067 char tmpComment[200];
3068 fin.getline(tmpComment,190);
3069 fin.getline(tmpComment,190);
3071 fin >> min >> step >> max;
3076 VFN_DEBUG_MESSAGE(
"PowderPattern::ImportPowderPatternPSI_DMC() :"\
3077 <<
" 2Theta min=" << min*RAD2DEG <<
" 2Theta max=" << max*RAD2DEG \
3083 fin.getline(tmpComment,100);
3094 sprintf(buf,
"Imported powder pattern: %d points, 2theta=%7.3f -> %7.3f, step=%6.3f",
3095 (
int)
mNbPoint,min*RAD2DEG,max*RAD2DEG,step*RAD2DEG);
3096 (*fpObjCrystInformUser)((string)buf);
3098 VFN_DEBUG_MESSAGE(
"PowderPattern::ImportPowderPatternPSI_DMC():finished",5)
3103 VFN_DEBUG_MESSAGE(
"PowderPattern::ImportPowderPatternILL_D1AD2B() : \
3104 from file : "+filename,5)
3105 ifstream fin(filename.c_str());
3109 Error opening file for input:"+filename);
3112 char tmpComment[200];
3113 fin.getline(tmpComment,190);
3114 fin.getline(tmpComment,190);
3115 fin.getline(tmpComment,190);
3118 fin.getline(tmpComment,190);
3124 VFN_DEBUG_MESSAGE(
"PowderPattern::ImportPowderPatternILL_D1AD2B() :"\
3125 <<
" 2Theta min=" << min*RAD2DEG <<
" 2Theta max=" << min*RAD2DEG+
mNbPoint*step*RAD2DEG \
3141 sprintf(buf,
"Imported powder pattern: %d points, 2theta=%7.3f -> %7.3f, step=%6.3f",
3144 (*fpObjCrystInformUser)((string)buf);
3146 VFN_DEBUG_MESSAGE(
"PowderPattern::ImportPowderPatternILL_D1AD2B():finished",5)
3151 VFN_DEBUG_MESSAGE(
"PowderPattern::ImportPowderPatternXdd():from file :" \
3153 ifstream fin (filename.c_str());
3157 Error opening file for input:"+filename);
3159 char tmpComment[200];
3160 fin.getline(tmpComment,100);
3162 REAL min,max,step,tmp;
3163 fin >> min >> step >> max;
3185 sprintf(buf,
"Imported powder pattern: %d points, 2theta=%7.3f -> %7.3f, step=%6.3f",
3186 (
int)
mNbPoint,min*RAD2DEG,max*RAD2DEG,step*RAD2DEG);
3187 (*fpObjCrystInformUser)((string)buf);
3189 VFN_DEBUG_MESSAGE(
"DiffractionDataPowder::ImportXddPattern() :finished",5)
3194 VFN_DEBUG_ENTRY(
"PowderPattern::ImportPowderPatternSietronicsCPI():from file :" \
3196 ifstream fin (filename.c_str());
3200 Error opening file for input:"+filename);
3202 char tmpComment[300];
3203 fin.getline(tmpComment,100);
3204 VFN_DEBUG_MESSAGE(
" ->Discarded comment :"<<tmpComment,1)
3206 fin >> min >> max >> step;
3222 VFN_DEBUG_MESSAGE(
" ->Read :"<<str,1)
3223 }
while (
"SCANDATA"!=str);
3233 sprintf(buf,
"Imported powder pattern: %d points, 2theta=%7.3f -> %7.3f, step=%6.3f",
3234 (
int)
mNbPoint,min*RAD2DEG,max*RAD2DEG,step*RAD2DEG);
3235 (*fpObjCrystInformUser)((string)buf);
3237 VFN_DEBUG_EXIT(
"DiffractionDataPowder::ImportPowderPatternSietronicsCPI()",5)
3242 VFN_DEBUG_MESSAGE(
"DiffractionDataPowder::ImportPowderPattern2ThetaObsSigma():from:" \
3244 ifstream fin (filename.c_str());
3248 Error opening file for input:"+filename);
3251 char tmpComment[200];
3252 for(
int i=0;i<nbSkip;i++) fin.getline(tmpComment,150);
3275 }
while(fin.eof() ==
false);
3290 sprintf(buf,
"Imported powder pattern: %d points, 2theta=%7.3f -> %7.3f, step=%6.3f",
3294 (*fpObjCrystInformUser)((string)buf);
3296 VFN_DEBUG_MESSAGE(
"PowderPattern::ImportPowderPattern2ThetaObsSigma()\
3297 :finished: "<<
mNbPoint<<
" points",5)
3302 VFN_DEBUG_MESSAGE(
"PowderPattern::ImportPowderPattern2ThetaObs():from:" \
3304 ifstream fin (filename.c_str());
3308 Error opening file for input:"+filename);
3311 char tmpComment[200];
3312 for(
int i=0;i<nbSkip;i++) fin.getline(tmpComment,150);
3333 }
while(fin.eof() ==
false);
3349 sprintf(buf,
"Imported powder pattern: %d points, 2theta=%7.3f -> %7.3f, step=%6.3f",
3353 (*fpObjCrystInformUser)((string)buf);
3355 VFN_DEBUG_MESSAGE(
"PowderPattern::ImportPowderPattern2ThetaObs():finished",5)
3375 VFN_DEBUG_MESSAGE(
"PowderPattern::ImportPowderPatternMultiDetectorLLBG42() : \
3376 from file : "+filename,5)
3377 ifstream fin(filename.c_str());
3380 throw ObjCrystException(
"PowderPattern::ImportPowderPatternMultiDetectorLLBG42() : \
3381 Error opening file for input:"+filename);
3388 fin >>junk>>junk>>step>>junk>>junk>>junk>>min>>junk>>junk>>junk>>junk;
3391 VFN_DEBUG_MESSAGE(
"PowderPattern::ImportPowderPatternMultiDetectorLLBG42() :"\
3392 <<
" 2Theta min=" << min*RAD2DEG <<
" 2Theta step=" << step*RAD2DEG,5)
3405 sscanf(str.c_str(),
"%f",&tmp);
3407 const unsigned int nb=str.length()/8;
3408 for(
unsigned int i=0;i<nb;i++)
3415 sub=str.substr(i*8,8);
3416 sscanf(sub.c_str(),
"%2f%6f",&ct,&iobs);
3433 sprintf(buf,
"Imported powder pattern: %d points, 2theta=%7.3f -> %7.3f, step=%6.3f",
3435 (*fpObjCrystInformUser)((string)buf);
3437 VFN_DEBUG_MESSAGE(
"PowderPattern::ImportPowderPatternMultiDetectorLLBG42():finished:"<<
mNbPoint<<
" points",5)
3447 VFN_DEBUG_MESSAGE(
"PowderPattern::ImportPowderPatternFullprof4() : \
3448 from file : "+filename,5)
3449 ifstream fin(filename.c_str());
3453 Error opening file for input:"+filename);
3456 fin >> min >> step >> max;
3461 VFN_DEBUG_MESSAGE(
"PowderPattern::ImportPowderPatternFullprof4() :"\
3462 <<
" 2Theta min=" << min*RAD2DEG <<
" 2Theta max=" << max*RAD2DEG \
3477 for(
unsigned int i=0;i<str.size();i++)
if(
' '==str[i]) str[i]=
'0';
3478 sscanf(str.c_str(),
"%8f%8f%8f%8f%8f%8f%8f%8f%8f%8f",
3479 line+0,line+1,line+2,line+3,line+4,line+5,line+6,line+7,line+8,line+9);
3480 for(
unsigned int j=0;j<10;j++)
3483 for(
unsigned int i=0;i<str.size();i++)
if(
' '==str[i]) str[i]=
'0';
3484 sscanf(str.c_str(),
"%8f%8f%8f%8f%8f%8f%8f%8f%8f%8f",
3485 line+0,line+1,line+2,line+3,line+4,line+5,line+6,line+7,line+8,line+9);
3486 for(
unsigned int j=0;j<10;j++)
3495 sprintf(buf,
"Imported powder pattern: %d points, 2theta=%7.3f -> %7.3f, step=%6.3f",
3496 (
int)
mNbPoint,min*RAD2DEG,max*RAD2DEG,
3498 (*fpObjCrystInformUser)((string)buf);
3500 VFN_DEBUG_MESSAGE(
"PowderPattern::ImportFullProfPattern4():finished:"<<
mNbPoint<<
" points",5)
3505 VFN_DEBUG_MESSAGE(
"DiffractionDataPowder::ImportPowderPatternTOF_ISIS_XYSigma():from:" \
3507 ifstream fin (filename.c_str());
3511 Error opening file for input:"+filename);
3515 fin.getline(junk,150);
3538 }
while(fin.eof() ==
false);
3548 for(
unsigned long i=0;i<(
mNbPoint/2);i++)
3571 sprintf(buf,
"Imported TOF powder pattern: %d points, TOF=%7.3f -> %7.3f",
3574 (*fpObjCrystInformUser)((string)buf);
3576 VFN_DEBUG_MESSAGE(
"PowderPattern::ImportPowderPatternTOF_ISIS_XYSigma()\
3577 :finished: "<<
mNbPoint<<
" points",5)
3582 VFN_DEBUG_ENTRY(
"PowderPattern::ImportPowderPatternGSAS():file:"<<filename,5)
3583 ifstream fin (filename.c_str());
3587 Error opening file for input:"+filename);
3592 while(isprint(fin.peek())==
false)
3594 if(fin.eof())
break;
3597 cout<<
"Title:"<<title<<endl;
3602 int numBank,nbRecords;
3603 string binType, type;
3610 fin.getline(line,80);
3611 while(isprint(fin.peek())==
false)
3613 if(fin.eof())
break;
3616 sscanf(line,
"%4s",bank);
3619 Could not find BANK statement !! In file: "+filename);
3621 while(
string(bank)!=
string(
"BANK"));
3625 char binTypeC[20],typeC[20];
3626 sscanf(line,
"%4s%d %ld %d %s %f %f %f %f %s",bank,&numBank,&
mNbPoint,&nbRecords,
3627 binTypeC,&bcoeff[0],&bcoeff[1],&bcoeff[2],&bcoeff[3],typeC);
3631 if(binType==
"CONST") binType=
"CONS";
3632 if((type!=
"ALT")&&(type!=
"ESD")) type=
"STD";
3634 cout<<
"BANK #"<<numBank<<endl;
3635 cout<<
"Number of data points:"<<
mNbPoint<<endl;
3636 cout<<
"Number of records:"<<nbRecords<<endl;
3637 cout<<
"BinType:"<<binType<<endl;
3638 cout<<
"BCoeff[1-4]:"<<bcoeff[0]<<
","<<bcoeff[1]<<
","<<bcoeff[2]<<
","<<bcoeff[3]<<endl;
3639 cout<<
"Type:"<<type<<endl;
3644 bool importOK=
false;
3645 if((binType==
"CONS") && (type==
"ESD"))
3649 unsigned long point=0;
3652 for(
long i=0;i<nbRecords;i++)
3656 while(isprint(fin.peek())==
false)
3658 if(fin.eof())
break;
3661 for(
unsigned int j=0;j<5;j++)
3667 substr=string(line).substr(j*16+0 ,8);
3668 istringstream(substr) >> iobs;
3669 substr=string(line).substr(j*16+8 ,8);
3670 istringstream(substr) >> isig;
3680 if((binType==
"CONS") && (type==
"STD"))
3683 unsigned long point=0;
3687 for(
long i=0;i<nbRecords;i++)
3691 while(isprint(fin.peek())==
false)
3693 if(fin.eof())
break;
3696 for(
unsigned int j=0;j<10;j++)
3707 substr=string(line).substr(j*8+0 ,2);
3708 if(substr==
" ") nc=1;
3709 else sscanf(substr.c_str(),
"%d",&nc);
3710 substr=string(line).substr(j*8+2 ,6);
3711 istringstream(substr) >> iobs;
3720 if((binType==
"RALF") && (type==
"ALT"))
3726 unsigned long point=0;
3727 REAL x,iobs,iobssigma;
3729 for(
long i=0;i<nbRecords;i++)
3733 while(isprint(fin.peek())==
false)
3735 if(fin.eof())
break;
3738 for(
unsigned int j=0;j<4;j++)
3744 substr=string(line).substr(j*20+0 ,8);
3745 istringstream(substr) >> x;
3746 substr=string(line).substr(j*20+8 ,7);
3747 istringstream(substr) >> iobs;
3748 substr=string(line).substr(j*20+15,5);
3749 istringstream(substr) >> iobssigma;
3759 for(
unsigned long i=0;i<(
mNbPoint/2);i++)
3783 this type of format is not handled yet (send an example file to the Fox author)!:"+filename);
3790 if(this->
GetRadiation().GetWavelengthType()==WAVELENGTH_TOF)
3793 sprintf(buf,
"Imported powder pattern: %d points, tof=%7.3f us-> %7.3f us",
3796 (*fpObjCrystInformUser)((string)buf);
3801 sprintf(buf,
"Imported powder pattern: %d points, 2theta=%7.3f -> %7.3f, step=%6.3f",
3805 (*fpObjCrystInformUser)((string)buf);
3807 VFN_DEBUG_EXIT(
"PowderPattern::ImportPowderPatternGSAS():file:"<<filename,5)
3812 VFN_DEBUG_ENTRY(
"PowderPattern::ImportPowderPatternCIF():file:",5)
3813 for(map<string,CIFData>::const_iterator pos=cif.
mvData.begin();pos!=cif.
mvData.end();++pos)
3814 if(pos->second.mPowderPatternObs.size()>10)
3816 mNbPoint=pos->second.mPowderPatternObs.size();
3821 if(pos->second.mDataType==WAVELENGTH_TOF)
3828 for(
unsigned long i=0;i<
mNbPoint;++i)
3831 mX(i)=pos->second.mPowderPatternX[i];
3837 VFN_DEBUG_EXIT(
"PowderPattern::ImportPowderPatternCIF():file:",5)
3843 VFN_DEBUG_MESSAGE(
"PowderPattern::SetPowderPatternObs()",5)
3844 if((
unsigned long)obs.numElements() !=
mNbPoint)
3847 supplied vector of observed intensities does not have the expected number of points!"));
3856 mClockIntegratedFactorsPrep.Reset();
3859 sprintf(buf,
"Changed powder pattern: %d points, 2theta=%7.3f -> %7.3f, step=%6.3f",
3862 (*fpObjCrystInformUser)((string)buf);
3868 VFN_DEBUG_MESSAGE(
"PowderPattern::SetPowderPatternObsSigma()",5)
3869 if((
unsigned long)sigma.numElements() !=
mNbPoint)
3872 supplied vector of sigma does not have the expected number of points!"));
3878 mClockIntegratedFactorsPrep.Reset();
3882 VFN_DEBUG_MESSAGE(
"PowderPattern::SavePowderPattern",5)
3884 ofstream out(filename.c_str());
3885 CrystVector_REAL ttheta;
3887 if(this->
GetRadiation().GetWavelengthType()!=WAVELENGTH_TOF) ttheta *= RAD2DEG;
3889 CrystVector_REAL diff;
3892 out <<
"# 2Theta/TOF Iobs ICalc Iobs-Icalc Weight Comp0" << endl;
3893 out << FormatVertVector<REAL>(ttheta,
3897 mPowderPatternComponentRegistry.GetObj(0).mPowderPatternCalc,16,8);
3899 VFN_DEBUG_MESSAGE(
"DiffractionDataPowder::SavePowderPattern:End",3)
3904 VFN_DEBUG_MESSAGE(
"DiffractionDataPowder::PrintObsCalcData()",5);
3905 CrystVector_REAL ttheta;
3907 if(this->
GetRadiation().GetWavelengthType()!=WAVELENGTH_TOF) ttheta *= RAD2DEG;
3908 os <<
"PowderPattern : " <<
mName <<endl;
3909 os <<
" 2Theta/TOF Obs Sigma Calc Weight" <<endl;
3923 TAU_PROFILE(
"PowderPattern::GetR()",
"void ()",TAU_DEFAULT);
3928 unsigned long maxPoints=mNbPointUsed;
3929 if( (
true==mStatisticsExcludeBackground)
3932 const REAL *p1, *p2, *p3;
3939 VFN_DEBUG_MESSAGE(
"PowderPattern::GetR():Exclude Backgd",4);
3940 for(
unsigned long i=0;i<maxPoints;i++)
3942 tmp1 += ((*p1)-(*p2)) * ((*p1)-(*p2));
3943 tmp2 += ((*p2)-(*p3)) * ((*p2)-(*p3));
3949 VFN_DEBUG_MESSAGE(
"PowderPattern::GetR():Exclude Backgd,Exclude regions",4);
3950 unsigned long min,max;
3952 for(
int j=0;j<nbExclude;j++)
3956 if(min>maxPoints)
break;
3957 if(max>maxPoints)max=maxPoints;
3960 tmp1 += ((*p1)-(*p2)) * ((*p1)-(*p2));
3961 tmp2 += ((*p2)-(*p3)) * ((*p2)-(*p3));
3969 for(;i<maxPoints;i++)
3971 tmp1 += ((*p1)-(*p2)) * ((*p1)-(*p2));
3972 tmp2 += ((*p2)-(*p3)) * ((*p2)-(*p3));
3980 const REAL *p1, *p2;
3986 VFN_DEBUG_MESSAGE(
"PowderPattern::GetR()",4);
3987 for(
unsigned long i=0;i<maxPoints;i++)
3989 tmp1 += ((*p1)-(*p2))*((*p1)-(*p2));
3990 tmp2 += (*p2) * (*p2);
3997 VFN_DEBUG_MESSAGE(
"PowderPattern::GetR(),Exclude regions",4);
3998 unsigned long min,max;
4000 for(
int j=0;j<nbExclude;j++)
4004 if(min>maxPoints)
break;
4005 if(max>maxPoints)max=maxPoints;
4008 tmp1 += ((*p1)-(*p2))*((*p1)-(*p2));
4009 tmp2 += (*p2) * (*p2);
4016 for(;i<maxPoints;i++)
4018 tmp1 += ((*p1)-(*p2))*((*p1)-(*p2));
4019 tmp2 += (*p2) * (*p2);
4025 VFN_DEBUG_MESSAGE(
"PowderPattern::GetR()="<<sqrt(tmp1/tmp2),4);
4029 return sqrt(tmp1/tmp2);
4032 REAL PowderPattern::GetIntegratedR()
const
4041 VFN_DEBUG_ENTRY(
"PowderPattern::GetIntegratedR()",4);
4042 TAU_PROFILE(
"PowderPattern::GetIntegratedR()",
"void ()",TAU_DEFAULT);
4046 const long numInterval=mIntegratedPatternMin.numElements();
4047 if( (
true==mStatisticsExcludeBackground)
4050 const REAL *p1, *p2, *p3;
4051 CrystVector_REAL integratedCalc(numInterval);
4053 CrystVector_REAL backgdCalc(numInterval);
4055 REAL *pp1=integratedCalc.data();
4056 REAL *pp2=backgdCalc.data();
4057 for(
int i=0;i<numInterval;i++)
4059 const long max=mIntegratedPatternMax(i);
4061 for(
int j=mIntegratedPatternMin(i);j<=max;j++) *pp1 += *p1++;
4064 for(
int j=mIntegratedPatternMin(i);j<=max;j++) *pp2 += *p1++;
4068 p1=integratedCalc.data();
4069 p2=mIntegratedObs.data();
4070 p3=backgdCalc.data();
4071 VFN_DEBUG_MESSAGE(
"PowderPattern::GetIntegratedR():Exclude Backgd",2);
4072 for(
long i=0;i<numInterval;i++)
4074 tmp1 += ((*p1)-(*p2)) * ((*p1)-(*p2));
4075 tmp2 += ((*p2)-(*p3)) * ((*p2)-(*p3));
4081 const REAL *p1, *p2;
4082 CrystVector_REAL integratedCalc(numInterval);
4084 REAL *pp1=integratedCalc.data();
4085 for(
int i=0;i<numInterval;i++)
4087 const long max=mIntegratedPatternMax(i);
4089 for(
int j=mIntegratedPatternMin(i);j<=max;j++) *pp1 += *p1++;
4092 p1=integratedCalc.data();
4093 p2=mIntegratedObs.data();
4094 VFN_DEBUG_MESSAGE(
"PowderPattern::GetIntegratedR()",2);
4095 for(
long i=0;i<numInterval;i++)
4097 tmp1 += ((*p1)-(*p2))*((*p1)-(*p2));
4098 tmp2 += (*p2) * (*p2);
4104 VFN_DEBUG_EXIT(
"PowderPattern::GetIntegratedR()="<<sqrt(tmp1/tmp2),4);
4105 return sqrt(tmp1/tmp2);
4116 TAU_PROFILE(
"PowderPattern::GetRw()",
"void ()",TAU_DEFAULT);
4117 VFN_DEBUG_MESSAGE(
"PowderPattern::GetRw()",3);
4126 unsigned long maxPoints=mNbPointUsed;
4128 if( (
true==mStatisticsExcludeBackground)
4131 VFN_DEBUG_MESSAGE(
"PowderPattern::GetRw():Exclude Backgd",3);
4132 const REAL *p1, *p2, *p3, *p4;
4140 for(
unsigned long i=0;i<maxPoints;i++)
4142 tmp1 += *p4 * ((*p1)-(*p2)) * ((*p1)-(*p2));
4143 tmp2 += *p4++ * ((*p2)-(*p3)) * ((*p2)-(*p3));
4149 unsigned long min,max;
4151 for(
int j=0;j<nbExclude;j++)
4155 if(min>maxPoints)
break;
4156 if(max>maxPoints)max=maxPoints;
4159 tmp1 += *p4 * ((*p1)-(*p2)) * ((*p1)-(*p2));
4160 tmp2 += *p4++ * ((*p2)-(*p3)) * ((*p2)-(*p3));
4169 for(;i<maxPoints;i++)
4171 tmp1 += *p4 * ((*p1)-(*p2)) * ((*p1)-(*p2));
4172 tmp2 += *p4++ * ((*p2)-(*p3)) * ((*p2)-(*p3));
4180 VFN_DEBUG_MESSAGE(
"PowderPattern::GetRw()",3);
4181 const REAL *p1, *p2, *p4;
4188 for(
unsigned long i=0;i<maxPoints;i++)
4190 tmp1 += *p4 * ((*p1)-(*p2))*((*p1)-(*p2));
4191 tmp2 += *p4++ * (*p2) * (*p2);
4197 unsigned long min,max;
4199 for(
int j=0;j<nbExclude;j++)
4203 if(min>maxPoints)
break;
4204 if(max>maxPoints)max=maxPoints;
4207 tmp1 += *p4 * ((*p1)-(*p2))*((*p1)-(*p2));
4208 tmp2 += *p4++ * (*p2) * (*p2);
4216 for(;i<maxPoints;i++)
4218 tmp1 += *p4 * ((*p1)-(*p2))*((*p1)-(*p2));
4219 tmp2 += *p4++ * (*p2) * (*p2);
4224 VFN_DEBUG_MESSAGE(
"PowderPattern::GetRw()="<<sqrt(tmp1/tmp2),3);
4225 return sqrt(tmp1/tmp2);
4227 REAL PowderPattern::GetIntegratedRw()
const
4236 TAU_PROFILE(
"PowderPattern::GetIntegratedRw()",
"void ()",TAU_DEFAULT);
4240 const long numInterval=mIntegratedPatternMin.numElements();
4241 if( (
true==mStatisticsExcludeBackground)
4244 const REAL *p1, *p2, *p3, *p4;
4245 CrystVector_REAL integratedCalc(numInterval);
4247 CrystVector_REAL backgdCalc(numInterval);
4249 REAL *pp1=integratedCalc.data();
4250 REAL *pp2=backgdCalc.data();
4251 for(
int i=0;i<numInterval;i++)
4253 const long max=mIntegratedPatternMax(i);
4255 for(
int j=mIntegratedPatternMin(i);j<=max;j++) *pp1 += *p1++;
4258 for(
int j=mIntegratedPatternMin(i);j<=max;j++) *pp2 += *p1++;
4262 p1=integratedCalc.data();
4263 p2=mIntegratedObs.data();
4264 p3=backgdCalc.data();
4265 if(mIntegratedWeight.numElements()==0) p4=mIntegratedWeightObs.data();
4266 else p4=mIntegratedWeight.data();
4267 VFN_DEBUG_MESSAGE(
"PowderPattern::GetIntegratedRw():Exclude Backgd",4);
4268 for(
long i=0;i<numInterval;i++)
4270 tmp1 += *p4 * ((*p1)-(*p2)) * ((*p1)-(*p2));
4271 tmp2 += *p4++ * ((*p2)-(*p3)) * ((*p2)-(*p3));
4279 const REAL *p1, *p2, *p4;
4280 CrystVector_REAL integratedCalc(numInterval);
4282 REAL *pp1=integratedCalc.data();
4283 for(
int i=0;i<numInterval;i++)
4285 const long max=mIntegratedPatternMax(i);
4287 for(
int j=mIntegratedPatternMin(i);j<=max;j++) *pp1 += *p1++;
4290 p1=integratedCalc.data();
4291 p2=mIntegratedObs.data();
4292 if(mIntegratedWeight.numElements()==0) p4=mIntegratedWeightObs.data();
4293 else p4=mIntegratedWeight.data();
4294 VFN_DEBUG_MESSAGE(
"PowderPattern::GetIntegratedRw()",4);
4295 for(
long i=0;i<numInterval;i++)
4297 tmp1 += *p4 * ((*p1)-(*p2))*((*p1)-(*p2));
4298 tmp2 += *p4++ * (*p2) * (*p2);
4304 VFN_DEBUG_MESSAGE(
"PowderPattern::GetIntegratedRw()="<<sqrt(tmp1/tmp2),4);
4308 return sqrt(tmp1/tmp2);
4329 TAU_PROFILE(
"PowderPattern::GetChi2()",
"void ()",TAU_DEFAULT);
4331 VFN_DEBUG_ENTRY(
"PowderPattern::GetChi2()",3);
4333 const unsigned long maxPoints=mNbPointUsed;
4337 VFN_DEBUG_MESSAGE(
"PowderPattern::GetChi2()Integrated profiles",3);
4338 const REAL * RESTRICT p1, * RESTRICT p2, * RESTRICT p3;
4345 for(
unsigned long i=0;i<maxPoints;i++)
4347 mChi2 += *p3 * ((*p1)-(*p2))*((*p1)-(*p2));
4349 else mChi2LikeNorm -= log(*p3++);
4355 unsigned long min,max;
4357 for(
int j=0;j<nbExclude;j++)
4361 if(min>maxPoints)
break;
4362 if(max>maxPoints)max=maxPoints;
4365 mChi2 += *p3 * ((*p1)-(*p2))*((*p1)-(*p2));
4367 else mChi2LikeNorm -= log(*p3++);
4375 for(;i<maxPoints;i++)
4377 mChi2 += *p3 * ((*p1)-(*p2))*((*p1)-(*p2));
4379 else mChi2LikeNorm -= log(*p3++);
4384 VFN_DEBUG_MESSAGE(
"Chi^2="<<mChi2<<
", log(norm)="<<mChi2LikeNorm,3)
4386 VFN_DEBUG_EXIT(
"PowderPattern::GetChi2()="<<mChi2,3);
4396 return mIntegratedChi2;
4399 if(mClockIntegratedChi2>
mClockMaster)
return mIntegratedChi2;
4407 this->FitScaleFactorForIntegratedRw();
4409 TAU_PROFILE(
"PowderPattern::GetIntegratedChi2()",
"void ()",TAU_DEFAULT);
4411 VFN_DEBUG_ENTRY(
"PowderPattern::GetIntegratedChi2()",3);
4414 mIntegratedChi2LikeNorm=0.;
4415 VFN_DEBUG_MESSAGE(
"PowderPattern::GetIntegratedChi2()",3);
4416 const REAL * RESTRICT p1, * RESTRICT p2, * RESTRICT p3;
4418 p2=mIntegratedObs.data();
4419 if(mIntegratedWeight.numElements()==0) p3=mIntegratedWeightObs.data();
4420 else p3=mIntegratedWeight.data();
4421 double weightProd=1.;
4422 VFN_DEBUG_MESSAGE(
"PowderPattern::GetIntegratedIntegratedRw()",4);
4423 for(
unsigned long i=0;i<mNbIntegrationUsed;)
4427 for(
unsigned long j=0;j<32;++j)
4429 mIntegratedChi2 += *p3 * ((*p1)-(*p2))*((*p1)-(*p2));
4430 if(*p3>0) weightProd *= *p3;
4432 if(++i == mNbIntegrationUsed)
break;
4434 mIntegratedChi2LikeNorm -= log(weightProd);
4437 mIntegratedChi2LikeNorm/=2;
4438 VFN_DEBUG_MESSAGE(
"Chi^2="<<mIntegratedChi2<<
", log(norm)="<<mIntegratedChi2LikeNorm,3)
4439 mClockIntegratedChi2.Click();
4440 VFN_DEBUG_EXIT(
"PowderPattern::GetChi2()="<<mIntegratedChi2,3);
4441 return mIntegratedChi2;
4458 TAU_PROFILE(
"PowderPattern::FitScaleFactorForR()",
"void ()",TAU_DEFAULT);
4459 VFN_DEBUG_ENTRY(
"PowderPattern::FitScaleFactorForR()",3);
4461 mScalableComponentIndex.resize(mPowderPatternComponentRegistry.GetNb());
4463 for(
int i=0;i<mPowderPatternComponentRegistry.GetNb();i++)
4465 if(mPowderPatternComponentRegistry.GetObj(i).IsScalable())
4466 mScalableComponentIndex(nbScale++)=i;
4468 VFN_DEBUG_MESSAGE(
"-> Number of Scale Factors:"<<nbScale<<
":Index:"<<endl<<mScalableComponentIndex,3);
4471 VFN_DEBUG_EXIT(
"PowderPattern::FitScaleFactorForR(): No scalable component!",3);
4474 mScalableComponentIndex.resizeAndPreserve(nbScale);
4476 mFitScaleFactorM.resize(nbScale,nbScale);
4477 mFitScaleFactorB.resize(nbScale,1);
4478 mFitScaleFactorX.resize(nbScale,1);
4483 for(
int i=0;i<nbScale;i++)
4485 for(
int j=i;j<nbScale;j++)
4489 const REAL *p1=mPowderPatternComponentRegistry.GetObj(mScalableComponentIndex(i))
4490 .mPowderPatternCalc.data();
4491 const REAL *p2=mPowderPatternComponentRegistry.GetObj(mScalableComponentIndex(j))
4492 .mPowderPatternCalc.data();
4494 for(
unsigned long k=0;k<mNbPointUsed;k++) m += *p1++ * *p2++;
4495 mFitScaleFactorM(i,j)=m;
4496 mFitScaleFactorM(j,i)=m;
4499 for(
int i=0;i<nbScale;i++)
4502 const REAL *p2=mPowderPatternComponentRegistry.GetObj(mScalableComponentIndex(i))
4503 .mPowderPatternCalc.data();
4506 for(
unsigned long k=0;k<mNbPointUsed;k++) b += *p1++ * *p2++;
4510 for(
unsigned long k=0;k<mNbPointUsed;k++) b += (*p1++ - *p3++) * *p2++;
4512 mFitScaleFactorB(i,0) =b;
4517 unsigned long min,max;
4518 for(
int i=0;i<nbScale;i++)
4520 for(
int j=i;j<nbScale;j++)
4524 const REAL *p1=mPowderPatternComponentRegistry.GetObj(mScalableComponentIndex(i))
4525 .mPowderPatternCalc.data();
4526 const REAL *p2=mPowderPatternComponentRegistry.GetObj(mScalableComponentIndex(j))
4527 .mPowderPatternCalc.data();
4530 for(
int k=0;k<nbExclude;k++)
4534 if(min>mNbPointUsed)
break;
4535 if(max>mNbPointUsed)max=mNbPointUsed;
4537 for(;l<min;l++) m += *p1++ * *p2++;
4542 for(;l<mNbPointUsed;l++) m += *p1++ * *p2++;
4543 mFitScaleFactorM(i,j)=m;
4544 mFitScaleFactorM(j,i)=m;
4547 for(
int i=0;i<nbScale;i++)
4550 const REAL *p2=mPowderPatternComponentRegistry
4551 .GetObj(mScalableComponentIndex(i))
4552 .mPowderPatternCalc.data();
4557 for(
int k=0;k<nbExclude;k++)
4561 if(min>mNbPointUsed)
break;
4562 if(max>mNbPointUsed)max=mNbPointUsed;
4564 for(;l<min;l++) b += *p1++ * *p2++;
4569 for(;l<mNbPointUsed;l++) b += *p1++ * *p2++;
4574 for(
int k=0;k<nbExclude;k++)
4578 if(min>mNbPointUsed)
break;
4579 if(max>mNbPointUsed)max=mNbPointUsed;
4581 for(;l<min;l++) b += (*p1++ - *p3++) * *p2++;
4586 for(;l<mNbPointUsed;l++) b += (*p1++ - *p3++) * *p2++;
4588 mFitScaleFactorB(i,0) =b;
4591 if(1==nbScale) mFitScaleFactorX=mFitScaleFactorB(0)/mFitScaleFactorM(0);
4593 mFitScaleFactorX=product(InvertMatrix(mFitScaleFactorM),mFitScaleFactorB);
4594 VFN_DEBUG_MESSAGE(
"B, M, X"<<endl<<mFitScaleFactorB<<endl<<mFitScaleFactorM<<endl<<mFitScaleFactorX,2)
4595 for(
int i=0;i<nbScale;i++)
4597 const REAL * p1=mPowderPatternComponentRegistry.GetObj(mScalableComponentIndex(i))
4598 .mPowderPatternCalc.data();
4600 const REAL s = mFitScaleFactorX(i)
4601 -mScaleFactor(mScalableComponentIndex(i));
4604 (*fpObjCrystInformUser)(
"Warning:FitScaleFactorForR: working around NaN scale factor...");
4607 for(
unsigned long j=0;j<mNbPointUsed;j++) *p0++ += s * *p1++;
4608 VFN_DEBUG_MESSAGE(
"-> Old:"<<mScaleFactor(mScalableComponentIndex(i)) <<
" Change:"<<mFitScaleFactorX(i),2);
4609 mScaleFactor(mScalableComponentIndex(i)) = mFitScaleFactorX(i);
4613 VFN_DEBUG_EXIT(
"PowderPattern::FitScaleFactorForR():End",3);
4616 void PowderPattern::FitScaleFactorForIntegratedR()
const
4625 VFN_DEBUG_ENTRY(
"PowderPattern::FitScaleFactorForIntegratedR()",3);
4626 TAU_PROFILE(
"PowderPattern::FitScaleFactorForIntegratedR()",
"void ()",TAU_DEFAULT);
4628 mScalableComponentIndex.resize(mPowderPatternComponentRegistry.GetNb());
4630 for(
int i=0;i<mPowderPatternComponentRegistry.GetNb();i++)
4632 if(mPowderPatternComponentRegistry.GetObj(i).IsScalable())
4633 mScalableComponentIndex(nbScale++)=i;
4635 VFN_DEBUG_MESSAGE(
"-> Number of Scale Factors:"<<nbScale<<
":Index:"<<endl<<mScalableComponentIndex,2);
4638 VFN_DEBUG_EXIT(
"PowderPattern::FitScaleFactorForIntegratedR(): No scalable component!",3);
4641 mScalableComponentIndex.resizeAndPreserve(nbScale);
4643 mFitScaleFactorM.resize(nbScale,nbScale);
4644 mFitScaleFactorB.resize(nbScale,1);
4645 mFitScaleFactorX.resize(nbScale,1);
4647 VFN_DEBUG_MESSAGE(
"PowderPattern::FitScaleFactorForIntegratedR():1",2);
4648 const long numInterval=mIntegratedPatternMin.numElements();
4649 CrystVector_REAL *integratedCalc=
new CrystVector_REAL[nbScale];
4650 for(
int i=0;i<nbScale;i++)
4652 integratedCalc[i].resize(numInterval);
4656 const REAL *p1=mPowderPatternComponentRegistry.GetObj(mScalableComponentIndex(i))
4657 .mPowderPatternCalc.data();
4659 REAL *p2=integratedCalc[i].data();
4660 for(
int j=0;j<numInterval;j++)
4662 const long max=mIntegratedPatternMax(j);
4663 p1=mPowderPatternComponentRegistry.GetObj(mScalableComponentIndex(i))
4664 .mPowderPatternCalc.data()+mIntegratedPatternMin(j);
4666 for(
int k=mIntegratedPatternMin(j);k<=max;k++) *p2 += *p1++;
4673 VFN_DEBUG_MESSAGE(
"PowderPattern::FitScaleFactorForIntegratedR():2",2);
4674 CrystVector_REAL backdIntegrated(numInterval);
4678 REAL *p2=backdIntegrated.data();
4679 for(
int j=0;j<numInterval;j++)
4681 const long max=mIntegratedPatternMax(j);
4684 for(
int k=mIntegratedPatternMin(j);k<=max;k++) *p2 += *p1++;
4696 VFN_DEBUG_MESSAGE(
"PowderPattern::FitScaleFactorForIntegratedR():3",2);
4697 for(
int i=0;i<nbScale;i++)
4699 for(
int j=i;j<nbScale;j++)
4701 const REAL *p1=integratedCalc[i].data();
4702 const REAL *p2=integratedCalc[j].data();
4704 for(
long k=0;k<numInterval;k++)
4709 mFitScaleFactorM(i,j)=m;
4710 mFitScaleFactorM(j,i)=m;
4713 VFN_DEBUG_MESSAGE(
"PowderPattern::FitScaleFactorForIntegratedR():4",2);
4714 for(
int i=0;i<nbScale;i++)
4716 const REAL *p1=mIntegratedObs.data();
4717 const REAL *p2=integratedCalc[i].data();
4720 for(
long k=0;k<numInterval;k++)
4727 const REAL *p3=backdIntegrated.data();
4728 for(
long k=0;k<numInterval;k++)
4733 b += (*p1++ - *p3++) * *p2++;
4736 mFitScaleFactorB(i,0) =b;
4738 VFN_DEBUG_MESSAGE(
"PowderPattern::FitScaleFactorForIntegratedR():5",2);
4740 if(1==nbScale) mFitScaleFactorX=mFitScaleFactorB(0)/mFitScaleFactorM(0);
4742 mFitScaleFactorX=product(InvertMatrix(mFitScaleFactorM),mFitScaleFactorB);
4743 VFN_DEBUG_MESSAGE(
"B, M, X"<<endl<<mFitScaleFactorB<<endl<<mFitScaleFactorM<<endl<<mFitScaleFactorX,3)
4744 for(
int i=0;i<nbScale;i++)
4746 const REAL * p1=mPowderPatternComponentRegistry.GetObj(mScalableComponentIndex(i))
4747 .mPowderPatternCalc.data();
4749 const REAL s = mFitScaleFactorX(i)
4750 -mScaleFactor(mScalableComponentIndex(i));
4753 (*fpObjCrystInformUser)(
"Warning:FitScaleFactorForIntegratedR: working around NaN scale factor...");
4756 for(
unsigned long j=0;j<mNbPointUsed;j++) *p0++ += s * *p1++;
4757 VFN_DEBUG_MESSAGE(
"-> Old:"<<mScaleFactor(mScalableComponentIndex(i)) <<
" New:"<<mFitScaleFactorX(i),3);
4758 mScaleFactor(mScalableComponentIndex(i)) = mFitScaleFactorX(i);
4762 delete[] integratedCalc;
4763 VFN_DEBUG_EXIT(
"PowderPattern::FitScaleFactorForIntegratedR():End",3);
4773 TAU_PROFILE(
"PowderPattern::FitScaleFactorForRw()",
"void ()",TAU_DEFAULT);
4774 VFN_DEBUG_ENTRY(
"PowderPattern::FitScaleFactorForRw()",3);
4778 mScalableComponentIndex.resize(mPowderPatternComponentRegistry.GetNb());
4780 for(
int i=0;i<mPowderPatternComponentRegistry.GetNb();i++)
4782 if(mPowderPatternComponentRegistry.GetObj(i).IsScalable())
4783 mScalableComponentIndex(nbScale++)=i;
4785 VFN_DEBUG_MESSAGE(
"-> Number of Scale Factors:"<<nbScale<<
":Index:"<<endl<<mScalableComponentIndex,2);
4788 VFN_DEBUG_EXIT(
"PowderPattern::FitScaleFactorForRw(): No scalable component!",3);
4791 mScalableComponentIndex.resizeAndPreserve(nbScale);
4793 mFitScaleFactorM.resize(nbScale,nbScale);
4794 mFitScaleFactorB.resize(nbScale,1);
4795 mFitScaleFactorX.resize(nbScale,1);
4800 for(
int i=0;i<nbScale;i++)
4802 for(
int j=i;j<nbScale;j++)
4806 const REAL *p1=mPowderPatternComponentRegistry.GetObj(mScalableComponentIndex(i))
4807 .mPowderPatternCalc.data();
4808 const REAL *p2=mPowderPatternComponentRegistry.GetObj(mScalableComponentIndex(j))
4809 .mPowderPatternCalc.data();
4812 for(
unsigned long k=0;k<mNbPointUsed;k++) m += *p1++ * *p2++ * *p3++;
4813 mFitScaleFactorM(i,j)=m;
4814 mFitScaleFactorM(j,i)=m;
4817 for(
int i=0;i<nbScale;i++)
4820 const REAL *p2=mPowderPatternComponentRegistry.GetObj(mScalableComponentIndex(i))
4821 .mPowderPatternCalc.data();
4825 for(
unsigned long k=0;k<mNbPointUsed;k++) b += *p1++ * *p2++ * *p3++;
4829 for(
unsigned long k=0;k<mNbPointUsed;k++)
4830 b += (*p1++ - *p4++) * *p2++ * *p3++;
4832 mFitScaleFactorB(i,0) =b;
4837 unsigned long min,max;
4838 for(
int i=0;i<nbScale;i++)
4840 for(
int j=i;j<nbScale;j++)
4844 const REAL *p1=mPowderPatternComponentRegistry.GetObj(mScalableComponentIndex(i))
4845 .mPowderPatternCalc.data();
4846 const REAL *p2=mPowderPatternComponentRegistry.GetObj(mScalableComponentIndex(j))
4847 .mPowderPatternCalc.data();
4851 for(
int k=0;k<nbExclude;k++)
4855 if(min>mNbPointUsed)
break;
4856 if(max>mNbPointUsed)max=mNbPointUsed;
4858 for(;l<min;l++) m += *p1++ * *p2++ * *p3++;
4864 for(;l<mNbPointUsed;l++) m += *p1++ * *p2++ * *p3++;
4865 mFitScaleFactorM(i,j)=m;
4866 mFitScaleFactorM(j,i)=m;
4869 for(
int i=0;i<nbScale;i++)
4872 const REAL *p2=mPowderPatternComponentRegistry
4873 .GetObj(mScalableComponentIndex(i))
4874 .mPowderPatternCalc.data();
4880 for(
int k=0;k<nbExclude;k++)
4884 if(min>mNbPointUsed)
break;
4885 if(max>mNbPointUsed)max=mNbPointUsed;
4887 for(;l<min;l++) b += *p1++ * *p2++ * *p3++;
4893 for(;l<mNbPointUsed;l++) b += *p1++ * *p2++ * *p3++;
4898 for(
int k=0;k<nbExclude;k++)
4902 if(min>mNbPointUsed)
break;
4903 if(max>mNbPointUsed)max=mNbPointUsed;
4905 for(;l<min;l++) b += (*p1++ - *p4++) * *p2++ * *p3++;
4911 for(;l<mNbPointUsed;l++) b += (*p1++ - *p4++) * *p2++ * *p3++;
4913 mFitScaleFactorB(i,0) =b;
4916 if(1==nbScale) mFitScaleFactorX=mFitScaleFactorB(0)/mFitScaleFactorM(0);
4918 mFitScaleFactorX=product(InvertMatrix(mFitScaleFactorM),mFitScaleFactorB);
4919 VFN_DEBUG_MESSAGE(
"B, M, X"<<endl<<mFitScaleFactorB<<endl<<mFitScaleFactorM<<endl<<mFitScaleFactorX,2)
4920 for(
int i=0;i<nbScale;i++)
4922 const REAL * p1=mPowderPatternComponentRegistry.GetObj(mScalableComponentIndex(i))
4923 .mPowderPatternCalc.data();
4925 const REAL s = mFitScaleFactorX(i)
4926 -mScaleFactor(mScalableComponentIndex(i));
4929 (*fpObjCrystInformUser)(
"Warning:FitScaleFactorForRw working around NaN scale factor...");
4932 for(
unsigned long j=0;j<mNbPointUsed;j++) *p0++ += s * *p1++;
4933 VFN_DEBUG_MESSAGE(
"-> Old:"<<mScaleFactor(mScalableComponentIndex(i)) <<
" Change:"<<mFitScaleFactorX(i),3);
4934 mScaleFactor(mScalableComponentIndex(i)) = mFitScaleFactorX(i);
4938 VFN_DEBUG_EXIT(
"PowderPattern::FitScaleFactorForRw():End",3);
4941 void PowderPattern::FitScaleFactorForIntegratedRw()
const
4950 VFN_DEBUG_ENTRY(
"PowderPattern::FitScaleFactorForIntegratedRw()",3);
4951 TAU_PROFILE(
"PowderPattern::FitScaleFactorForIntegratedRw()",
"void ()",TAU_DEFAULT);
4953 mScalableComponentIndex.resize(mPowderPatternComponentRegistry.GetNb());
4954 CrystVector_REAL vScalableVarianceIndex(mPowderPatternComponentRegistry.GetNb());
4957 for(
int i=0;i<mPowderPatternComponentRegistry.GetNb();i++)
4959 if(mPowderPatternComponentRegistry.GetObj(i).IsScalable())
4960 mScalableComponentIndex(nbScale++)=i;
4961 if(mPowderPatternComponentRegistry.GetObj(i).HasPowderPatternCalcVariance())
4963 if(0!=mPowderPatternComponentRegistry.GetObj(i)
4964 .GetPowderPatternIntegratedCalcVariance().first->numElements())
4965 vScalableVarianceIndex(nbVarCalc++)=i;
4968 VFN_DEBUG_MESSAGE(
"-> Number of Scale Factors:"<<nbScale<<
"("<<nbVarCalc
4969 <<
"with variance). Index:"<<endl<<mScalableComponentIndex,2);
4972 VFN_DEBUG_EXIT(
"PowderPattern::FitScaleFactorForIntegratedRw(): No scalable component!",3);
4976 unsigned int ctagain=0;
4977 RECALC_SCALE_FACTOR_VARIANCE_FitScaleFactorForIntegratedRw:
4981 double a=0.,b=0.,c=0.;
4984 const REAL * RESTRICT p1=mIntegratedObs.data();
4985 const REAL * RESTRICT p2=mPowderPatternComponentRegistry.GetObj(mScalableComponentIndex(0))
4986 .GetPowderPatternIntegratedCalc().first->data();
4987 const REAL * RESTRICT p1v=mIntegratedVarianceObs.data();
4988 const REAL * RESTRICT p2v=mPowderPatternComponentRegistry.GetObj(mScalableComponentIndex(0))
4989 .GetPowderPatternIntegratedCalcVariance().first->data();
4992 for(
unsigned long k=mNbIntegrationUsed;k>0;k--)
4994 a += *p2v * *p1 * *p2;
4995 b += *p2 * *p2 * *p1v - *p1 * *p1 * *p2v;
4996 c -= *p1 * *p2 * *p1v;
4997 p1++;p2++;p1v++;p2v++;
5003 for(
unsigned long k=mNbIntegrationUsed;k>0;k--)
5005 a += *p2v * (*p1 - *p3) * *p2;
5006 b += *p2 * *p2 * *p1v - (*p1 - *p3) * (*p1 - *p3) * *p2v;
5007 c -= (*p1 - *p3) * *p2 * *p1v;
5008 p1++;p2++;p1v++;p2v++;p3++;
5012 newscale=(REAL)((-b+sqrt(b*b-4*a*c))/(2*a));
5015 const REAL s = newscale-mScaleFactor(mScalableComponentIndex(0));
5016 const REAL s2 = newscale*newscale
5017 -mScaleFactor(mScalableComponentIndex(0))
5018 *mScaleFactor(mScalableComponentIndex(0));
5020 const REAL * RESTRICT p1=mPowderPatternComponentRegistry.GetObj(mScalableComponentIndex(0))
5021 .GetPowderPatternIntegratedCalc().first->data();
5023 const REAL * RESTRICT p1v=mPowderPatternComponentRegistry.GetObj(mScalableComponentIndex(0))
5024 .GetPowderPatternIntegratedCalcVariance().first->data();
5025 REAL * RESTRICT p0w = mIntegratedWeight.data();
5026 for(
unsigned long j=mNbIntegrationUsed;j>0;j--)
5029 *p0v += s2 * *p1v++;
5030 if(*p0v <=0) {*p0w++ =0;p0v++;}
else *p0w++ = 1. / *p0v++;
5033 mScaleFactor(mScalableComponentIndex(0)) = newscale;
5040 mScalableComponentIndex.resizeAndPreserve(nbScale);
5042 mFitScaleFactorM.resize(nbScale,nbScale);
5043 mFitScaleFactorB.resize(nbScale,1);
5044 mFitScaleFactorX.resize(nbScale,1);
5046 VFN_DEBUG_MESSAGE(
"PowderPattern::FitScaleFactorForIntegratedRw():1",2);
5047 vector<const CrystVector_REAL*> integratedCalc;
5048 for(
int i=0;i<nbScale;i++)
5050 integratedCalc.push_back(mPowderPatternComponentRegistry.
5051 GetObj(mScalableComponentIndex(i))
5052 .GetPowderPatternIntegratedCalc().first);
5054 VFN_DEBUG_MESSAGE(
"PowderPattern::FitScaleFactorForIntegratedRw():3",2);
5055 for(
int i=0;i<nbScale;i++)
5057 for(
int j=i;j<nbScale;j++)
5059 const REAL * RESTRICT p1=integratedCalc[i]->data();
5060 const REAL * RESTRICT p2=integratedCalc[j]->data();
5061 const REAL * RESTRICT p3;
5062 if(mIntegratedWeight.numElements()==0)
5064 p3=mIntegratedWeightObs.data();
5065 if(ctagain>5) VFN_DEBUG_MESSAGE(
"ctagain="<<ctagain<<
", using mIntegratedWeightObs",5);
5069 p3=mIntegratedWeight.data();
5070 if(ctagain>5) VFN_DEBUG_MESSAGE(
"ctagain="<<ctagain<<
", using mIntegratedWeight",5);
5074 for(
unsigned long k=mNbIntegrationUsed;k>0;k--)
5075 {m += *p1 * *p1 * *p3++; p1++;}
5077 for(
unsigned long k=mNbIntegrationUsed;k>0;k--)
5078 m += *p1++ * *p2++ * *p3++;
5079 mFitScaleFactorM(i,j)=m;
5080 mFitScaleFactorM(j,i)=m;
5083 VFN_DEBUG_MESSAGE(
"PowderPattern::FitScaleFactorForIntegratedRw():4",2);
5084 for(
int i=0;i<nbScale;i++)
5086 const REAL * RESTRICT p1=mIntegratedObs.data();
5087 const REAL * RESTRICT p2=integratedCalc[i]->data();
5088 const REAL * RESTRICT p4;
5089 if(mIntegratedWeight.numElements()==0) p4=mIntegratedWeightObs.data();
5090 else p4=mIntegratedWeight.data();
5094 for(
unsigned long k=mNbIntegrationUsed;k>0;k--)
5096 b += *p1++ * *p2++ * *p4++;
5103 for(
unsigned long k=mNbIntegrationUsed;k>0;k--)
5108 b += (*p1++ - *p3++) * *p2++ * *p4++;
5111 mFitScaleFactorB(i,0) =b;
5113 VFN_DEBUG_MESSAGE(
"PowderPattern::FitScaleFactorForIntegratedRw():5",2);
5115 if(1==nbScale) mFitScaleFactorX=mFitScaleFactorB(0)/mFitScaleFactorM(0);
5117 mFitScaleFactorX=product(InvertMatrix(mFitScaleFactorM),mFitScaleFactorB);
5118 VFN_DEBUG_MESSAGE(
"B, M, X"<<endl<<mFitScaleFactorB<<endl<<mFitScaleFactorM<<endl<<mFitScaleFactorX,3)
5123 for(
int i=0;i<nbScale;i++)
5125 if(mPowderPatternComponentRegistry.GetObj(mScalableComponentIndex(i)).HasPowderPatternCalcVariance())
5127 if(0!=mPowderPatternComponentRegistry.GetObj(mScalableComponentIndex(i))
5128 .GetPowderPatternIntegratedCalcVariance().first->numElements())
5130 const REAL * RESTRICT p1=mPowderPatternComponentRegistry.GetObj(mScalableComponentIndex(i))
5131 .GetPowderPatternIntegratedCalcVariance().first->data();
5135 const REAL s2 = mFitScaleFactorX(i)*mFitScaleFactorX(i)
5136 -mScaleFactor(mScalableComponentIndex(i))
5137 *mScaleFactor(mScalableComponentIndex(i));
5138 for(
unsigned long j=mNbIntegrationUsed;j>0;j--) *p0++ += s2 * *p1++;
5144 REAL * RESTRICT p0 = mIntegratedWeight.data();
5146 for(
unsigned long j=mNbIntegrationUsed;j>0;j--)
5147 if(*p1 <=0) {*p0++ =0;p1++;}
5148 else *p0++ = 1. / *p1++;
5152 for(
int i=0;i<nbScale;i++)
5154 const REAL * RESTRICT p1=mPowderPatternComponentRegistry.GetObj(mScalableComponentIndex(i))
5155 .GetPowderPatternIntegratedCalc().first->data();
5157 const REAL s = mFitScaleFactorX(i)
5158 -mScaleFactor(mScalableComponentIndex(i));
5161 (*fpObjCrystInformUser)(
"Warning:FitScaleFactorForIntegratedRw: working around NaN scale factor...");
5166 if(abs(s/mFitScaleFactorX(i))>0.001)
5172 if((!again)&&(ctagain>0))
5174 VFN_DEBUG_MESSAGE(
"log(scale) :"<<log(mScaleFactor(mScalableComponentIndex(i)))
5175 <<
"->"<<log(mFitScaleFactorX(i))<<
" Again="<<ctagain,5);
5178 for(
unsigned long j=mNbIntegrationUsed;j>0;j--) *p0++ += s * *p1++;
5181 VFN_DEBUG_MESSAGE(
"->log(scale) Old :"<<log(mScaleFactor(mScalableComponentIndex(i))) <<
" New:"<<log(mFitScaleFactorX(i)),10);
5183 mScaleFactor(mScalableComponentIndex(i)) = mFitScaleFactorX(i);
5189 if(again && (ctagain<20))
5194 VFN_DEBUG_MESSAGE(
"PowderPattern::FitScaleFactorForIntegratedRw(), scaling again #"<<ctagain,10)
5196 goto RECALC_SCALE_FACTOR_VARIANCE_FitScaleFactorForIntegratedRw;
5198 VFN_DEBUG_EXIT(
"PowderPattern::FitScaleFactorForIntegratedRw():End",3);
5203 VFN_DEBUG_MESSAGE(
"PowderPattern::SetSigmaToSqrtIobs()",5);
5209 VFN_DEBUG_MESSAGE(
"PowderPattern::SetWeightToInvSigmaSq()",5);
5223 VFN_DEBUG_MESSAGE(
"PowderPattern::SetWeightToSinTheta()",5);
5229 const REAL minRelatIobs)
5231 VFN_DEBUG_MESSAGE(
"PowderPattern::SetWeightPolynomial()",5);
5244 const bool enableRestraints)
5247 if(0 == mOptProfileIntegration.GetChoice()) this->FitScaleFactorForIntegratedRw();
5269 VFN_DEBUG_MESSAGE(
"PowderPattern::AddExcludedRegion()",5)
5285 CrystVector_long subs;
5287 CrystVector_REAL tmp1,tmp2;
5296 VFN_DEBUG_MESSAGE(
"PowderPattern::Add2ThetaExcludedRegion():End",5)
5302 if(mOptProfileIntegration.GetChoice()==0) tmp+=mIntegratedChi2LikeNorm;
5303 else tmp+=mChi2LikeNorm;
5309 const CrystVector_REAL&
5312 TAU_PROFILE(
"PowderPattern::GetLSQCalc()",
"void ()",TAU_DEFAULT);
5331 const CrystVector_REAL&
5334 TAU_PROFILE(
"PowderPattern::GetLSQObs()",
"void ()",TAU_DEFAULT);
5353 const CrystVector_REAL&
5356 TAU_PROFILE(
"PowderPattern::GetLSQWeight()",
"void ()",TAU_DEFAULT);
5364 if(mIntegratedWeight.numElements()==0)
5381 TAU_PROFILE(
"PowderPattern::GetLSQ_FullDeriv()",
"void ()",TAU_DEFAULT);
5385 this->CalcPowderPatternIntegrated_FullDeriv(vPar);
5387 std::map<RefinablePar*, CrystVector_REAL> fullderiv_old;
5388 std::vector<const CrystVector_REAL*> v;
5391 cout<<
"PowderPattern::GetLSQ_FullDeriv(integrated):parameters:"<<endl;
5392 for(std::set<RefinablePar*>::iterator par=vPar.begin();par!=vPar.end();++par)
5394 v.push_back(&(mPowderPatternIntegrated_FullDeriv[*par]));
5395 fullderiv_old[*par]=this->
GetLSQDeriv(idx,*(*par));
5396 v.push_back(&(fullderiv_old[*par]));
5397 cout<<(*par)->GetName()<<
":"<<mPowderPatternIntegrated_FullDeriv[*par].size()<<
","<<fullderiv_old[*par].size()<<endl;
5400 cout<<
"PowderPattern::GetLSQ_FullDeriv(integrated):"<<endl<<FormatVertVector<REAL>(v,12,1,20)<<endl;
5403 return mPowderPatternIntegrated_FullDeriv;
5405 mPowderPattern_FullDeriv=this->GetPowderPattern_FullDeriv(vPar);
5406 return mPowderPattern_FullDeriv;
5411 VFN_DEBUG_MESSAGE(
"PowderPattern::Prepare()",5);
5412 for(
int i=0;i<mPowderPatternComponentRegistry.GetNb();i++)
5415 mPowderPatternComponentRegistry.GetObj(i).Prepare();
5419 CrystVector_uint & groupIndex,
5420 unsigned int &first)
const
5423 unsigned int scaleIndex=0;
5424 unsigned int thetaIndex=0;
5425 VFN_DEBUG_MESSAGE(
"PowderPattern::GetGeneGroup()",4)
5427 for(
long j=0;j<this->
GetNbPar();j++)
5432 if(scaleIndex==0) scaleIndex=first++;
5433 groupIndex(i)=scaleIndex;
5437 if(thetaIndex==0) thetaIndex=first++;
5438 groupIndex(i)=thetaIndex;
5446 for(
int i=0;i<mPowderPatternComponentRegistry.GetNb();i++)
5457 return mIntegratedPatternMin;
5463 return mIntegratedPatternMax;
5468 return mClockIntegratedFactorsPrep;
5474 if(this->
GetRadiation().GetWavelengthType()==WAVELENGTH_TOF)
5476 if(stol>0) x = 1.0/(2*stol);
else return 0;
5477 x = mDIFC*x+mDIFA*x*x;
5478 VFN_DEBUG_MESSAGE(
"PowderPattern::STOL2X("<<stol<<
","<<1.0/(2*stol)<<
")="<<x,2)
5483 if(abs(x)<1.0) x=2*asin(x);
else x=2*M_PI;
5491 if(this->
GetRadiation().GetWavelengthType()==WAVELENGTH_TOF)
5493 if(abs(mDIFA)>abs(mDIFC*1e-6))
5495 const REAL delta=mDIFC*mDIFC+4.0*mDIFA*x;
5496 stol = (-mDIFC+sqrt(delta))/(2.0*mDIFA);
5497 stol = 1/(2.0*stol);
5499 else stol=mDIFC/(2.0*x);
5503 VFN_DEBUG_MESSAGE(
"PowderPattern::X2STOL("<<x<<
")="<<stol,2)
5517 if(this->
GetRadiation().GetWavelengthType()!=WAVELENGTH_TOF)
5520 for(finish=0;finish<nb;++finish)
5529 for(start=nb-1;start>=0;--start)
5536 unsigned int width_golay=10;
5538 CrystVector_REAL obs;
5539 const int nbwidth=9;
5540 CrystVector_long width(nbwidth);
5550 if(max>=obs.numElements()) max = obs.numElements();
5551 for(
long j=min;j<max;j++) obs(j) = 0;
5553 const long nb=obs.numElements();
5554 for(
int j=0;j<nbwidth;j++)
5556 const long imax=obs.imax(nb/10,nb-1);
5557 const REAL iobs_max=obs(imax);
5558 REAL thres=iobs_max;
5560 for(i=imax-100;i<(imax+100);++i)
5562 if(i<0){i=0;
continue;}
5564 if(obs(i)<thres) thres=obs(i);
5566 thres=(iobs_max+thres)/2;
5568 while(obs(i)>=thres)
5576 while(obs(i)>=thres)
5583 cout<<endl<<
" => "<<width(j)<<endl;
5584 for(i=imax-width(j)*2;i<(imax+width(j)*2);++i)
5591 cout<<
"Width of "<<nbwidth<<
" strongest peaks:"<<endl<<width;
5592 width_golay=width(SortSubs(width)(nbwidth/2));
5593 cout<<
"median width:"<<width_golay<<endl;
5594 if(width_golay<=4)width_golay=4;
5595 if(width_golay>=16)width_golay=16;
5599 CrystVector_REAL obsd2;
5608 if(max>=obsd2.numElements()) max = obsd2.numElements();
5609 for(
long j=min;j<max;j++) obsd2(j) = 0;
5611 const float norm=-obsd2.min();
5618 CrystVector_REAL tmp;
5620 tmp.resizeAndPreserve(tmp.numElements()/4);
5622 QuickSortSubs(tmp,sub,tmp.numElements()-1,0);
5623 min_iobs=5*(tmp(tmp.numElements()/2)-tmp(tmp.numElements()/4));
5630 const long imax=obsd2.imax(start,finish);
5631 REAL iobs=obsd2(imax);
5633 REAL xmax=
mX(imax)*iobs;
5636 REAL lastiobs=obsd2(i);
5637 const REAL iobs_max=lastiobs;
5642 if(obsd2(--i)>=lastiobs)
break;
5646 xmax+=
mX(i)*lastiobs;nbav++;
5647 if(lastiobs<=0)
break;
5650 float dleft=
mX(i+1);
5655 if(i>=(nb-2))
break;
5656 if(obsd2(++i)>=lastiobs)
break;
5660 xmax+=
mX(i)*lastiobs;nbav++;
5661 if(lastiobs<=0)
break;
5664 float dright=
mX(i-1);
5667 REAL dmax=this->
X2STOL(xmax)*2;
5668 dright=this->
X2STOL(dright)*2;
5669 dleft =this->
X2STOL(dleft)*2;
5674 min_iobs=iobs/nbav*maxratio;
5680 if(nbav_min<3)nbav_min=3;
5686 cout<<
"Peak #"<<pl.
GetPeakList().size()<<
"imax="<<imax<<
", x="<<xmax<<
",d="<<1/dmax<<
", d2iobs_max="<<iobs_max
5687 <<
", d2Iobs="<<iobs<<
", nbav="<<nbav<<
", min_iobs="<<min_iobs<<
",sigma="<<sigma<<endl;
5688 if( ((nbav>=nbav_min) &&(iobs_max>min_iobs)&&((iobs/nbav)>min_iobs))
5689 ||((nbav>=nbav_min) &&(iobs_max>min_iobs)&&((iobs/nbav)>min_iobs*.2)&&((iobs/nbav)>3*sigma))
5690 ||((nbav>=nbav_min/2)&&(iobs_max>min_iobs)&&((iobs/nbav)>min_iobs*2 )&&((iobs/nbav)>6*sigma)))
5692 pl.
AddPeak(dmax,iobs,abs(dright-dleft)*.25);
5693 if((pl.
GetPeakList().size()==1)&&(maxratio<0)&&(min_iobs<0.005*iobs/nbav)) min_iobs=0.005*iobs/nbav;
5707 exportAtom(
string n,REAL X, REAL Y, REAL Z,REAL b,REAL o,
const ScatteringPower *pow):
5708 name(n),x(X),y(Y),z(Z),biso(b),occ(o),occMult(1),mpScattPow(pow){}
5710 REAL x,y,z,biso,occ;
5713 const ScatteringPower *mpScattPow;
5718 exportBond(
const string &a1,
const string &a2, REAL d, REAL s):
5719 at1(a1),at2(a2),dist(d),sigma(s){}
5726 exportAngle(
const string &a1,
const string &a2,
const string &a3, REAL a, REAL s):
5727 at1(a1),at2(a2),at3(a3),ang(a),sigma(s){}
5736 vector<const PowderPatternDiffraction*> vDiff;
5744 if((pBackground==0)||vDiff.size()==0)
return;
5747 ofstream dat((prefix+
".dat").c_str());
5749 <<
"INTER 1.0 1.0 0"<<endl<<endl<<endl<<endl<<endl;
5751 CrystVector_REAL ttheta;
5753 if(this->
GetRadiation().GetWavelengthType()!=WAVELENGTH_TOF) ttheta *= RAD2DEG;
5758 ofstream pcr((prefix+
".pcr").c_str());
5762 pcr<<
"Fox/ObjCryst exported file:"<<this->
GetName()<<endl;
5764 pcr<<
"NPATT 1"<<endl;
5766 pcr<<
"W_PAT 1.0"<<endl;
5768 pcr<<
"! Nph Dum Ias Nre Cry Opt Aut"<<endl;
5769 pcr<<
" "<<vDiff.size()
5770 <<
" 0 0 0 0 1 1 "<<endl;
5777 pcr<<
"! Job Npr Nba Nex Nsc Nor Iwg Ilo Res Ste Uni Cor Anm"<<endl
5779 <<
" 5 "<<pBackground->GetInterpPoints().first->numElements()
5780 <<
" 0 0 1 0 0 0 1 0 0 0"<<endl;
5785 std::string::size_type idx =prefix.rfind(
"/");
5786 std::string::size_type idx2=prefix.rfind(
"\\");
5787 std::string::size_type idx3=prefix.rfind(
":");
5788 if(((
long)idx2!=(
long)string::npos)&&((
long)idx2>(
long)idx))idx=idx2;
5789 if(((
long)idx3!=(
long)string::npos)&&((
long)idx3>(
long)idx))idx=idx3;
5790 if(idx==string::npos)
5793 shortName=prefix.substr(idx+1);
5795 pcr<<
"! File names of data files"<<endl;
5796 pcr<<shortName<<
".dat"<<endl;
5798 pcr<<
"! Mat Pcr Syo Rpa Sym Sho"<<endl
5799 <<
" 1 1 0 -1 0 0 "<<endl;
5801 pcr<<
"! Ipr Ppl Ioc Ls1 Ls2 Ls3 Prf Ins Hkl Fou Ana"<<endl
5802 <<
" 0 0 0 0 0 0 1 10 0 0 1 "<<endl;
5806 pcr<<
"!lambda1 lambda2 Ratio Bkpos Wdt Cthm muR AsyLim Rpolarz -> Patt #1"<<endl
5809 <<
" 0 0 0 0.95"<<endl;
5811 pcr<<
"!NCY Eps R_at R_an R_pr R_gl"<<endl
5812 <<
" 5 0.2 1.0 1.0 1.0 1.0"<<endl;
5814 pcr<<
"! Thmin Step Thmax PSD Sent0 -> Patt #1"<<endl
5815 <<
" 0 0 0 0 0"<<endl;
5817 pcr<<
"!2Theta Background for Pattern #1"<<endl;
5818 for(
long i=0;i<pBackground->GetInterpPoints().first->numElements();i++)
5819 pcr<<(*(pBackground->GetInterpPoints().first))(i)*RAD2DEG<<
" "
5820 <<(*(pBackground->GetInterpPoints().second))(i)<<
" 0.0"<<endl;
5822 pcr<<
"!"<<endl<<
"!"<<endl<<
"1 !Number of refined parameters"<<endl;
5824 pcr<<
"! Zero Code Sycos Code Sysin Code Lambda Code More -> Patt #1"<<endl;
5825 pcr<<
" "<<
mXZero*RAD2DEG <<
" 0.0 "
5826 <<m2ThetaDisplacement*RAD2DEG <<
" 0.0 "
5827 <<m2ThetaTransparency*RAD2DEG <<
" 0.0 "
5828 <<
"0.000 0.0 0"<<endl;
5830 for(
unsigned int i=0;i<vDiff.size();++i)
5832 pcr<<
"!-------------------------------------------------------------------------------"<<endl
5833 <<
"! Data for PHASE number: "<<i<<
" ==> Current R_Bragg for Pattern# 1: 0.00 "<<endl
5834 <<
"!-------------------------------------------------------------------------------"<<endl;
5836 pcr<<vDiff[i]->GetCrystal().GetName()<<endl;
5839 map<int,exportAtom> vExportAtom;
5840 list<exportBond> vExportBond;
5841 list<exportAngle> vExportAngle;
5843 CrystMatrix_REAL minDistTable;
5844 minDistTable=vDiff[i]->GetCrystal().GetMinDistanceTable(-1.);
5849 for(
int s=0;s<vDiff[i]->GetCrystal().GetScattererRegistry().GetNb();s++)
5856 if(vDiff[i]->GetCrystal().GetScatt(s).GetClassName()==
"Molecule")
5857 pMol=
dynamic_cast<const Molecule*
>(&(vDiff[i]->GetCrystal().GetScatt(s)));
5858 map<const MolAtom*,string> vMolAtomName;
5862 if(0==list(j).mpScattPow)
continue;
5863 bool redundant=
false;
5864 for(
unsigned long l=0;l<k;++l)
5865 if(abs(minDistTable(l,k))<0.5)
5867 map<int,exportAtom>::iterator pos=vExportAtom.find(l);
5868 if(pos!=vExportAtom.end()) pos->second.occMult+=1;
5875 name<<list(j).mpScattPow->GetName()<<k+1;
5876 vExportAtom.insert(make_pair(k,exportAtom(name.str(),
5877 list(j).
mX,list(j).mY,list(j).mZ,
5878 list(j).mpScattPow->GetBiso(),
5879 list(j).mOccupancy*list0(k).mDynPopCorr,
5880 list(j).mpScattPow)));
5881 if(pMol!=0) vMolAtomName.insert(make_pair(pMol->GetAtomList()[j],name.str()));
5887 for(vector<MolBond*>::const_iterator pos=pMol->GetBondList().begin();
5888 pos!=pMol->GetBondList().end();++pos)
5890 map<const MolAtom*,string>::const_iterator p1,p2;
5891 p1=vMolAtomName.find(&((*pos)->GetAtom1()));
5892 p2=vMolAtomName.find(&((*pos)->GetAtom2()));
5893 if( (p1!=vMolAtomName.end()) && (p2!=vMolAtomName.end()))
5894 vExportBond.push_back(exportBond(p1->second, p2->second,
5895 (*pos)->GetLength0(),(*pos)->GetLengthSigma()));
5898 for(vector<MolBondAngle*>::const_iterator pos=pMol->GetBondAngleList().begin();
5899 pos!=pMol->GetBondAngleList().end();++pos)
5901 map<const MolAtom*,string>::const_iterator p1,p2,p3;
5902 p1=vMolAtomName.find(&((*pos)->GetAtom1()));
5903 p2=vMolAtomName.find(&((*pos)->GetAtom2()));
5904 p3=vMolAtomName.find(&((*pos)->GetAtom3()));
5905 if( (p1!=vMolAtomName.end()) && (p2!=vMolAtomName.end()) && (p3!=vMolAtomName.end()))
5906 vExportAngle.push_back(exportAngle(p1->second, p2->second,p3->second,
5907 (*pos)->GetAngle0(),(*pos)->GetAngleSigma()));
5918 pcr<<
"!Nat Dis Ang Jbt Isy Str Furth ATZ Nvk More"<<endl
5919 << vExportAtom.size() <<
" "<<vExportBond.size()<<
" "<<vExportAngle.size()
5920 <<
" 0 0 0 0 1.0 0 1"<<endl;
5921 pcr<<
"!Jvi Jdi Hel Sol Mom Ter N_Domains"<<endl
5922 <<
" 0 3 0 0 0 0 0"<<endl;
5924 pcr<<
"!Contributions (0/1) of this phase to the patterns"<<endl
5927 pcr<<
"!Irf Npr Jtyp Nsp_Ref Ph_Shift for Pattern#"<<i<<endl
5928 <<
" 0 0 0 0 0"<<endl;
5929 pcr<<
"! Pr1 Pr2 Pr3 Brind. Rmua Rmub Rmuc for Pattern#"<<i<<endl
5930 <<
" 1.0 1.0 1.0 1.0 0.0 0.0 0.0"<<endl;
5933 pcr<<
"! Max_dst(dist) (angles) Bond-Valence Calc."<<endl
5934 <<
" 2.7000 1.5000 0"<<endl;
5937 pcr<<vDiff[i]->GetCrystal().GetSpaceGroup().GetCCTbxSpg().match_tabulated_settings().hermann_mauguin()
5938 <<
" <- Space Group Symbol"<<endl;
5940 pcr<<
"!Atom Typ X Y Z Biso Occ In Fin N_t Spc / Codes"<<endl;
5941 for(map<int,exportAtom>::const_iterator pos=vExportAtom.begin();pos!=vExportAtom.end();++pos)
5943 pcr<<pos->second.name
5944 <<
" "<<pos->second.mpScattPow->GetSymbol()<<
" "
5945 <<pos->second.x<<
" "<<pos->second.y<<
" "<<pos->second.z<<
" "
5946 <<pos->second.biso<<
" "
5947 <<pos->second.occ*pos->second.occMult<<
" 0 0 0 0"<<endl
5948 <<
" 0 0 0 0 0"<<endl;
5951 REAL eta0=vDiff[0]->GetProfile().GetPar(
"Eta0").GetHumanValue();
5952 if(eta0<.01) eta0=.01;
5953 else if(eta0>.99) eta0=.99;
5954 pcr<<
"!Scale Shape1 Bov Str1 Str2 Str3 Strain-Model"<<endl
5956 <<
" 0.0 0.0 0.0 0.0 0"<<endl
5957 <<
" 1.0 0.0 0.0 0.0 0.0 0.0 0"<<endl;
5960 pcr<<
"! U V W X Y GauSiz LorSiz Size-Model"<<endl;
5963 pcr<<vDiff[i]->GetProfile().GetPar(
"U").GetHumanValue()<<
" ";
5964 pcr<<vDiff[i]->GetProfile().GetPar(
"V").GetHumanValue()<<
" ";
5965 pcr<<vDiff[i]->GetProfile().GetPar(
"W").GetHumanValue()<<
" ";
5966 pcr<<
" 0.0 0.0 0.0 0.0 "<<endl
5967 <<
" 0.0 0.0 0.0 0.0 0.0 0.0 0.0 "<<endl;
5969 pcr<<
"! a b c alpha beta gamma #Cell Info"<<endl
5970 <<vDiff[i]->GetCrystal().GetLatticePar(0)<<
" "
5971 <<vDiff[i]->GetCrystal().GetLatticePar(1)<<
" "
5972 <<vDiff[i]->GetCrystal().GetLatticePar(2)<<
" "
5973 <<vDiff[i]->GetCrystal().GetLatticePar(3)*RAD2DEG<<
" "
5974 <<vDiff[i]->GetCrystal().GetLatticePar(4)*RAD2DEG<<
" "
5975 <<vDiff[i]->GetCrystal().GetLatticePar(5)*RAD2DEG<<endl
5976 <<
" 0.0 0.0 0.0 0.0 0.0 0.0"<<endl;
5977 pcr<<
"! Pref1 Pref2 alpha0 beta0 alpha1 beta1 ?"<<endl
5978 <<
" 0.0 0.0 0.0 0.0 0.0 0.0"<<endl
5979 <<
" 0.0 0.0 0.0 0.0 0.0 0.0"<<endl;
5983 if(vExportBond.size()>0)
5985 pcr<<
"!Soft distance constraints"<<endl;
5986 for(list<exportBond>::const_iterator pos=vExportBond.begin();pos!=vExportBond.end();++pos)
5988 pcr<<pos->at1<<
" "<<pos->at2<<
" 1 0 0 0 "<<pos->dist<<
" "<<pos->sigma<<endl;
5991 if(vExportBond.size()>0)
5993 pcr<<
"!Soft angle constraints"<<endl;
5994 for(list<exportAngle>::const_iterator pos=vExportAngle.begin();pos!=vExportAngle.end();++pos)
5996 pcr<<pos->at1<<
" "<<pos->at2<<
" "<<pos->at3<<
" 1 1 0 0 0 0 0 0 "
5997 <<pos->ang*RAD2DEG<<
" "<<pos->sigma*RAD2DEG<<endl;
6013 TAU_PROFILE(
"PowderPattern::CalcPowderPattern()",
"void ()",TAU_DEFAULT);
6014 VFN_DEBUG_ENTRY(
"PowderPattern::CalcPowderPattern()",3);
6015 if(mPowderPatternComponentRegistry.GetNb()==0)
6026 for(
unsigned long j=0;j<
mNbPoint;j++)
6028 if(*p0 <=0) {*p1 =0.;}
6029 else *p1 = 1. / *p0;
6032 VFN_DEBUG_EXIT(
"PowderPattern::CalcPowderPattern():no components!",3);
6035 TAU_PROFILE_TIMER(timer1,
"PowderPattern::CalcPowderPattern1()Calc components",
"", TAU_FIELD);
6036 TAU_PROFILE_TIMER(timer2,
"PowderPattern::CalcPowderPattern2(Add spectrums-scaled)"\
6038 TAU_PROFILE_TIMER(timer3,
"PowderPattern::CalcPowderPattern2(Add spectrums-backgd1)"\
6040 TAU_PROFILE_TIMER(timer4,
"PowderPattern::CalcPowderPattern2(Add spectrums-backgd2)"\
6042 TAU_PROFILE_TIMER(timer5,
"PowderPattern::CalcPowderPattern3(Variance)",
"", TAU_FIELD);
6043 TAU_PROFILE_START(timer1);
6044 for(
int i=0;i<mPowderPatternComponentRegistry.GetNb();i++)
6045 mPowderPatternComponentRegistry.GetObj(i).CalcPowderPattern();
6046 TAU_PROFILE_STOP(timer1);
6047 VFN_DEBUG_MESSAGE(
"PowderPattern::CalcPowderPattern():Calculated components..",3);
6054 for(
int i=0;i<mPowderPatternComponentRegistry.GetNb();i++)
6056 mPowderPatternComponentRegistry.GetObj(i).GetClockPowderPatternCalc())
6064 VFN_DEBUG_EXIT(
"PowderPattern::CalcPowderPattern():no need to recalc",3);
6069 for(
int i=0;i<mPowderPatternComponentRegistry.GetNb();i++)
6071 VFN_DEBUG_MESSAGE(
"PowderPattern::CalcPowderPattern():Adding "<< mPowderPatternComponentRegistry.GetObj(i).GetName(),3);
6072 if(
true==mPowderPatternComponentRegistry.GetObj(i).IsScalable())
6074 TAU_PROFILE_START(timer2);
6077 const REAL * p1=mPowderPatternComponentRegistry.GetObj(i)
6078 .mPowderPatternCalc.data();
6080 const REAL s = mScaleFactor(i);
6081 for(
unsigned long j=0;j<mNbPointUsed;j++) *p0++ = s * *p1++;
6086 const REAL * p1=mPowderPatternComponentRegistry.GetObj(i)
6087 .mPowderPatternCalc.data();
6089 const REAL s = mScaleFactor(i);
6090 for(
unsigned long j=0;j<mNbPointUsed;j++) *p0++ += s * *p1++;
6092 TAU_PROFILE_STOP (timer2);
6096 TAU_PROFILE_START(timer3);
6099 const REAL * p1=mPowderPatternComponentRegistry.GetObj(i)
6100 .mPowderPatternCalc.data();
6102 for(
unsigned long j=0;j<mNbPointUsed;j++) *p0++ = *p1++;
6107 const REAL * p1=mPowderPatternComponentRegistry.GetObj(i)
6108 .mPowderPatternCalc.data();
6110 for(
unsigned long j=0;j<mNbPointUsed;j++) *p0++ += *p1++;
6113 TAU_PROFILE_STOP(timer3);
6114 TAU_PROFILE_START(timer4);
6120 const REAL *p1=mPowderPatternComponentRegistry.GetObj(i)
6121 .mPowderPatternCalc.data();
6122 for(
unsigned long j=0;j<mNbPointUsed;j++) *p0++ = *p1++;
6128 const REAL *p1=mPowderPatternComponentRegistry.GetObj(i)
6129 .mPowderPatternCalc.data();
6130 for(
unsigned long j=0;j<mNbPointUsed;j++) *p0++ += *p1++;
6133 TAU_PROFILE_STOP(timer4);
6138 TAU_PROFILE_START(timer5);
6141 VFN_DEBUG_MESSAGE(
"PowderPattern::CalcPowderPattern():variance",2);
6144 for(
int i=0;i<mPowderPatternComponentRegistry.GetNb();i++)
6146 if(mPowderPatternComponentRegistry.GetObj(i).HasPowderPatternCalcVariance())
6148 if(0==mPowderPatternComponentRegistry.GetObj(i).GetPowderPatternCalcVariance()
6149 .numElements())
break;
6152 const REAL *p1=mPowderPatternComponentRegistry.GetObj(i)
6153 .GetPowderPatternCalcVariance().data();
6155 if(
true==mPowderPatternComponentRegistry.GetObj(i).IsScalable())
6157 const REAL s2 = mScaleFactor(i) * mScaleFactor(i);
6158 for(
unsigned long j=0;j<mNbPointUsed;j++) *p0++ += s2 * *p1++;
6160 else for(
unsigned long j=0;j<mNbPointUsed;j++) *p0++ += *p1++;
6165 for(
unsigned long j=0;j<mNbPointUsed;j++)
6166 if(*p1 <=0) {*p0++ =0;p1++;}
6167 else *p0++ = 1. / *p1++;
6170 TAU_PROFILE_STOP(timer5);
6171 VFN_DEBUG_EXIT(
"PowderPattern::CalcPowderPattern():End",3);
6174 void PowderPattern::CalcPowderPattern_FullDeriv(std::set<RefinablePar*> &vPar)
6176 TAU_PROFILE(
"PowderPattern::CalcPowderPattern_FullDeriv()",
"void ()",TAU_DEFAULT);
6178 mPowderPattern_FullDeriv.clear();
6179 if(mPowderPatternComponentRegistry.GetNb()==0)
return;
6180 std::vector<std::map<RefinablePar*,CrystVector_REAL>*> comps;
6181 for(
int i=0;i<mPowderPatternComponentRegistry.GetNb();i++)
6182 comps.push_back(&(mPowderPatternComponentRegistry.GetObj(i).GetPowderPattern_FullDeriv(vPar)));
6184 for(std::set<RefinablePar *>::iterator par=vPar.begin();par!=vPar.end();++par)
6186 if(*par==0)
continue;
6187 for(
int i=0;i<mPowderPatternComponentRegistry.GetNb();i++)
6189 if((*par)->GetPointer()==mScaleFactor.data()+i)
6191 mPowderPattern_FullDeriv[*par]=mPowderPatternComponentRegistry.GetObj(i).GetPowderPatternCalc();
6196 if((*(comps[i]))[*par].size()==0)
continue;
6198 if(
true==mPowderPatternComponentRegistry.GetObj(i).IsScalable())
6200 if(mPowderPattern_FullDeriv[*par].size()==0)
6202 mPowderPattern_FullDeriv[*par].resize(
mNbPoint);
6203 const REAL * p1=(*comps[i])[*par].data();
6204 REAL * p0 = mPowderPattern_FullDeriv[*par].data();
6205 const REAL s = mScaleFactor(i);
6206 for(
unsigned long j=0;j<mNbPointUsed;j++) *p0++ = s * *p1++;
6207 for(
unsigned long j=mNbPointUsed;j<
mNbPoint;j++) *p0++ = 0;
6211 const REAL * p1=(*comps[i])[*par].data();
6212 REAL * p0 = mPowderPattern_FullDeriv[*par].data();
6213 const REAL s = mScaleFactor(i);
6214 for(
unsigned long j=0;j<mNbPointUsed;j++) *p0++ += s * *p1++;
6219 if(mPowderPattern_FullDeriv[*par].size()==0)
6221 mPowderPattern_FullDeriv[*par].resize(
mNbPoint);
6222 const REAL * p1=(*comps[i])[*par].data();
6223 REAL * p0 = mPowderPattern_FullDeriv[*par].data();
6224 for(
unsigned long j=0;j<mNbPointUsed;j++) *p0++ = *p1++;
6225 for(
unsigned long j=mNbPointUsed;j<
mNbPoint;j++) *p0++ = 0;
6229 const REAL * p1=(*comps[i])[*par].data();
6230 REAL * p0 = mPowderPattern_FullDeriv[*par].data();
6231 for(
unsigned long j=0;j<mNbPointUsed;j++) *p0++ += *p1++;
6238 std::map<RefinablePar*, CrystVector_REAL> oldDeriv;
6239 std::vector<const CrystVector_REAL*> v;
6241 for(std::map<RefinablePar*, CrystVector_REAL>::reverse_iterator pos=mPowderPattern_FullDeriv.rbegin();pos!=mPowderPattern_FullDeriv.rend();++pos)
6243 if(pos->first==0)
continue;
6244 if(pos->second.size()==0)
continue;
6246 const REAL step=pos->first->GetDerivStep();
6247 pos->first->Mutate(step);
6250 pos->first->Mutate(-2*step);
6253 oldDeriv[pos->first]/=2*step;
6254 pos->first->Mutate(step);
6256 v.push_back(&(pos->second));
6257 v.push_back(&(oldDeriv[pos->first]));
6258 cout<<pos->first->GetName()<<
":"<<pos->second.size()<<
","<<oldDeriv[pos->first].size()<<endl;
6261 cout<<FormatVertVector<REAL>(v,16)<<endl;
6272 TAU_PROFILE(
"PowderPattern::CalcPowderPatternIntegrated()",
"void ()",TAU_DEFAULT);
6273 VFN_DEBUG_ENTRY(
"PowderPattern::CalcPowderPatternIntegrated()",4);
6274 if(mPowderPatternComponentRegistry.GetNb()==0)
6278 VFN_DEBUG_EXIT(
"PowderPattern::CalcPowderPatternIntegrated():no components!",4);
6281 TAU_PROFILE_TIMER(timer1,
"PowderPattern::CalcPowderPatternIntegrated()1:Calc components",\
6283 TAU_PROFILE_TIMER(timer2,
"PowderPattern::CalcPowderPatternIntegrated()2:Add comps-scaled"\
6285 TAU_PROFILE_TIMER(timer3,
"PowderPattern::CalcPowderPatternIntegrated()2:Add backgd1"\
6287 TAU_PROFILE_TIMER(timer4,
"PowderPattern::CalcPowderPatternIntegrated()2:Add backgd2"\
6289 TAU_PROFILE_TIMER(timer5,
"PowderPattern::CalcPowderPatternIntegrated()3:Variance"
6291 TAU_PROFILE_START(timer1);
6292 vector< pair<const CrystVector_REAL*,const RefinableObjClock*> > comps;
6293 for(
int i=0;i<mPowderPatternComponentRegistry.GetNb();i++)
6295 comps.push_back(mPowderPatternComponentRegistry.GetObj(i).
6296 GetPowderPatternIntegratedCalc());
6298 TAU_PROFILE_STOP(timer1);
6305 for(vector< pair<const CrystVector_REAL*,const RefinableObjClock*> >::iterator
6306 pos=comps.begin();pos!=comps.end();++pos)
6315 VFN_DEBUG_EXIT(
"PowderPattern::CalcPowderPatternIntegrated():no need to recalc",4);
6318 VFN_DEBUG_MESSAGE(
"PowderPattern::CalcPowderPatternIntegrated():Recalc",3);
6321 for(
int i=0;i<mPowderPatternComponentRegistry.GetNb();i++)
6323 VFN_DEBUG_MESSAGE(
"PowderPattern::CalcPowderPatternIntegrated():Adding "\
6324 << mPowderPatternComponentRegistry.GetObj(i).GetName(),3);
6325 if(
true==mPowderPatternComponentRegistry.GetObj(i).IsScalable())
6327 TAU_PROFILE_START(timer2);
6330 const REAL * RESTRICT p1= comps[i].first->data();
6332 const REAL s = mScaleFactor(i);
6333 for(
unsigned long j=mNbIntegrationUsed;j>0;j--) *p0++ = s * *p1++;
6337 const REAL * RESTRICT p1= comps[i].first->data();
6339 const REAL s = mScaleFactor(i);
6340 for(
unsigned long j=mNbIntegrationUsed;j>0;j--) *p0++ += s * *p1++;
6342 TAU_PROFILE_STOP (timer2);
6346 TAU_PROFILE_START(timer3);
6349 const REAL * RESTRICT p1= comps[i].first->data();
6351 for(
unsigned long j=mNbIntegrationUsed;j>0;j--) *p0++ = *p1++;
6355 const REAL * RESTRICT p1= comps[i].first->data();
6357 for(
unsigned long j=mNbIntegrationUsed;j>0;j--) *p0++ += *p1++;
6360 TAU_PROFILE_STOP(timer3);
6361 TAU_PROFILE_START(timer4);
6366 const REAL * RESTRICT p1= comps[i].first->data();
6368 for(
unsigned long j=mNbIntegrationUsed;j>0;j--) *p0++ = *p1++;
6372 const REAL * RESTRICT p1= comps[i].first->data();
6374 for(
unsigned long j=mNbIntegrationUsed;j>0;j--) *p0++ += *p1++;
6377 TAU_PROFILE_STOP(timer4);
6380 TAU_PROFILE_START(timer5);
6383 bool useCalcVariance=
false;
6384 for(
int i=0;i<mPowderPatternComponentRegistry.GetNb();i++)
6385 if(mPowderPatternComponentRegistry.GetObj(i).HasPowderPatternCalcVariance())
6386 if(mPowderPatternComponentRegistry.GetObj(i)
6387 .GetPowderPatternIntegratedCalcVariance().first->numElements() !=0)
6388 useCalcVariance=
true;
6391 VFN_DEBUG_MESSAGE(
"PowderPattern::CalcPowderPatternIntegrated():variance",3);
6394 mIntegratedWeight.resize(mNbIntegrationUsed);
6395 const REAL * RESTRICT p1= mIntegratedVarianceObs.data();
6397 for(
unsigned long j=mNbIntegrationUsed;j>0;j--) *p0++ = *p1++;
6401 for(
int i=0;i<mPowderPatternComponentRegistry.GetNb();i++)
6403 if(mPowderPatternComponentRegistry.GetObj(i).HasPowderPatternCalcVariance())
6405 if(0==mPowderPatternComponentRegistry.GetObj(i)
6406 .GetPowderPatternIntegratedCalcVariance().first->numElements())
break;
6408 const REAL * RESTRICT p1= mPowderPatternComponentRegistry.GetObj(i)
6409 .GetPowderPatternIntegratedCalcVariance().first->data();
6414 if(
true==mPowderPatternComponentRegistry.GetObj(i).IsScalable())
6416 const REAL s2 = mScaleFactor(i) * mScaleFactor(i);
6417 for(
unsigned long j=mNbIntegrationUsed;j>0;j--) *p0++ += s2 * *p1++;
6419 else for(
unsigned long j=mNbIntegrationUsed;j>0;j--) *p0++ += *p1++;
6424 REAL *p0 = mIntegratedWeight.data();
6426 for(
unsigned long j=mNbIntegrationUsed;j>0;j--)
6427 if(*p1 <=0) {*p0++ =0;p1++;}
6428 else *p0++ = 1. / *p1++;
6430 else mIntegratedWeight.resize(0);
6432 TAU_PROFILE_STOP(timer5);
6458 VFN_DEBUG_EXIT(
"PowderPattern::CalcPowderPatternIntegrated():End",4);
6461 void PowderPattern::CalcPowderPatternIntegrated_FullDeriv(std::set<RefinablePar *> &vPar)
6463 TAU_PROFILE(
"PowderPattern::CalcPowderPatternIntegrated_FullDeriv()",
"void ()",TAU_DEFAULT);
6468 mPowderPatternUsed_FullDeriv.clear();
6469 if(mPowderPatternComponentRegistry.GetNb()==0)
return;
6470 std::vector<map<RefinablePar*,CrystVector_REAL>*> comps;
6471 for(
int i=0;i<mPowderPatternComponentRegistry.GetNb();i++)
6473 comps.push_back(&(mPowderPatternComponentRegistry.GetObj(i).
6474 GetPowderPatternIntegrated_FullDeriv(vPar)));
6478 mPowderPatternIntegrated_FullDeriv.clear();
6479 for(std::set<RefinablePar *>::iterator par=vPar.begin();par!=vPar.end();++par)
6481 if(*par==0)
continue;
6482 for(
int i=0;i<mPowderPatternComponentRegistry.GetNb();i++)
6484 if((*par)->GetPointer()==mScaleFactor.data()+i)
6487 mPowderPatternIntegrated_FullDeriv[*par]=*(mPowderPatternComponentRegistry.GetObj(i).GetPowderPatternIntegratedCalc().first);
6493 if((*(comps[i]))[*par].size()==0)
continue;
6495 if(
true==mPowderPatternComponentRegistry.GetObj(i).IsScalable())
6497 if(mPowderPatternIntegrated_FullDeriv[*par].size()==0)
6499 mPowderPatternIntegrated_FullDeriv[*par].resize(mNbIntegrationUsed);
6500 const REAL * RESTRICT p1= (*(comps[i]))[*par].data();
6501 REAL * RESTRICT p0 = mPowderPatternIntegrated_FullDeriv[*par].data();
6502 const REAL s = mScaleFactor(i);
6503 for(
unsigned long j=mNbIntegrationUsed;j>0;j--)
6509 if((j==mNbIntegrationUsed)&&((*par)->GetName()==
"Cimetidine_C11_x")) cout<<__FILE__<<
":"<<__LINE__<<
":"<<*p0<<
","<<s<<
"*"<<*p1<<endl;
6516 const REAL * RESTRICT p1= (*(comps[i]))[*par].data();
6517 REAL * RESTRICT p0 = mPowderPatternIntegrated_FullDeriv[*par].data();
6518 const REAL s = mScaleFactor(i);
6519 for(
unsigned long j=mNbIntegrationUsed;j>0;j--)
6525 if((j==mNbIntegrationUsed)&&((*par)->GetName()==
"Cimetidine_C11_x")) cout<<__FILE__<<
":"<<__LINE__<<
":"<<*p0<<
","<<s<<
"*"<<*p1<<endl;
6533 if(mPowderPatternIntegrated_FullDeriv[*par].size()==0)
6535 mPowderPatternIntegrated_FullDeriv[*par].resize(mNbIntegrationUsed);
6536 const REAL * RESTRICT p1= (*(comps[i]))[*par].data();
6537 REAL * RESTRICT p0 = mPowderPatternIntegrated_FullDeriv[*par].data();
6538 for(
unsigned long j=mNbIntegrationUsed;j>0;j--)
6544 if((j==mNbIntegrationUsed)&&((*par)->GetName()==
"Cimetidine_C11_x")) cout<<__FILE__<<
":"<<__LINE__<<
":"<<*p0<<
","<<*p1<<endl;
6551 const REAL * RESTRICT p1= (*(comps[i]))[*par].data();
6552 REAL * RESTRICT p0 = mPowderPatternIntegrated_FullDeriv[*par].data();
6553 for(
unsigned long j=mNbIntegrationUsed;j>0;j--)
6559 if((j==mNbIntegrationUsed)&&((*par)->GetName()==
"Cimetidine_C11_x")) cout<<__FILE__<<
":"<<__LINE__<<
":"<<*p0<<
","<<*p1<<endl;
6566 if((*par)->GetName()==
"Cimetidine_C11_x")
6567 cout<<
"PowderPattern::CalcPowderPatternIntegrated_FullDeriv():"
6568 <<(*par)->GetName()<<
":s="<<mScaleFactor(i)<<
", d[0]="<<(*(comps[i]))[*par](0)
6569 <<
", integ="<<mPowderPatternIntegrated_FullDeriv[*par](0)<<endl;
6583 REFPAR_DERIV_STEP_ABSOLUTE,
true,
true,
true,
false,RAD2DEG);
6590 REFPAR_DERIV_STEP_ABSOLUTE,
true,
true,
true,
false,RAD2DEG);
6597 REFPAR_DERIV_STEP_ABSOLUTE,
true,
true,
true,
false,RAD2DEG);
6604 REFPAR_DERIV_STEP_ABSOLUTE,
true,
true,
true,
false,1.0);
6611 REFPAR_DERIV_STEP_ABSOLUTE,
true,
true,
true,
false,1.0);
6618 REFPAR_DERIV_STEP_ABSOLUTE,
true,
true,
true,
false,1.0);
6626 bool needPrep=
false;
6627 for(
int i=0;i<mPowderPatternComponentRegistry.GetNb();i++)
6629 mPowderPatternComponentRegistry.GetObj(i).GetBraggLimits();
6630 if(mPowderPatternComponentRegistry.GetObj(i).GetClockBraggLimits()
6631 >mClockIntegratedFactorsPrep)
6640 if(mClockIntegratedFactorsPrep<mClockNbPointUsed) needPrep=
true;
6642 if(
false==needPrep)
return;
6643 VFN_DEBUG_ENTRY(
"PowderPattern::PrepareIntegratedRfactor()",3);
6644 TAU_PROFILE(
"PowderPattern::PrepareIntegratedRfactor()",
"void ()",TAU_DEFAULT);
6650 for(
int i=0;i<mPowderPatternComponentRegistry.GetNb();i++)
6652 const CrystVector_long vLim=mPowderPatternComponentRegistry.GetObj(i).GetBraggLimits();
6653 for(
int j=0;j<vLim.numElements();j++) vLimits.push_back(vLim(j));
6655 if(vLimits.size()<2)
6657 mIntegratedPatternMin.resize(0);
6658 mIntegratedPatternMax.resize(0);
6659 mNbIntegrationUsed=0;
6660 mClockIntegratedFactorsPrep.Click();
6662 VFN_DEBUG_EXIT(
"PowderPattern::PrepareIntegratedRfactor(): no intervals !",3);
6665 if(*(vLimits.begin())<0)
6667 vLimits.push_back(0);
6670 for(list<long>::iterator pos=vLimits.begin();pos!=vLimits.end();)
6672 if( (*pos<0) || (*pos>=
long(mNbPointUsed)) ) pos=vLimits.erase(pos);
6677 list<long> vLimits2;
6678 list<long>::iterator pos1=vLimits.begin();
6679 list<long>::iterator pos2=pos1;pos2++;
6680 for(;pos2!=vLimits.end();)
6682 const long pix1=*pos1;
6684 vLimits2.push_back(pix1);
6688 if(pos2==vLimits.end())
break;
6689 if(*pos2>(pix1+8))
break;
6692 vLimits2.push_back(*pos1);
6695 pos1=vLimits2.begin();
6697 for(;pos2!=vLimits2.end();)
6700 if( *pos2<=((*pos1)+2))
6703 pos2=vLimits2.erase(pos2);
6706 else {pos1++;pos2++;}
6710 list<pair<long,long> > vLimits3;
6711 pos1=vLimits2.begin();
6713 for(;pos2!=vLimits2.end();)
6715 if(*pos2!=(
long(mNbPointUsed)-1)) vLimits3.push_back(make_pair(*pos1++,*pos2++-1));
6716 else vLimits3.push_back(make_pair(*pos1++,*pos2++));
6720 mIntegratedPatternMin.resize(vLimits3.size());
6721 mIntegratedPatternMax.resize(vLimits3.size());
6723 for(list<pair<long,long> >::iterator pos=vLimits3.begin();pos!=vLimits3.end();++pos)
6725 mIntegratedPatternMin(i)=pos->first;
6726 mIntegratedPatternMax(i++)=pos->second;
6729 long numInterval=mIntegratedPatternMin.numElements();
6730 CrystVector_bool keep(numInterval);
6737 VFN_DEBUG_MESSAGE(
"PowderPattern::PrepareIntegratedRfactor():5:Excluded regions("<<nbExclude<<
")",3);
6739 long minExcl,maxExcl;
6742 for(
int i=0;i<nbExclude;i++)
6744 while(mIntegratedPatternMax(j)<minExcl)
6747 if(j>=numInterval)
break;
6749 if(j>=numInterval)
break;
6750 while(mIntegratedPatternMin(j)<maxExcl)
6752 if( (mIntegratedPatternMin(j)>minExcl) &&(mIntegratedPatternMax(j)<maxExcl))
6754 if( (mIntegratedPatternMin(j)<minExcl) &&(mIntegratedPatternMax(j)<maxExcl))
6755 mIntegratedPatternMax(j)=minExcl;
6756 if( (mIntegratedPatternMin(j)>minExcl) &&(mIntegratedPatternMax(j)>maxExcl))
6757 mIntegratedPatternMin(j)=maxExcl;
6758 if(j==(numInterval-1))
break;
6765 while(mIntegratedPatternMax(j)>=minExcl)
6773 VFN_DEBUG_MESSAGE(
"PowderPattern::PrepareIntegratedRfactor():6",3);
6775 for(
int i=0;i<numInterval;i++)
6779 mIntegratedPatternMin(j )=mIntegratedPatternMin(i);
6780 mIntegratedPatternMax(j++)=mIntegratedPatternMax(i);
6784 mIntegratedPatternMax.resizeAndPreserve(numInterval);
6785 mIntegratedPatternMin.resizeAndPreserve(numInterval);
6787 VFN_DEBUG_MESSAGE(
"PowderPattern::PrepareIntegratedRfactor():intervals"<<endl\
6790 mIntegratedObs.resize(numInterval);
6791 mIntegratedVarianceObs.resize(numInterval);
6792 mIntegratedVarianceObs=0;
6794 mIntegratedWeightObs.resize(numInterval);
6795 for(
int i=0;i<numInterval;i++)
6797 for(
int j=mIntegratedPatternMin(i);j<=mIntegratedPatternMax(i);j++)
6802 if(mIntegratedVarianceObs(i) <= 0) mIntegratedWeightObs(i)=0;
6803 else mIntegratedWeightObs(i)=1./mIntegratedVarianceObs(i);
6810 mNbIntegrationUsed=mIntegratedPatternMin.numElements();
6811 mClockIntegratedFactorsPrep.Click();
6812 VFN_DEBUG_EXIT(
"PowderPattern::PrepareIntegratedRfactor()",3);
6819 if(this->
GetRadiation().GetWavelengthType()==WAVELENGTH_TOF)
6826 if(1>fabs(sinth)) tmp=(
unsigned long)(this->
X2PixelCorr(2*asin(sinth)));
else tmp=
mNbPoint;
6829 if(tmp !=mNbPointUsed)
6832 mClockNbPointUsed.Click();
6833 VFN_DEBUG_MESSAGE(
"PowderPattern::CalcNbPointUsed():"<<mNbPointUsed<<
" max(sin(theta)/lambda)="<<
mMaxSinThetaOvLambda, 3)
6840 VFN_DEBUG_MESSAGE(
"PowderPattern::InitOptions()",5)
6841 static string OptProfileIntegrationName;
6842 static string OptProfileIntegrationChoices[2];
6844 static bool needInitNames=
true;
6845 if(
true==needInitNames)
6847 OptProfileIntegrationName=
"Use Integrated Profiles";
6848 OptProfileIntegrationChoices[0]=
"Yes (recommended)";
6849 OptProfileIntegrationChoices[1]=
"No";
6851 needInitNames=
false;
6853 mOptProfileIntegration.Init(2,&OptProfileIntegrationName,OptProfileIntegrationChoices);
6854 this->
AddOption(&mOptProfileIntegration);
6857 #ifdef __WX__CRYST__
6862 return mpWXCrystObj;
6875 SPGScore::SPGScore(
const string &s,
const REAL r,
const REAL g,
const unsigned int nbextinct,
const REAL ngof,
const unsigned int nbrefl):
6876 hm(s),rw(r),gof(g),ngof(ngof),nbextinct446(nbextinct),nbreflused(nbrefl)
6882 if(s1.ngof > 0.00001 && s2.ngof >0.00001)
return s1.ngof < s2.ngof;
6883 return s1.gof < s2.gof;
6891 std::vector<bool> fingerprint(5*5*7-1+6);
6899 for(
long h=0;h<5;++h)
6900 for(
long k=0;k<5;++k)
6901 for (
long l=0;l<7;++l)
6903 if((h+k+l)==0)
continue;
6904 cctbx::miller::index<long> hkl(scitbx::vec3<long>(h,k,l));
6905 if(i>=fingerprint.size()) cout<<
"WHOOOOOOOOOOOOOPS"<<endl;
6906 fingerprint[i++] =spg.is_sys_absent(hkl);
6914 SpaceGroupExplorer::SpaceGroupExplorer(PowderPatternDiffraction *pd):
6917 SPGScore SpaceGroupExplorer::Run(
const string &spgId,
const bool fitprofile,
6918 const bool verbose,
const bool restore_orig,
6919 const bool update_display,
const REAL relative_length_tolerance,
6920 const REAL absolute_angle_tolerance_degree)
6922 cctbx::sgtbx::space_group sg;
6926 cctbx::sgtbx::space_group_symbols sgs=cctbx::sgtbx::space_group_symbols(spgId);
6927 sg = cctbx::sgtbx::space_group(sgs);
6929 catch(exception ex1)
6934 sg = cctbx::sgtbx::space_group(spgId);
6936 catch(exception ex2)
6938 string emsg =
"Space group symbol '" + spgId +
"' not recognized";
6939 throw ObjCrystException(emsg);
6942 return this->Run(sg, fitprofile, verbose, restore_orig, update_display,
6943 relative_length_tolerance, absolute_angle_tolerance_degree);
6946 SPGScore SpaceGroupExplorer::Run(
const cctbx::sgtbx::space_group &spg,
const bool fitprofile,
const bool verbose,
6947 const bool restore_orig,
const bool update_display,
6948 const REAL relative_length_tolerance,
const REAL absolute_angle_tolerance_degree)
6950 TAU_PROFILE(
"SpaceGroupExplorer::Run()",
"void (wxCommandEvent &)",TAU_DEFAULT);
6951 TAU_PROFILE_TIMER(timer1,
"SpaceGroupExplorer::Run()LSQ-P1",
"", TAU_FIELD);
6952 TAU_PROFILE_TIMER(timer2,
"SpaceGroupExplorer::Run()LSQ1",
"", TAU_FIELD);
6953 TAU_PROFILE_TIMER(timer3,
"SpaceGroupExplorer::Run()LSQ2",
"", TAU_FIELD);
6954 Crystal *pCrystal=&(mpDiff->GetCrystal());
6956 const REAL a=pCrystal->GetLatticePar(0),
6957 b=pCrystal->GetLatticePar(1),
6958 c=pCrystal->GetLatticePar(2),
6959 d=pCrystal->GetLatticePar(3),
6960 e=pCrystal->GetLatticePar(4),
6961 f=pCrystal->GetLatticePar(5);
6962 const string spghm=pCrystal->GetSpaceGroup().GetCCTbxSpg().match_tabulated_settings().hermann_mauguin();
6963 const string name=pCrystal->GetName();
6965 const cctbx::sgtbx::space_group_symbols s = spg.match_tabulated_settings();
6966 const string hm=s.universal_hermann_mauguin();
6967 const cctbx::uctbx::unit_cell uc(scitbx::af::double6(a,b,c,d*RAD2DEG,e*RAD2DEG,f*RAD2DEG));
6968 if(!spg.is_compatible_unit_cell(uc,relative_length_tolerance,absolute_angle_tolerance_degree))
6970 throw ObjCrystException(
"Spacegroup is not compatible with unit cell.");
6972 mpDiff->GetCrystal().Init(a,b,c,d,e,f,hm,name);
6973 mpDiff->SetExtractionMode(
true,
true);
6974 unsigned int nbcycle=1;
6975 if(update_display) mpDiff->GetParentPowderPattern().UpdateDisplay();
6977 unsigned int nbfreepar=mpDiff->GetProfileFitNetNbObs();
6978 if(nbfreepar<1) nbfreepar=1;
6982 lsq.SetRefinedObj(mpDiff->GetParentPowderPattern(),0,
true,
true);
6983 lsq.PrepareRefParList(
true);
6984 const unsigned int saved_par = lsq.GetCompiledRefinedObj().CreateParamSet(
"SpaceGroupExplorer saved parameters");
6985 lsq.GetCompiledRefinedObj().SaveParamSet(saved_par);
6988 for(
unsigned int j=0;j<nbcycle;j++)
6991 mpDiff->SetExtractionMode(
true,
true);
6992 const float t0=chrono.seconds();
6993 if(verbose) cout<<
"Doing Le Bail, t="<<
FormatFloat(t0,6,2)<<
"s";
6994 mpDiff->ExtractLeBail(5);
6995 if(verbose) cout<<
", dt="<<
FormatFloat(chrono.seconds()-t0,6,2)<<
"s"<<endl;
6999 TAU_PROFILE_START(timer2);
7000 lsq.SetParIsFixed(gpRefParTypeObjCryst,
true);
7001 lsq.SetParIsFixed(gpRefParTypeScattDataScale,
false);
7002 std::list<RefinablePar*> vnewpar;
7003 std::list<const RefParType*> vnewpartype;
7005 if(s.number()==1) vnewpar.push_back(&lsq.GetCompiledRefinedObj().GetPar(
"Zero"));
7006 lsq.SetParIsFixed(gpRefParTypeUnitCell,
false);
7007 lsq.SafeRefine(vnewpar, vnewpartype,1.01,2,
true,
true);
7009 TAU_PROFILE_STOP(timer2);
7012 TAU_PROFILE_START(timer1);
7013 vnewpar.push_back(&lsq.GetCompiledRefinedObj().GetPar(
"2ThetaDispl"));
7014 vnewpar.push_back(&lsq.GetCompiledRefinedObj().GetPar(
"2ThetaTransp"));
7015 lsq.SafeRefine(vnewpar, vnewpartype,1.01,2,
true,
true);
7017 lsq.SetParIsFixed(gpRefParTypeScattDataBackground,
false);
7019 const unsigned int nbcomp= mpDiff->GetParentPowderPattern().GetNbPowderPatternComponent();
7020 for(
unsigned int i=0;i<nbcomp;++i)
7021 if(mpDiff->GetParentPowderPattern().GetPowderPatternComponent(i).GetClassName()==
"PowderPatternBackground")
7023 PowderPatternBackground *pback=
dynamic_cast<PowderPatternBackground *
>
7024 (&(mpDiff->GetParentPowderPattern().GetPowderPatternComponent(i)));
7025 pback->FixParametersBeyondMaxresolution(lsq.GetCompiledRefinedObj());
7027 for(
unsigned int i=0; i<lsq.GetCompiledRefinedObj().GetNbPar();i++)
7028 if( (lsq.GetCompiledRefinedObj().GetPar(i).IsFixed()==
false)
7030 vnewpar.push_back(&lsq.GetCompiledRefinedObj().GetPar(i));
7031 lsq.SafeRefine(vnewpar,vnewpartype,1.01,2,
true,
true);
7033 TAU_PROFILE_STOP(timer1);
7036 mpDiff->SetExtractionMode(
true,
true);
7037 mpDiff->ExtractLeBail(5);
7038 TAU_PROFILE_START(timer3);
7039 lsq.SafeRefine(vnewpar,vnewpartype,1.01,3,
true,
true);
7040 TAU_PROFILE_STOP(timer3);
7044 if(update_display) mpDiff->GetParentPowderPattern().UpdateDisplay();
7045 const REAL rw=mpDiff->GetParentPowderPattern().GetRw()*100;
7046 const REAL gof=mpDiff->GetParentPowderPattern().GetChi2()/mpDiff->GetParentPowderPattern().GetNbPointUsed();
7047 if(verbose) cout << boost::format(
" (cycle #%u)\n Rwp=%5.2f%%\n GoF=%9.2f") % j % rw % gof <<endl;
7049 const REAL rw=mpDiff->GetParentPowderPattern().GetRw()*100;
7050 const REAL gof=mpDiff->GetParentPowderPattern().GetChi2()/nbfreepar;
7051 const REAL ngof = this->GetP1IntegratedGoF();
7052 unsigned int nbextinct446=0;
7054 for(
unsigned int i=6;i<fgp.size();++i) nbextinct446+=(
unsigned int)(fgp[i]);
7055 const unsigned int nbrefl = mpDiff->GetNbReflBelowMaxSinThetaOvLambda();
7056 if(verbose>0) cout << boost::format(
" Rwp= %5.2f%% GoF=%9.2f nGoF =%9.2f (%3u reflections, %3u extinct)") % rw % gof % ngof % nbrefl % nbextinct446<<endl;
7060 mpDiff->GetCrystal().Init(a,b,c,d,e,f,spghm,name);
7061 lsq.GetCompiledRefinedObj().RestoreParamSet(saved_par);
7064 mpDiff->GetCrystal().UpdateDisplay();
7065 this->mpDiff->GetParentPowderPattern().UpdateDisplay();
7069 return SPGScore(hm.c_str(),rw,gof,nbextinct446, ngof, nbrefl);
7072 void SpaceGroupExplorer::RunAll(
const bool fitprofile_all,
const bool verbose,
const bool keep_best,
7073 const bool update_display,
const bool fitprofile_p1,
7074 const REAL relative_length_tolerance,
const REAL absolute_angle_tolerance_degree)
7076 Crystal *pCrystal=&(mpDiff->GetCrystal());
7079 const REAL a=pCrystal->GetLatticePar(0),
7080 b=pCrystal->GetLatticePar(1),
7081 c=pCrystal->GetLatticePar(2),
7082 d=pCrystal->GetLatticePar(3),
7083 e=pCrystal->GetLatticePar(4),
7084 f=pCrystal->GetLatticePar(5);
7085 const cctbx::uctbx::unit_cell uc(scitbx::af::double6(a,b,c,d*RAD2DEG,e*RAD2DEG,f*RAD2DEG));
7086 const string spghm=pCrystal->GetSpaceGroup().GetCCTbxSpg().match_tabulated_settings().hermann_mauguin();
7087 const string name=pCrystal->GetName();
7089 cctbx::sgtbx::space_group_symbol_iterator it=cctbx::sgtbx::space_group_symbol_iterator();
7091 unsigned int nbspg=0;
7094 cctbx::sgtbx::space_group_symbols s=it.next();
7095 if(s.number()==0)
break;
7096 cctbx::sgtbx::space_group spg(s);
7097 if(spg.is_compatible_unit_cell(uc,relative_length_tolerance, absolute_angle_tolerance_degree)) nbspg++;
7100 if(verbose) cout << boost::format(
"Beginning spacegroup exploration... %u to go...\n") % nbspg;
7103 mvSPGExtinctionFingerprint.clear();
7106 unsigned int nb_refl_p1=1;
7108 it=cctbx::sgtbx::space_group_symbol_iterator();
7113 cctbx::sgtbx::space_group_symbols s=it.next();
7114 if(s.number()==0)
break;
7115 cctbx::sgtbx::space_group spg(s);
7116 bool compat=spg.is_compatible_unit_cell(uc,relative_length_tolerance,absolute_angle_tolerance_degree);
7120 const string hm=s.universal_hermann_mauguin();
7122 pCrystal->Init(a,b,c,d,e,f,hm,name);
7125 std::map<std::vector<bool>,SPGScore>::iterator posfgp=mvSPGExtinctionFingerprint.find(fgp);
7126 if(posfgp!=mvSPGExtinctionFingerprint.end())
7128 pCrystal->Init(a,b,c,d,e,f,hm,name);
7129 mpDiff->SetExtractionMode(
true,
true);
7130 unsigned int nbrefl = mpDiff->GetNbReflBelowMaxSinThetaOvLambda();
7131 REAL ngof = (posfgp->second.ngof * nbrefl) / posfgp->second.nbreflused;
7132 mvSPG.push_back(SPGScore(hm.c_str(),posfgp->second.rw,posfgp->second.gof,posfgp->second.nbextinct446, ngof, nbrefl));
7133 if(verbose) cout<<boost::format(
" (#%3d) %-14s: Rwp= %5.2f%% GoF=%9.2f nGoF=%9.2f (%3u reflections, %3u extinct)")
7134 % s.number() % hm.c_str() % mvSPG.back().rw % mvSPG.back().gof % mvSPG.back().ngof % mvSPG.back().nbreflused % mvSPG.back().nbextinct446
7135 <<
" [same extinctions as:"<<posfgp->second.hm<<
"]\n";
7139 if(((s.number()==1) && fitprofile_p1) || fitprofile_all) mvSPG.push_back(this->Run(spg,
true,
false,
false, update_display,
7140 relative_length_tolerance, absolute_angle_tolerance_degree));
7141 else mvSPG.push_back(this->Run(spg,
false,
false,
true, update_display,relative_length_tolerance,absolute_angle_tolerance_degree));
7142 if(s.number() == 1) nb_refl_p1 = mvSPG.back().nbreflused;
7143 mvSPG.back().ngof *= mpDiff->GetNbReflBelowMaxSinThetaOvLambda() / (float)nb_refl_p1;
7144 mvSPGExtinctionFingerprint.insert(make_pair(fgp, mvSPG.back()));
7146 if(verbose) cout<<boost::format(
" (#%3d) %-14s: Rwp= %5.2f%% GoF=%9.2f nGoF=%9.2f (%3u reflections, %3u extinct)\n")
7147 % s.number() % hm.c_str() % mvSPG.back().rw % mvSPG.back().gof % mvSPG.back().ngof % mvSPG.back().nbreflused % mvSPG.back().nbextinct446;
7151 mvSPG.sort(compareSPGScore);
7154 if(verbose) cout<<
"Restoring best spacegroup: "<<mvSPG.front().hm<<endl;
7155 pCrystal->ChangeSpaceGroup(mvSPG.front().hm);
7160 pCrystal->Init(a,b,c,d,e,f,spghm,name);
7162 mpDiff->SetExtractionMode(
true,
true);
7163 mpDiff->ExtractLeBail(5);
7166 pCrystal->UpdateDisplay();
7167 this->mpDiff->GetParentPowderPattern().UpdateDisplay();
7172 const list<SPGScore>& SpaceGroupExplorer::GetScores()
const
7177 REAL SpaceGroupExplorer::GetP1IntegratedGoF()
7179 if(mpDiff->GetCrystal().GetSpaceGroup().GetSpaceGroupNumber()==1)
7181 mpDiff->GetPowderPatternIntegratedCalc();
7182 mP1IntegratedProfileMin = mpDiff->GetParentPowderPattern().GetIntegratedProfileMin();
7183 mP1IntegratedProfileMax = mpDiff->GetParentPowderPattern().GetIntegratedProfileMax();
7187 else if (mP1IntegratedProfileMin.size()==0)
return 0;
7190 REAL integratedChi2=0.;
7191 REAL integratedChi2LikeNorm=0.;
7192 const REAL * RESTRICT p1, * RESTRICT p2, * RESTRICT p3;
7193 CrystVector_REAL
const* pcalc = &(mpDiff->GetParentPowderPattern().GetPowderPatternCalc());
7194 CrystVector_REAL
const* pobs = &(mpDiff->GetParentPowderPattern().GetPowderPatternObs());
7195 CrystVector_REAL
const* psigma = &(mpDiff->GetParentPowderPattern().GetPowderPatternObsSigma());
7196 const unsigned int jmax = mpDiff->GetParentPowderPattern().GetNbPointUsed();
7198 unsigned int nbpoint = 0;
7199 for(
unsigned long i=0;i<mP1IntegratedProfileMin.size();i++)
7201 if(mP1IntegratedProfileMin(i) > jmax)
break;
7202 if(mP1IntegratedProfileMax(i) < 0)
continue;
7204 for(
unsigned long j=mP1IntegratedProfileMin(i); j<=mP1IntegratedProfileMax(i); j++)
7207 if(j >= jmax)
break;
7211 v += (*psigma)(j)*(*psigma)(j);
7213 if(v>0) chi2 += (c-o)*(c-o)/v;
7215 return chi2 / nbpoint;
The namespace which includes all objects (crystallographic and algorithmic) in ObjCryst++.
ReflectionProfileType
Profile type for powder (could it be used fopr single crystals on 2D detectors ?)
const RefParType * gpRefParTypeScattDataScale
Type for scattering data scale factors.
std::vector< bool > spgExtinctionFingerprint(Crystal &c, const cctbx::sgtbx::space_group &spg)
Function which determines the unique extinction fingerprint of a spacegroup, by calculating all prese...
ObjRegistry< PowderPatternComponent > gPowderPatternComponentRegistry("List of all PowderPattern Components")
Global registry for all PowderPatternComponent objects.
RadiationType
Type of radiation used.
const RefParType * gpRefParTypeScattDataCorrIntAbsorp
Parameter type for absorption correction.
const RefParType * gpRefParTypeScattDataCorrPos
Parameter type for correction to peak positions.
const RefParType * gpRefParTypeScattDataProfile
Type for reflection profile.
bool ISNAN_OR_INF(REAL r)
Test if the value is a NaN.
const RefParType * gpRefParTypeScattDataBackground
Parameter type for background intensity.
const RefParType * gpRefParTypeObjCryst
Top RefParType for the ObjCryst++ library.
ObjRegistry< RefinableObj > gTopRefinableObjRegistry("Global Top RefinableObj registry")
This is a special registry for 'top' object for an optimization.
Main CIF class - parses the stream and separates data blocks, comments, items, loops.
std::map< std::string, CIFData > mvData
The data blocks, after parsing. The key is the name of the data block.
Crystal class: Unit cell, spacegroup, scatterers.
DiffractionData object for Single Crystal analysis.
Exception class for ObjCryst++ library.
Class to store positions of observed reflections.
vector< hkl > & GetPeakList()
Get peak list.
void AddPeak(const float d, const float iobs=1.0, const float dobssigma=0.0, const float iobssigma=0.0, const int h=0, const int k=0, const int l=0, const float d2calc=0)
Add one peak.
Molecule : class for complex scatterer descriptions using cartesian coordinates with bond length/angl...
virtual const string & GetName() const
Get the name of this object.
virtual void CalcCorr() const
Do the computation of corrected intensities.
CylinderAbsCorr(const PowderPatternDiffraction &data)
Constructor.
virtual const string & GetClassName() const
Get the name of the class.
Generic class to compute components (eg the contribution of a given phase, or background) of a powder...
CrystVector_REAL mPowderPatternCalcVariance
The variance associated to each point of the calculated powder pattern.
virtual const string & GetClassName() const
Name for this class ("RefinableObj", "Crystal",...).
bool IsScalable() const
Is this component scalable ?
virtual const CrystVector_REAL & GetPowderPatternCalc() const =0
Get the calculated powder pattern for this component.
RefinableObjClock mClockPowderPatternIntegratedVarianceCalc
When was the 'integrated' powder pattern variance last computed ?
RefinableObjClock mClockPowderPatternVarianceCalc
When was the powder pattern variance last computed ?
virtual void SetParentPowderPattern(PowderPattern &)=0
Set the PowderPattern object which uses this component.
RefinableObjClock mClockBraggLimits
Get last time the Bragg Limits were changed.
const PowderPattern & GetParentPowderPattern() const
Get the PowderPattern object which uses this component.
virtual void CalcPowderPatternIntegrated_FullDeriv(std::set< RefinablePar * > &vPar)
Calc the integrated powder pattern.
CrystVector_REAL mPowderPatternIntegratedCalc
The calculated powder pattern, integrated.
bool mIsScalable
Scalable ? (crystal phase = scalable, background= not scalable)
RefinableObjClock mClockPowderPatternCalc
When was the powder pattern last computed ?
list< pair< const REAL,const string > > mvLabel
The labels associated to different points of the pattern.
virtual pair< const CrystVector_REAL *, const RefinableObjClock * > GetPowderPatternIntegratedCalc() const =0
Get the integrated values of the powder pattern.
CrystVector_long mIntegratedReflLimits
Interval limits around each reflection, for integrated R-factors.
const RefinableObjClock & GetClockPowderPatternCalc() const
Last time the powder pattern was calculated.
PowderPattern * mpParentPowderPattern
The PowderPattern object in which this component is included.
RefinableObjClock mClockPowderPatternIntegratedCalc
When was the 'integrated' powder pattern last computed ?
CrystVector_REAL mPowderPatternCalc
The calculated component of a powder pattern.
CrystVector_REAL mPowderPatternIntegratedCalcVariance
The variance associated to each point of the calculated powder pattern, integrated.
const RefinableObjClock & GetClockBraggLimits() const
Get last time the Bragg Limits were changed.
const list< pair< const REAL,const string > > & GetPatternLabelList() const
Get a list of labels for the pattern (usually reflection indexes).
Phase to compute a background contribution to a powder pattern using an interpolation.
virtual void CalcPowderPatternIntegrated_FullDeriv(std::set< RefinablePar * > &vPar)
Calc the integrated powder pattern.
void OptimizeBayesianBackground()
Optimize the background using a Bayesian approach.
virtual bool HasPowderPatternCalcVariance() const
Does this component have a variance associated with each calculated point ? i.e., do we use maximum l...
CrystVector_long mPointOrder
Subscript of the points, sorted the correct order, taking into account the type of radiation (monochr...
CrystVector_REAL mvSplinePixel
Vector of pixel values between each interval, for faster CubicSpline calculations.
int mBackgroundNbPoint
Number of fitting points for background.
virtual void TagNewBestConfig() const
During a global optimization, tells the object that the current config is the latest "best" config.
virtual void SetMaxSinThetaOvLambda(const REAL max)
Set the maximum value for sin(theta)/lambda.
virtual void BeginOptimization(const bool allowApproximations=false, const bool enableRestraints=false)
This should be called by any optimization class at the begining of an optimization.
CubicSpline mvSpline
Spline used for interpolation.
CrystVector_REAL mBackgroundInterpPointX
Vector of 2theta values for the fitting points of the background.
virtual void GetGeneGroup(const RefinableObj &obj, CrystVector_uint &groupIndex, unsigned int &firstGroup) const
Get the gene group assigned to each parameter.
virtual void CalcPowderPattern() const
Calc the powder pattern.
virtual void SetParentPowderPattern(PowderPattern &)
Set the PowderPattern object which uses this component.
virtual const CrystVector_REAL & GetPowderPatternCalcVariance() const
Get the variance associated to each point of the calculated powder pattern, for this component.
void FixParametersBeyondMaxresolution(RefinableObj &obj)
Fix parameters corresponding to points of the pattern that are not actually calculated.
virtual const string & GetClassName() const
Name for this class ("RefinableObj", "Crystal",...).
virtual pair< const CrystVector_REAL *, const RefinableObjClock * > GetPowderPatternIntegratedCalcVariance() const
Get the variance associated to each point of the calculated powder pattern, for this component (integ...
RefObjOpt mInterpolationModel
Type of interpolation performed: linear or cubic spline.
CrystVector_REAL mBackgroundInterpPointIntensity
Values of background at interpolating points.
REAL mMaxSinThetaOvLambda
Maximum sin(theta)/lambda for all calculations (10 by default).
virtual pair< const CrystVector_REAL *, const RefinableObjClock * > GetPowderPatternIntegratedCalc() const
Get the integrated values of the powder pattern.
RefinableObjClock mClockBackgroundPoint
Modification of the interpolated points.
void ImportUserBackground(const string &filename)
Import background points from a file (with two columns 2theta (or tof), intensity)
virtual const CrystVector_long & GetBraggLimits() const
Get the pixel positions separating the integration intervals around reflections.
virtual const CrystVector_REAL & GetPowderPatternCalc() const
Get the calculated powder pattern for this component.
RefinableObjClock mClockSpline
Initialization of the spline.
REAL mModelVariance
Constant error (sigma) on the calculated pattern, due to an incomplete model.
Class to compute the contribution to a powder pattern from a crystalline phase.
void SetExtractionMode(const bool extract=true, const bool init=false)
Prepare intensity extraction (Le Bail or Pawley)
const ReflectionProfile & GetProfile() const
Get reflection profile.
virtual pair< const CrystVector_REAL *, const RefinableObjClock * > GetPowderPatternIntegratedCalcVariance() const
Get the variance associated to each point of the calculated powder pattern, for this component (integ...
unsigned int GetProfileFitNetNbObs() const
Get the 'net' number of observed intensities, minus the number of reflections, for a profile fit.
void SetFrozenLatticePar(const unsigned int i, REAL v)
Change one parameter in mFrozenLatticePar. This triggers a call to CalcLocalBMatrix() if the paramete...
REAL GetFrozenLatticePar(const unsigned int i) const
Access to one parameter in mFrozenLatticePar.
virtual const CrystVector_REAL & GetPowderPatternCalc() const
Get the calculated powder pattern for this component.
virtual pair< const CrystVector_REAL *, const RefinableObjClock * > GetPowderPatternIntegratedCalc() const
Get the integrated values of the powder pattern.
const CrystVector_REAL & GetFhklObsSq() const
Get the extracted structure factors modulus (squared), e.g.
virtual void CalcPowderPatternIntegrated_FullDeriv(std::set< RefinablePar * > &vPar)
Calc the integrated powder pattern.
virtual void CalcPowderPattern() const
Calc the powder pattern.
virtual void SetCrystal(Crystal &crystal)
Set the crystal for this experiment.
void CalcPowderReflProfile_FullDeriv(std::set< RefinablePar * > &vPar)
virtual void SetMaxSinThetaOvLambda(const REAL max)
Set the maximum value for sin(theta)/lambda.
virtual bool HasPowderPatternCalcVariance() const
Does this component have a variance associated with each calculated point ? i.e., do we use maximum l...
virtual const CrystVector_REAL & GetPowderPatternCalcVariance() const
Get the variance associated to each point of the calculated powder pattern, for this component.
virtual void GetGeneGroup(const RefinableObj &obj, CrystVector_uint &groupIndex, unsigned int &firstGroup) const
Get the gene group assigned to each parameter.
void CalcIntensityCorr() const
void SetReflectionProfilePar(const ReflectionProfileType prof, const REAL fwhmCagliotiW, const REAL fwhmCagliotiU=0, const REAL fwhmCagliotiV=0, const REAL eta0=0.5, const REAL eta1=0.)
Set reflection profile parameters.
virtual long GetNbReflBelowMaxSinThetaOvLambda() const
Recalc, and get the number of reflections which should be actually used, due to the maximuml sin(thet...
virtual void CalcIhkl() const
bool GetExtractionMode() const
Return true if in extraction mode, i.e. using extracted intensities instead of computed structure fac...
virtual const CrystVector_long & GetBraggLimits() const
Get the pixel positions separating the integration intervals around reflections.
void ExtractLeBail(unsigned int nbcycle=1)
Extract intensities using Le Bail method.
virtual const string & GetClassName() const
Name for this class ("RefinableObj", "Crystal",...).
virtual void EndOptimization()
This should be called by any optimization class at the end of an optimization.
virtual void SetParentPowderPattern(PowderPattern &)
Set the PowderPattern object which uses this component.
void CalcPowderReflProfile() const
void SetProfile(ReflectionProfile *prof)
Assign a new profile.
virtual PowderPatternDiffraction * CreateCopy() const
So-called virtual copy constructor.
bool HasFhklObsSq() const
Return true if there are extracted (le Bail) squared structure factors, false otherwise.
virtual const CrystMatrix_REAL & GetBMatrix() const
This can use either locally stored lattice parameters from mLocalLatticePar, or the Crystal's,...
void CalcFrozenBMatrix() const
Calculate the local BMatrix, used if mFreezeLatticePar is true.
virtual void SetApproximationFlag(const bool allow)
Enable or disable numerical approximations.
bool FreezeLatticePar() const
Do we use local cell parameters ? (see mFrozenLatticePar)
virtual const Radiation & GetRadiation() const
Get the radiation object for this data.
Powder pattern class, with an observed pattern and several calculated components to modelize the patt...
REAL X2Pixel(const REAL x) const
Get the pixel number on the experimental pattern, corresponding to a given (experimental) x coordinat...
void AddPowderPatternComponent(PowderPatternComponent &)
Add a component (phase, backround) to this pattern.
REAL X2STOL(const REAL x) const
Convert X (either 2theta or TOF) to sin(theta)/lambda, depending on the type of radiation.
CrystVector_REAL mPowderPatternCalc
The calculated powder pattern.
Radiation mRadiation
The Radiation corresponding to this experiment.
REAL STOL2X(const REAL stol) const
Convert sin(theta)/lambda to X (i.e.
virtual const CrystVector_REAL & GetLSQObs(const unsigned int) const
Get the observed values for the LSQ function.
const RefinableObjClock & GetClockNbPointUsed() const
Clock corresponding to the last time the number of points used was changed.
CrystVector_REAL mChi2Cumul
The cumulative Chi^2 (integrated or not, depending on the option)
void PrepareIntegratedRfactor() const
Prepare the calculation of the integrated R-factors.
REAL GetIntegratedChi2() const
Return integrated Chi^2.
REAL GetWavelength() const
wavelength of the experiment (in Angstroems)
CrystVector_REAL mPowderPatternBackgroundIntegratedCalc
The calculated powder pattern part which corresponds to 'background' (eg non-scalable components),...
unsigned long mNbPoint
Number of points in the pattern.
REAL GetPowderPatternXMax() const
Get the maximum 2theta.
void ImportPowderPattern2ThetaObs(const string &fileName, const int nbSkip=0)
Import file with 2 columns 2Theta Iobs.
void ImportPowderPatternPSI_DMC(const string &filename)
Import powder pattern, format DMC from PSI.
const CrystVector_long & GetIntegratedProfileMin() const
Get the list of first pixels for the integration intervals.
void SetSigmaToSqrtIobs()
Set sigma=sqrt(Iobs)
void CalcPowderPatternIntegrated() const
Calc the integrated powder pattern.
void ImportPowderPatternFullprof4(const string &fileName)
Import diffraction data from a file, with the first line has 2ThetaMin, step, 2thetaMax,...
REAL GetR() const
Unweighted R-factor.
virtual const CrystVector_REAL & GetLSQWeight(const unsigned int) const
Get the weight values for the LSQ function.
const RefinableObjClock & GetClockPowderPatternPar() const
When were the pattern parameters (2theta range, step) changed ?
const RefinableObjClock & GetIntegratedProfileLimitsClock() const
When were the integration intervals last changed ?
void SetWeightPolynomial(const REAL a, const REAL b, const REAL c, const REAL minRelatIobs=1e-3)
Set w = 1/(a+ Iobs + b*Iobs^2+c*Iobs^3)
const PowderPatternComponent & GetPowderPatternComponent(const string &name) const
Access to a component of the powder pattern.
const CrystVector_REAL & GetPowderPatternWeight() const
Get the weight for each point of the powder pattern.
void CalcNbPointUsed() const
Calculate the number of points of the pattern actually used, from the maximum value of sin(theta)/lam...
void ImportPowderPatternFullprof(const string &fullprofFileName)
Import fullprof-style diffraction data.
virtual std::map< RefinablePar *, CrystVector_REAL > & GetLSQ_FullDeriv(const unsigned int, std::set< RefinablePar * > &vPar)
Get the first derivative for the LSQ function for each parameter supplied in a list.
virtual void SetMaxSinThetaOvLambda(const REAL max)
Set the maximum value for sin(theta)/lambda.
CrystVector_REAL mPowderPatternUsedCalc
The calculated powder pattern. Cropped to the maximum sin(theta)/lambda for LSQ.
void ImportPowderPatternXdd(const string &fileName)
Import *.xdd diffraction data (Topas,...).
virtual unsigned int GetNbLSQFunction() const
Number of LSQ functions.
const CrystVector_REAL & GetPowderPatternVariance() const
Get the variance (obs+model) for each point of the powder pattern.
void ImportPowderPatternCIF(const CIF &cif)
Import CIF powder pattern data.
CrystVector_REAL mPowderPatternUsedObs
The calculated powder pattern. Cropped to the maximum sin(theta)/lambda for LSQ.
REAL GetMuR() const
Get the $\mu R$ value.
void SetWeightToUnit()
Set w = 1.
REAL mMaxSinThetaOvLambda
Displacement correction : REAL m2ThetaDisplacement; Transparency correction : REAL m2ThetaTranspare...
void SetXZero(const REAL newZero)
Change Zero in x (2theta,tof)
void ExportFullprof(const std::string &prefix) const
Export powder pattern & crystal structure in Fullprof format.
void SetScaleFactor(const int i, REAL s)
Access to the scale factor of components (will be 1 for background components)
void SetPowderPatternX(const CrystVector_REAL &x)
Set the x coordinate of the powder pattern : either the 2theta or time-of-flight values for each reco...
CrystVector_REAL mExcludedRegionMinX
Min coordinate for for all excluded regions.
CrystVector_REAL mPowderPatternWeight
The weight for each point of the pattern.
unsigned long GetNbPoint() const
Number of points ?
REAL GetPowderPatternXStep() const
Get the average step in 2theta.
REAL GetRw() const
Get the weighted R-factor.
void Set2ThetaTransparency(const REAL transparency)
Change transparency correction .
virtual const CrystVector_REAL & GetLSQCalc(const unsigned int) const
Get the current calculated value for the LSQ function.
CrystVector_REAL mExcludedRegionMaxX
Max coordinate for 2theta for all excluded regions.
void ImportPowderPatternGSAS(const string &fileName)
Import GSAS standard powder pattern data (see GSAS manual).
void SetRadiation(const Radiation &radiation)
Set the radiation.
CrystVector_REAL mPowderPatternVarianceIntegrated
The complete variance associated to each point of the powder pattern, taking into account observation...
RadiationType GetRadiationType() const
Neutron or x-ray experiment ?
void RemovePowderPatternComponent(PowderPatternComponent &)
Remove a powder pattern component.
void SetPowderPatternObs(const CrystVector_REAL &obs)
Set observed powder pattern from vector array.
REAL X2PixelCorr(const REAL x) const
Get the pixel number on the experimental pattern, from the theoretical (uncorrected) x coordinate,...
const Radiation & GetRadiation() const
Neutron or x-ray experiment ?
RefinableObjClock mClockPowderPatternIntegratedCalc
When was the powder pattern (integrated) last computed ?
void ImportPowderPatternILL_D1A5(const string &filename)
Import powder pattern, format from ILL D1A/D2B (format without counter info)
const RefinableObjClock & GetClockPowderPatternAbsCorr() const
When were the absorption correction parameters (muR) last changed ?
void CalcPowderPattern() const
Calc the powder pattern.
const CrystVector_REAL & GetPowderPatternX() const
Get the vector of X (2theta or time-of-flight) coordinates.
unsigned long GetNbPointUsed() const
Number of points actually calculated (below the chosen max(sin(theta)/lambda)) ?
void AddExcludedRegion(const REAL min2Theta, const REAL max2theta)
Add an Exclusion region, in 2theta, which will be ignored when computing R's XMLInput values must be,...
virtual void GlobalOptRandomMove(const REAL mutationAmplitude, const RefParType *type=gpRefParTypeObjCryst)
Make a random move of the current configuration.
void ImportPowderPatternTOF_ISIS_XYSigma(const string &fileName)
Import TOF file (ISIS type, 3 columns t, Iobs, sigma(Iobs))
const CrystVector_REAL & GetChi2Cumul(const int mode=-1) const
Get the powder pattern cumulative Chi^2.
CrystVector_REAL mPowderPatternBackgroundCalc
The calculated powder pattern part which corresponds to 'background' (eg non-scalable components).
CrystVector_REAL mPowderPatternObs
The observed powder pattern.
virtual void GetGeneGroup(const RefinableObj &obj, CrystVector_uint &groupIndex, unsigned int &firstGroup) const
Get the gene group assigned to each parameter.
RefinableObjClock mClockPowderPatternCalc
When was the powder pattern last computed ?
CrystVector_REAL mPowderPatternObsSigma
The sigma of the observed pattern.
virtual REAL GetLogLikelihood() const
Get -log(likelihood) of the current configuration for the object.
CrystVector_REAL mPowderPatternUsedWeight
The weight for each point of the pattern. Cropped to the maximum sin(theta)/lambda for LSQ.
REAL STOL2Pixel(const REAL stol) const
Convert sin(theta)/lambda to pixel, depending on the type of radiation.
REAL GetPowderPatternXMin() const
Get the Minimum 2theta.
virtual void InitOptions()
Initialize options.
virtual void Init()
Init parameters and options.
void FitScaleFactorForRw() const
Fit the scale(s) factor of each component to minimize Rw.
REAL GetChi2() const
Return conventionnal Chi^2.
REAL GetMaxSinThetaOvLambda() const
Get the maximum value for sin(theta)/lambda.
RefinableObjClock mClockScaleFactor
Last modification of the scale factor.
CrystVector_REAL mPowderPatternIntegratedCalc
The calculated powder pattern, integrated.
void ImportPowderPatternMultiDetectorLLBG42(const string &fileName)
diffraction data in a multi-detector format (fullprof format #6).
void PrintObsCalcData(ostream &os=cout) const
Print to thee screen/console the observed and calculated pattern (long, mostly useful for debugging)
const CrystVector_REAL & GetPowderPatternObs() const
Get the observed powder pattern.
CrystVector_REAL mX
Vector of x coordinates (either 2theta or time-of-flight) for the pattern.
unsigned int GetNbPowderPatternComponent() const
Number of components.
CrystVector_REAL mPowderPatternVariance
The complete variance associated to each point of the powder pattern, taking into account observation...
RefinableObjClock mClockCorrAbs
Last modification of absorption correction parameters.
const RefinableObjClock & GetClockPowderPatternCalc() const
Last time the pattern was calculated.
REAL X2XCorr(const REAL x) const
Get the experimental x (2theta, tof) from the theoretical value, taking into account all corrections ...
void SetRadiationType(const RadiationType radiation)
Set the radiation type.
PeakList FindPeaks(const float dmin=2.0, const float maxratio=0.01, const unsigned int maxpeak=100)
Find peaks in the pattern.
void ImportPowderPatternSietronicsCPI(const string &fileName)
Import *.cpi Sietronics diffraction data.
void SetPowderPatternPar(const REAL min, const REAL step, unsigned long nbPoint)
\briefSet the powder pattern angular range & resolution parameter.
const CrystVector_long & GetIntegratedProfileMax() const
Get the list of last pixels for the integration intervals.
void SetWeightToInvSigmaSq(const REAL minRelatSigma=1e-3)
Set w = 1/sigma^2.
void FitScaleFactorForR() const
Fit the scale(s) factor of each component to minimize R.
const CrystVector_REAL & GetPowderPatternCalc() const
Get the calculated powder pattern.
virtual void BeginOptimization(const bool allowApproximations=false, const bool enableRestraints=false)
This should be called by any optimization class at the begining of an optimization.
RefinableObjClock mClockPowderPatternXCorr
Corrections to 2Theta.
bool mIsXAscending
Is the mX vector sorted in ascending order ? (true for 2theta, false for TOF)
RefinableObjClock mClockPowderPatternPar
When were the pattern parameters (2theta or time-of-flight range) changed ?
REAL mXZero
Zero correction : Thus mPowderPattern2ThetaMin=(mPowderPattern2ThetaMin-m2ThetaZero)
const RefinableObjClock & GetClockPowderPatternRadiation() const
When were the radiation parameter (radiation type, wavelength) changed ?
void SetMuR(const REAL muR)
Set the $\mu R$ value for cylinder absorption correction.
REAL GetChi2_Option() const
Return the conventionnal or integrated Chi^2, depending on the option.
void SavePowderPattern(const string &filename="powderPattern.out") const
Save powder pattern to one file, text format, 3 columns theta Iobs Icalc.
void SetWavelength(const REAL lambda)
Set the wavelength of the experiment (in Angstroems).
void Set2ThetaDisplacement(const REAL displacement)
Change displacement correction .
void ImportPowderPattern2ThetaObsSigma(const string &fileName, const int nbSkip=0)
Import file with 3 columns 2Theta Iobs Sigma.
void SetPowderPatternObsSigma(const CrystVector_REAL &sigma)
Set sigma for each point of the observed powder pattern.
virtual const string & GetClassName() const
Name for this class ("RefinableObj", "Crystal",...).
const CrystVector_REAL & GetScaleFactor() const
Access the scale factors (see PowderPattern::mScaleFactor)
const RefinableObjClock & GetClockPowderPatternXCorr() const
When were the parameters for 2theta/TOF correction (zero, transparency, displacement) last changed ?
RefinableObjClock mClockPowderPatternRadiation
When were the radiation parameter (radiation type, wavelength) changed ?
const CrystVector_REAL & GetPowderPatternObsSigma() const
Get the sigma for each point of the observed powder pattern.
Structure to hold the score corresponding to a given spacegroup.
SPGScore(const string &s, const REAL r, const REAL g, const unsigned int nbextinct, const REAL ngof=0, const unsigned int nbrefl=0)
Structure to hold the score corresponding to a given spacegroup.
This object is used to estimate the background in a powder pattern, using a Bayesian approach (David ...
Abstract base class for reflection profiles.
Pseudo-Voigt reflection profile.
void SetProfilePar(const REAL fwhmCagliotiW, const REAL fwhmCagliotiU=0, const REAL fwhmCagliotiV=0, const REAL eta0=0.5, const REAL eta1=0.)
Set reflection profile parameters.
Double-Exponential Pseudo-Voigt profile for TOF.
void SetUnitCell(const UnitCell &cell)
Set unit cell.
Base class to compute all kind of corrections to intensities: Lorentz, Polar, absorption,...
const ScatteringData * mpData
The associated ScatteringData object.
CrystVector_REAL mCorr
The vector of correction to intensities.
RefinableObjClock mClockCorrCalc
The clock marking the last time the correction was calculated.
Class to define the radiation (type, monochromaticity, wavelength(s)) of an experiment.
WavelengthType GetWavelengthType() const
Get the Wavelength type (monochromatic, Alpha1+Alpha2, Time Of Flight...)
const CrystVector_REAL & GetWavelength() const
Get the wavelength(s) in Angstroems.
void SetRadiationType(const RadiationType)
Set the radiation type (X-Rays, Neutron)
void SetWavelength(const REAL)
Set the (monochromatic) wavelength of the beam.
void SetWavelengthType(const WavelengthType &type)
Set the Wavelength type (monochromatic, Alpha1+Alpha2, Time Of Flight...)
RadiationType GetRadiationType() const
Get the radiation type (X-Rays, Neutron)
Class to compute structure factors for a set of reflections and a Crystal.
const CrystVector_REAL & GetTheta() const
Return an array with theta values for all reflections.
virtual void SetMaxSinThetaOvLambda(const REAL max)
Set the maximum value for sin(theta)/lambda.
virtual void SetCrystal(Crystal &crystal)
Set the crystal for this experiment.
const RefinableObjClock & GetClockTheta() const
Clock the last time the sin(theta)/lambda and theta arrays were re-computed.
CrystVector_REAL mX
reflection coordinates in an orthonormal base
virtual const CrystMatrix_REAL & GetBMatrix() const
Get access to the B matrix used to compute reflection positions.
long GetNbRefl() const
Return the number of reflections in this experiment.
virtual long GetNbReflBelowMaxSinThetaOvLambda() const
Recalc, and get the number of reflections which should be actually used, due to the maximuml sin(thet...
Abstract Base Class to describe the scattering power of any Scatterer component in a crystal.
list of scattering positions in a crystal, associated with the corresponding occupancy and a pointer ...
long GetNbComponent() const
Number of components.
Unit Cell class: Unit cell with spacegroup information.
Vector library (Blitz++ mimic) for ObjCryst++.
void Init(const CrystVector_REAL &x, const CrystVector_REAL &y, const REAL yp1, const REAL ypn)
Spline with given extremum derivatives.
void AddRefinableObj(RefinableObj &)
Add a refined object. All sub-objects are also added.
virtual REAL GetLogLikelihood() const
The optimized (minimized, actually) function.
(Quick & dirty) Least-Squares Refinement Object with Numerical derivatives
void Refine(int nbCycle=1, bool useLevenbergMarquardt=false, const bool silent=false, const bool callBeginEndOptimization=true, const float minChi2var=0.01)
Do the refinement.
void SetParIsUsed(const std::string &parName, const bool use)
Set a parameter to be used.
void SetRefinedObj(RefinableObj &obj, const unsigned int LSQFuncIndex=0, const bool init=true, const bool recursive=false)
Choose the object to refine.
void PrepareRefParList(const bool copy_param=false)
Prepare the full parameter list for the refinement.
class of refinable parameter types.
bool IsDescendantFromOrSameAs(const RefParType *type) const
Returns true if the parameter is a descendant of 'type'.
We need to record exactly when refinable objects have been modified for the last time (to avoid re-co...
void Reset()
Reset a Clock to 0, to force an update.
void AddChild(const RefinableObjClock &)
Add a 'child' clock.
void RemoveChild(const RefinableObjClock &)
remove a child clock. This also tells the child clock to remove the parent.
void Click()
Record an event for this clock (generally, the 'time' an object has been modified,...
Generic class for parameters of refinable objects.
void SetDerivStep(const REAL)
Fixed step to use to compute numerical derivative.
void SetGlobalOptimStep(const REAL)
Maximum step to use during Global Optimization algorithms.
bool IsUsed() const
Is the parameter used (if not, it is simply irrelevant in the model) ?
void AssignClock(RefinableObjClock &clock)
Generic Refinable Object.
virtual void EndOptimization()
This should be called by any optimization class at the end of an optimization.
void AddPar(const RefinablePar &newRefPar)
Add a refinable parameter.
ObjRegistry< RefinableObj > mSubObjRegistry
Registry of RefinableObject needed for this object (owned by this object or not)
bool IsBeingRefined() const
Is the object being refined ? (Can be refined by one algorithm at a time only.)
virtual void SetName(const string &name)
Name of the object.
virtual const CrystVector_REAL & GetLSQDeriv(const unsigned int, RefinablePar &)
Get the first derivative values for the LSQ function, for a given parameter.
virtual void RegisterClient(RefinableObj &) const
Register a new object using this object.
RefinablePar & GetPar(const long i)
Access all parameters in the order they were inputted.
vector< RefinablePar * >::iterator RemovePar(RefinablePar *refPar)
Remove a refinable parameter.
virtual const string & GetName() const
Name of the object.
long GetNbPar() const
Total number of refinable parameter in the object.
void ResetParList()
Re-init the list of refinable parameters, removing all parameters.
virtual void SetApproximationFlag(const bool allow)
Enable or disable numerical approximations.
virtual void BeginOptimization(const bool allowApproximations=false, const bool enableRestraints=false)
This should be called by any optimization class at the begining of an optimization.
const RefinableObjClock & GetClockMaster() const
This clocks records any change in the object. See refinableObj::mClockMaster.
void AddOption(RefObjOpt *opt)
int mOptimizationDepth
Is the object being refined or optimized ? if mOptimizationDepth=0, no optimization is taking place.
virtual void DeRegisterClient(RefinableObj &) const
Deregister an object (which not any more) using this object.
virtual void UpdateDisplay() const
If there is an interface, this should be automatically be called each time there is a 'new,...
RefinableObjClock mClockMaster
Master clock, which is changed whenever the object has been altered.
string mName
Name for this RefinableObject. Should be unique, at least in the same scope.+.
void AddSubRefObj(RefinableObj &)
void SetGlobalOptimStep(const RefParType *type, const REAL step)
Change the maximum step to use during Global Optimization algorithms.
void RemoveSubRefObj(RefinableObj &)
virtual void GlobalOptRandomMove(const REAL mutationAmplitude, const RefParType *type=gpRefParTypeObjCryst)
Make a random move of the current configuration.
Conjugate Gradient Algorithm object.
virtual void Optimize(long &nbSteps, const bool silent=false, const REAL finalcost=0, const REAL maxTime=-1)
Launch optimization (a single run) for N steps.
Simple chronometer class, with microsecond precision.
output a number as a formatted float:
output one or several vectors as (a) column(s):
Output vectors as column arrays, with the first 3 columns printed as integers.
Abstract base class for all objects in wxCryst.
WX Class for PowderPattern objects.