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);
3867 VFN_DEBUG_MESSAGE(
"PowderPattern::SavePowderPattern",5)
3869 ofstream out(filename.c_str());
3870 CrystVector_REAL ttheta;
3872 if(this->
GetRadiation().GetWavelengthType()!=WAVELENGTH_TOF) ttheta *= RAD2DEG;
3874 CrystVector_REAL diff;
3877 out <<
"# 2Theta/TOF Iobs ICalc Iobs-Icalc Weight Comp0" << endl;
3878 out << FormatVertVector<REAL>(ttheta,
3882 mPowderPatternComponentRegistry.GetObj(0).mPowderPatternCalc,16,8);
3884 VFN_DEBUG_MESSAGE(
"DiffractionDataPowder::SavePowderPattern:End",3)
3889 VFN_DEBUG_MESSAGE(
"DiffractionDataPowder::PrintObsCalcData()",5);
3890 CrystVector_REAL ttheta;
3892 if(this->
GetRadiation().GetWavelengthType()!=WAVELENGTH_TOF) ttheta *= RAD2DEG;
3893 os <<
"PowderPattern : " <<
mName <<endl;
3894 os <<
" 2Theta/TOF Obs Sigma Calc Weight" <<endl;
3908 TAU_PROFILE(
"PowderPattern::GetR()",
"void ()",TAU_DEFAULT);
3913 unsigned long maxPoints=mNbPointUsed;
3914 if( (
true==mStatisticsExcludeBackground)
3917 const REAL *p1, *p2, *p3;
3924 VFN_DEBUG_MESSAGE(
"PowderPattern::GetR():Exclude Backgd",4);
3925 for(
unsigned long i=0;i<maxPoints;i++)
3927 tmp1 += ((*p1)-(*p2)) * ((*p1)-(*p2));
3928 tmp2 += ((*p2)-(*p3)) * ((*p2)-(*p3));
3934 VFN_DEBUG_MESSAGE(
"PowderPattern::GetR():Exclude Backgd,Exclude regions",4);
3935 unsigned long min,max;
3937 for(
int j=0;j<nbExclude;j++)
3941 if(min>maxPoints)
break;
3942 if(max>maxPoints)max=maxPoints;
3945 tmp1 += ((*p1)-(*p2)) * ((*p1)-(*p2));
3946 tmp2 += ((*p2)-(*p3)) * ((*p2)-(*p3));
3954 for(;i<maxPoints;i++)
3956 tmp1 += ((*p1)-(*p2)) * ((*p1)-(*p2));
3957 tmp2 += ((*p2)-(*p3)) * ((*p2)-(*p3));
3965 const REAL *p1, *p2;
3971 VFN_DEBUG_MESSAGE(
"PowderPattern::GetR()",4);
3972 for(
unsigned long i=0;i<maxPoints;i++)
3974 tmp1 += ((*p1)-(*p2))*((*p1)-(*p2));
3975 tmp2 += (*p2) * (*p2);
3982 VFN_DEBUG_MESSAGE(
"PowderPattern::GetR(),Exclude regions",4);
3983 unsigned long min,max;
3985 for(
int j=0;j<nbExclude;j++)
3989 if(min>maxPoints)
break;
3990 if(max>maxPoints)max=maxPoints;
3993 tmp1 += ((*p1)-(*p2))*((*p1)-(*p2));
3994 tmp2 += (*p2) * (*p2);
4001 for(;i<maxPoints;i++)
4003 tmp1 += ((*p1)-(*p2))*((*p1)-(*p2));
4004 tmp2 += (*p2) * (*p2);
4010 VFN_DEBUG_MESSAGE(
"PowderPattern::GetR()="<<sqrt(tmp1/tmp2),4);
4014 return sqrt(tmp1/tmp2);
4017 REAL PowderPattern::GetIntegratedR()
const
4026 VFN_DEBUG_ENTRY(
"PowderPattern::GetIntegratedR()",4);
4027 TAU_PROFILE(
"PowderPattern::GetIntegratedR()",
"void ()",TAU_DEFAULT);
4031 const long numInterval=mIntegratedPatternMin.numElements();
4032 if( (
true==mStatisticsExcludeBackground)
4035 const REAL *p1, *p2, *p3;
4036 CrystVector_REAL integratedCalc(numInterval);
4038 CrystVector_REAL backgdCalc(numInterval);
4040 REAL *pp1=integratedCalc.data();
4041 REAL *pp2=backgdCalc.data();
4042 for(
int i=0;i<numInterval;i++)
4044 const long max=mIntegratedPatternMax(i);
4046 for(
int j=mIntegratedPatternMin(i);j<=max;j++) *pp1 += *p1++;
4049 for(
int j=mIntegratedPatternMin(i);j<=max;j++) *pp2 += *p1++;
4053 p1=integratedCalc.data();
4054 p2=mIntegratedObs.data();
4055 p3=backgdCalc.data();
4056 VFN_DEBUG_MESSAGE(
"PowderPattern::GetIntegratedR():Exclude Backgd",2);
4057 for(
long i=0;i<numInterval;i++)
4059 tmp1 += ((*p1)-(*p2)) * ((*p1)-(*p2));
4060 tmp2 += ((*p2)-(*p3)) * ((*p2)-(*p3));
4066 const REAL *p1, *p2;
4067 CrystVector_REAL integratedCalc(numInterval);
4069 REAL *pp1=integratedCalc.data();
4070 for(
int i=0;i<numInterval;i++)
4072 const long max=mIntegratedPatternMax(i);
4074 for(
int j=mIntegratedPatternMin(i);j<=max;j++) *pp1 += *p1++;
4077 p1=integratedCalc.data();
4078 p2=mIntegratedObs.data();
4079 VFN_DEBUG_MESSAGE(
"PowderPattern::GetIntegratedR()",2);
4080 for(
long i=0;i<numInterval;i++)
4082 tmp1 += ((*p1)-(*p2))*((*p1)-(*p2));
4083 tmp2 += (*p2) * (*p2);
4089 VFN_DEBUG_EXIT(
"PowderPattern::GetIntegratedR()="<<sqrt(tmp1/tmp2),4);
4090 return sqrt(tmp1/tmp2);
4101 TAU_PROFILE(
"PowderPattern::GetRw()",
"void ()",TAU_DEFAULT);
4102 VFN_DEBUG_MESSAGE(
"PowderPattern::GetRw()",3);
4111 unsigned long maxPoints=mNbPointUsed;
4113 if( (
true==mStatisticsExcludeBackground)
4116 VFN_DEBUG_MESSAGE(
"PowderPattern::GetRw():Exclude Backgd",3);
4117 const REAL *p1, *p2, *p3, *p4;
4125 for(
unsigned long i=0;i<maxPoints;i++)
4127 tmp1 += *p4 * ((*p1)-(*p2)) * ((*p1)-(*p2));
4128 tmp2 += *p4++ * ((*p2)-(*p3)) * ((*p2)-(*p3));
4134 unsigned long min,max;
4136 for(
int j=0;j<nbExclude;j++)
4140 if(min>maxPoints)
break;
4141 if(max>maxPoints)max=maxPoints;
4144 tmp1 += *p4 * ((*p1)-(*p2)) * ((*p1)-(*p2));
4145 tmp2 += *p4++ * ((*p2)-(*p3)) * ((*p2)-(*p3));
4154 for(;i<maxPoints;i++)
4156 tmp1 += *p4 * ((*p1)-(*p2)) * ((*p1)-(*p2));
4157 tmp2 += *p4++ * ((*p2)-(*p3)) * ((*p2)-(*p3));
4165 VFN_DEBUG_MESSAGE(
"PowderPattern::GetRw()",3);
4166 const REAL *p1, *p2, *p4;
4173 for(
unsigned long i=0;i<maxPoints;i++)
4175 tmp1 += *p4 * ((*p1)-(*p2))*((*p1)-(*p2));
4176 tmp2 += *p4++ * (*p2) * (*p2);
4182 unsigned long min,max;
4184 for(
int j=0;j<nbExclude;j++)
4188 if(min>maxPoints)
break;
4189 if(max>maxPoints)max=maxPoints;
4192 tmp1 += *p4 * ((*p1)-(*p2))*((*p1)-(*p2));
4193 tmp2 += *p4++ * (*p2) * (*p2);
4201 for(;i<maxPoints;i++)
4203 tmp1 += *p4 * ((*p1)-(*p2))*((*p1)-(*p2));
4204 tmp2 += *p4++ * (*p2) * (*p2);
4209 VFN_DEBUG_MESSAGE(
"PowderPattern::GetRw()="<<sqrt(tmp1/tmp2),3);
4210 return sqrt(tmp1/tmp2);
4212 REAL PowderPattern::GetIntegratedRw()
const
4221 TAU_PROFILE(
"PowderPattern::GetIntegratedRw()",
"void ()",TAU_DEFAULT);
4225 const long numInterval=mIntegratedPatternMin.numElements();
4226 if( (
true==mStatisticsExcludeBackground)
4229 const REAL *p1, *p2, *p3, *p4;
4230 CrystVector_REAL integratedCalc(numInterval);
4232 CrystVector_REAL backgdCalc(numInterval);
4234 REAL *pp1=integratedCalc.data();
4235 REAL *pp2=backgdCalc.data();
4236 for(
int i=0;i<numInterval;i++)
4238 const long max=mIntegratedPatternMax(i);
4240 for(
int j=mIntegratedPatternMin(i);j<=max;j++) *pp1 += *p1++;
4243 for(
int j=mIntegratedPatternMin(i);j<=max;j++) *pp2 += *p1++;
4247 p1=integratedCalc.data();
4248 p2=mIntegratedObs.data();
4249 p3=backgdCalc.data();
4250 if(mIntegratedWeight.numElements()==0) p4=mIntegratedWeightObs.data();
4251 else p4=mIntegratedWeight.data();
4252 VFN_DEBUG_MESSAGE(
"PowderPattern::GetIntegratedRw():Exclude Backgd",4);
4253 for(
long i=0;i<numInterval;i++)
4255 tmp1 += *p4 * ((*p1)-(*p2)) * ((*p1)-(*p2));
4256 tmp2 += *p4++ * ((*p2)-(*p3)) * ((*p2)-(*p3));
4264 const REAL *p1, *p2, *p4;
4265 CrystVector_REAL integratedCalc(numInterval);
4267 REAL *pp1=integratedCalc.data();
4268 for(
int i=0;i<numInterval;i++)
4270 const long max=mIntegratedPatternMax(i);
4272 for(
int j=mIntegratedPatternMin(i);j<=max;j++) *pp1 += *p1++;
4275 p1=integratedCalc.data();
4276 p2=mIntegratedObs.data();
4277 if(mIntegratedWeight.numElements()==0) p4=mIntegratedWeightObs.data();
4278 else p4=mIntegratedWeight.data();
4279 VFN_DEBUG_MESSAGE(
"PowderPattern::GetIntegratedRw()",4);
4280 for(
long i=0;i<numInterval;i++)
4282 tmp1 += *p4 * ((*p1)-(*p2))*((*p1)-(*p2));
4283 tmp2 += *p4++ * (*p2) * (*p2);
4289 VFN_DEBUG_MESSAGE(
"PowderPattern::GetIntegratedRw()="<<sqrt(tmp1/tmp2),4);
4293 return sqrt(tmp1/tmp2);
4314 TAU_PROFILE(
"PowderPattern::GetChi2()",
"void ()",TAU_DEFAULT);
4316 VFN_DEBUG_ENTRY(
"PowderPattern::GetChi2()",3);
4318 const unsigned long maxPoints=mNbPointUsed;
4322 VFN_DEBUG_MESSAGE(
"PowderPattern::GetChi2()Integrated profiles",3);
4323 const REAL * RESTRICT p1, * RESTRICT p2, * RESTRICT p3;
4330 for(
unsigned long i=0;i<maxPoints;i++)
4332 mChi2 += *p3 * ((*p1)-(*p2))*((*p1)-(*p2));
4334 else mChi2LikeNorm -= log(*p3++);
4340 unsigned long min,max;
4342 for(
int j=0;j<nbExclude;j++)
4346 if(min>maxPoints)
break;
4347 if(max>maxPoints)max=maxPoints;
4350 mChi2 += *p3 * ((*p1)-(*p2))*((*p1)-(*p2));
4352 else mChi2LikeNorm -= log(*p3++);
4360 for(;i<maxPoints;i++)
4362 mChi2 += *p3 * ((*p1)-(*p2))*((*p1)-(*p2));
4364 else mChi2LikeNorm -= log(*p3++);
4369 VFN_DEBUG_MESSAGE(
"Chi^2="<<mChi2<<
", log(norm)="<<mChi2LikeNorm,3)
4371 VFN_DEBUG_EXIT(
"PowderPattern::GetChi2()="<<mChi2,3);
4381 return mIntegratedChi2;
4384 if(mClockIntegratedChi2>
mClockMaster)
return mIntegratedChi2;
4392 this->FitScaleFactorForIntegratedRw();
4394 TAU_PROFILE(
"PowderPattern::GetIntegratedChi2()",
"void ()",TAU_DEFAULT);
4396 VFN_DEBUG_ENTRY(
"PowderPattern::GetIntegratedChi2()",3);
4399 mIntegratedChi2LikeNorm=0.;
4400 VFN_DEBUG_MESSAGE(
"PowderPattern::GetIntegratedChi2()",3);
4401 const REAL * RESTRICT p1, * RESTRICT p2, * RESTRICT p3;
4403 p2=mIntegratedObs.data();
4404 if(mIntegratedWeight.numElements()==0) p3=mIntegratedWeightObs.data();
4405 else p3=mIntegratedWeight.data();
4406 double weightProd=1.;
4407 VFN_DEBUG_MESSAGE(
"PowderPattern::GetIntegratedIntegratedRw()",4);
4408 for(
unsigned long i=0;i<mNbIntegrationUsed;)
4412 for(
unsigned long j=0;j<32;++j)
4414 mIntegratedChi2 += *p3 * ((*p1)-(*p2))*((*p1)-(*p2));
4415 if(*p3>0) weightProd *= *p3;
4417 if(++i == mNbIntegrationUsed)
break;
4419 mIntegratedChi2LikeNorm -= log(weightProd);
4422 mIntegratedChi2LikeNorm/=2;
4423 VFN_DEBUG_MESSAGE(
"Chi^2="<<mIntegratedChi2<<
", log(norm)="<<mIntegratedChi2LikeNorm,3)
4424 mClockIntegratedChi2.Click();
4425 VFN_DEBUG_EXIT(
"PowderPattern::GetChi2()="<<mIntegratedChi2,3);
4426 return mIntegratedChi2;
4443 TAU_PROFILE(
"PowderPattern::FitScaleFactorForR()",
"void ()",TAU_DEFAULT);
4444 VFN_DEBUG_ENTRY(
"PowderPattern::FitScaleFactorForR()",3);
4446 mScalableComponentIndex.resize(mPowderPatternComponentRegistry.GetNb());
4448 for(
int i=0;i<mPowderPatternComponentRegistry.GetNb();i++)
4450 if(mPowderPatternComponentRegistry.GetObj(i).IsScalable())
4451 mScalableComponentIndex(nbScale++)=i;
4453 VFN_DEBUG_MESSAGE(
"-> Number of Scale Factors:"<<nbScale<<
":Index:"<<endl<<mScalableComponentIndex,3);
4456 VFN_DEBUG_EXIT(
"PowderPattern::FitScaleFactorForR(): No scalable component!",3);
4459 mScalableComponentIndex.resizeAndPreserve(nbScale);
4461 mFitScaleFactorM.resize(nbScale,nbScale);
4462 mFitScaleFactorB.resize(nbScale,1);
4463 mFitScaleFactorX.resize(nbScale,1);
4468 for(
int i=0;i<nbScale;i++)
4470 for(
int j=i;j<nbScale;j++)
4474 const REAL *p1=mPowderPatternComponentRegistry.GetObj(mScalableComponentIndex(i))
4475 .mPowderPatternCalc.data();
4476 const REAL *p2=mPowderPatternComponentRegistry.GetObj(mScalableComponentIndex(j))
4477 .mPowderPatternCalc.data();
4479 for(
unsigned long k=0;k<mNbPointUsed;k++) m += *p1++ * *p2++;
4480 mFitScaleFactorM(i,j)=m;
4481 mFitScaleFactorM(j,i)=m;
4484 for(
int i=0;i<nbScale;i++)
4487 const REAL *p2=mPowderPatternComponentRegistry.GetObj(mScalableComponentIndex(i))
4488 .mPowderPatternCalc.data();
4491 for(
unsigned long k=0;k<mNbPointUsed;k++) b += *p1++ * *p2++;
4495 for(
unsigned long k=0;k<mNbPointUsed;k++) b += (*p1++ - *p3++) * *p2++;
4497 mFitScaleFactorB(i,0) =b;
4502 unsigned long min,max;
4503 for(
int i=0;i<nbScale;i++)
4505 for(
int j=i;j<nbScale;j++)
4509 const REAL *p1=mPowderPatternComponentRegistry.GetObj(mScalableComponentIndex(i))
4510 .mPowderPatternCalc.data();
4511 const REAL *p2=mPowderPatternComponentRegistry.GetObj(mScalableComponentIndex(j))
4512 .mPowderPatternCalc.data();
4515 for(
int k=0;k<nbExclude;k++)
4519 if(min>mNbPointUsed)
break;
4520 if(max>mNbPointUsed)max=mNbPointUsed;
4522 for(;l<min;l++) m += *p1++ * *p2++;
4527 for(;l<mNbPointUsed;l++) m += *p1++ * *p2++;
4528 mFitScaleFactorM(i,j)=m;
4529 mFitScaleFactorM(j,i)=m;
4532 for(
int i=0;i<nbScale;i++)
4535 const REAL *p2=mPowderPatternComponentRegistry
4536 .GetObj(mScalableComponentIndex(i))
4537 .mPowderPatternCalc.data();
4542 for(
int k=0;k<nbExclude;k++)
4546 if(min>mNbPointUsed)
break;
4547 if(max>mNbPointUsed)max=mNbPointUsed;
4549 for(;l<min;l++) b += *p1++ * *p2++;
4554 for(;l<mNbPointUsed;l++) b += *p1++ * *p2++;
4559 for(
int k=0;k<nbExclude;k++)
4563 if(min>mNbPointUsed)
break;
4564 if(max>mNbPointUsed)max=mNbPointUsed;
4566 for(;l<min;l++) b += (*p1++ - *p3++) * *p2++;
4571 for(;l<mNbPointUsed;l++) b += (*p1++ - *p3++) * *p2++;
4573 mFitScaleFactorB(i,0) =b;
4576 if(1==nbScale) mFitScaleFactorX=mFitScaleFactorB(0)/mFitScaleFactorM(0);
4578 mFitScaleFactorX=product(InvertMatrix(mFitScaleFactorM),mFitScaleFactorB);
4579 VFN_DEBUG_MESSAGE(
"B, M, X"<<endl<<mFitScaleFactorB<<endl<<mFitScaleFactorM<<endl<<mFitScaleFactorX,2)
4580 for(
int i=0;i<nbScale;i++)
4582 const REAL * p1=mPowderPatternComponentRegistry.GetObj(mScalableComponentIndex(i))
4583 .mPowderPatternCalc.data();
4585 const REAL s = mFitScaleFactorX(i)
4586 -mScaleFactor(mScalableComponentIndex(i));
4589 (*fpObjCrystInformUser)(
"Warning:FitScaleFactorForR: working around NaN scale factor...");
4592 for(
unsigned long j=0;j<mNbPointUsed;j++) *p0++ += s * *p1++;
4593 VFN_DEBUG_MESSAGE(
"-> Old:"<<mScaleFactor(mScalableComponentIndex(i)) <<
" Change:"<<mFitScaleFactorX(i),2);
4594 mScaleFactor(mScalableComponentIndex(i)) = mFitScaleFactorX(i);
4598 VFN_DEBUG_EXIT(
"PowderPattern::FitScaleFactorForR():End",3);
4601 void PowderPattern::FitScaleFactorForIntegratedR()
const
4610 VFN_DEBUG_ENTRY(
"PowderPattern::FitScaleFactorForIntegratedR()",3);
4611 TAU_PROFILE(
"PowderPattern::FitScaleFactorForIntegratedR()",
"void ()",TAU_DEFAULT);
4613 mScalableComponentIndex.resize(mPowderPatternComponentRegistry.GetNb());
4615 for(
int i=0;i<mPowderPatternComponentRegistry.GetNb();i++)
4617 if(mPowderPatternComponentRegistry.GetObj(i).IsScalable())
4618 mScalableComponentIndex(nbScale++)=i;
4620 VFN_DEBUG_MESSAGE(
"-> Number of Scale Factors:"<<nbScale<<
":Index:"<<endl<<mScalableComponentIndex,2);
4623 VFN_DEBUG_EXIT(
"PowderPattern::FitScaleFactorForIntegratedR(): No scalable component!",3);
4626 mScalableComponentIndex.resizeAndPreserve(nbScale);
4628 mFitScaleFactorM.resize(nbScale,nbScale);
4629 mFitScaleFactorB.resize(nbScale,1);
4630 mFitScaleFactorX.resize(nbScale,1);
4632 VFN_DEBUG_MESSAGE(
"PowderPattern::FitScaleFactorForIntegratedR():1",2);
4633 const long numInterval=mIntegratedPatternMin.numElements();
4634 CrystVector_REAL *integratedCalc=
new CrystVector_REAL[nbScale];
4635 for(
int i=0;i<nbScale;i++)
4637 integratedCalc[i].resize(numInterval);
4641 const REAL *p1=mPowderPatternComponentRegistry.GetObj(mScalableComponentIndex(i))
4642 .mPowderPatternCalc.data();
4644 REAL *p2=integratedCalc[i].data();
4645 for(
int j=0;j<numInterval;j++)
4647 const long max=mIntegratedPatternMax(j);
4648 p1=mPowderPatternComponentRegistry.GetObj(mScalableComponentIndex(i))
4649 .mPowderPatternCalc.data()+mIntegratedPatternMin(j);
4651 for(
int k=mIntegratedPatternMin(j);k<=max;k++) *p2 += *p1++;
4658 VFN_DEBUG_MESSAGE(
"PowderPattern::FitScaleFactorForIntegratedR():2",2);
4659 CrystVector_REAL backdIntegrated(numInterval);
4663 REAL *p2=backdIntegrated.data();
4664 for(
int j=0;j<numInterval;j++)
4666 const long max=mIntegratedPatternMax(j);
4669 for(
int k=mIntegratedPatternMin(j);k<=max;k++) *p2 += *p1++;
4681 VFN_DEBUG_MESSAGE(
"PowderPattern::FitScaleFactorForIntegratedR():3",2);
4682 for(
int i=0;i<nbScale;i++)
4684 for(
int j=i;j<nbScale;j++)
4686 const REAL *p1=integratedCalc[i].data();
4687 const REAL *p2=integratedCalc[j].data();
4689 for(
long k=0;k<numInterval;k++)
4694 mFitScaleFactorM(i,j)=m;
4695 mFitScaleFactorM(j,i)=m;
4698 VFN_DEBUG_MESSAGE(
"PowderPattern::FitScaleFactorForIntegratedR():4",2);
4699 for(
int i=0;i<nbScale;i++)
4701 const REAL *p1=mIntegratedObs.data();
4702 const REAL *p2=integratedCalc[i].data();
4705 for(
long k=0;k<numInterval;k++)
4712 const REAL *p3=backdIntegrated.data();
4713 for(
long k=0;k<numInterval;k++)
4718 b += (*p1++ - *p3++) * *p2++;
4721 mFitScaleFactorB(i,0) =b;
4723 VFN_DEBUG_MESSAGE(
"PowderPattern::FitScaleFactorForIntegratedR():5",2);
4725 if(1==nbScale) mFitScaleFactorX=mFitScaleFactorB(0)/mFitScaleFactorM(0);
4727 mFitScaleFactorX=product(InvertMatrix(mFitScaleFactorM),mFitScaleFactorB);
4728 VFN_DEBUG_MESSAGE(
"B, M, X"<<endl<<mFitScaleFactorB<<endl<<mFitScaleFactorM<<endl<<mFitScaleFactorX,3)
4729 for(
int i=0;i<nbScale;i++)
4731 const REAL * p1=mPowderPatternComponentRegistry.GetObj(mScalableComponentIndex(i))
4732 .mPowderPatternCalc.data();
4734 const REAL s = mFitScaleFactorX(i)
4735 -mScaleFactor(mScalableComponentIndex(i));
4738 (*fpObjCrystInformUser)(
"Warning:FitScaleFactorForIntegratedR: working around NaN scale factor...");
4741 for(
unsigned long j=0;j<mNbPointUsed;j++) *p0++ += s * *p1++;
4742 VFN_DEBUG_MESSAGE(
"-> Old:"<<mScaleFactor(mScalableComponentIndex(i)) <<
" New:"<<mFitScaleFactorX(i),3);
4743 mScaleFactor(mScalableComponentIndex(i)) = mFitScaleFactorX(i);
4747 delete[] integratedCalc;
4748 VFN_DEBUG_EXIT(
"PowderPattern::FitScaleFactorForIntegratedR():End",3);
4758 TAU_PROFILE(
"PowderPattern::FitScaleFactorForRw()",
"void ()",TAU_DEFAULT);
4759 VFN_DEBUG_ENTRY(
"PowderPattern::FitScaleFactorForRw()",3);
4763 mScalableComponentIndex.resize(mPowderPatternComponentRegistry.GetNb());
4765 for(
int i=0;i<mPowderPatternComponentRegistry.GetNb();i++)
4767 if(mPowderPatternComponentRegistry.GetObj(i).IsScalable())
4768 mScalableComponentIndex(nbScale++)=i;
4770 VFN_DEBUG_MESSAGE(
"-> Number of Scale Factors:"<<nbScale<<
":Index:"<<endl<<mScalableComponentIndex,2);
4773 VFN_DEBUG_EXIT(
"PowderPattern::FitScaleFactorForRw(): No scalable component!",3);
4776 mScalableComponentIndex.resizeAndPreserve(nbScale);
4778 mFitScaleFactorM.resize(nbScale,nbScale);
4779 mFitScaleFactorB.resize(nbScale,1);
4780 mFitScaleFactorX.resize(nbScale,1);
4785 for(
int i=0;i<nbScale;i++)
4787 for(
int j=i;j<nbScale;j++)
4791 const REAL *p1=mPowderPatternComponentRegistry.GetObj(mScalableComponentIndex(i))
4792 .mPowderPatternCalc.data();
4793 const REAL *p2=mPowderPatternComponentRegistry.GetObj(mScalableComponentIndex(j))
4794 .mPowderPatternCalc.data();
4797 for(
unsigned long k=0;k<mNbPointUsed;k++) m += *p1++ * *p2++ * *p3++;
4798 mFitScaleFactorM(i,j)=m;
4799 mFitScaleFactorM(j,i)=m;
4802 for(
int i=0;i<nbScale;i++)
4805 const REAL *p2=mPowderPatternComponentRegistry.GetObj(mScalableComponentIndex(i))
4806 .mPowderPatternCalc.data();
4810 for(
unsigned long k=0;k<mNbPointUsed;k++) b += *p1++ * *p2++ * *p3++;
4814 for(
unsigned long k=0;k<mNbPointUsed;k++)
4815 b += (*p1++ - *p4++) * *p2++ * *p3++;
4817 mFitScaleFactorB(i,0) =b;
4822 unsigned long min,max;
4823 for(
int i=0;i<nbScale;i++)
4825 for(
int j=i;j<nbScale;j++)
4829 const REAL *p1=mPowderPatternComponentRegistry.GetObj(mScalableComponentIndex(i))
4830 .mPowderPatternCalc.data();
4831 const REAL *p2=mPowderPatternComponentRegistry.GetObj(mScalableComponentIndex(j))
4832 .mPowderPatternCalc.data();
4836 for(
int k=0;k<nbExclude;k++)
4840 if(min>mNbPointUsed)
break;
4841 if(max>mNbPointUsed)max=mNbPointUsed;
4843 for(;l<min;l++) m += *p1++ * *p2++ * *p3++;
4849 for(;l<mNbPointUsed;l++) m += *p1++ * *p2++ * *p3++;
4850 mFitScaleFactorM(i,j)=m;
4851 mFitScaleFactorM(j,i)=m;
4854 for(
int i=0;i<nbScale;i++)
4857 const REAL *p2=mPowderPatternComponentRegistry
4858 .GetObj(mScalableComponentIndex(i))
4859 .mPowderPatternCalc.data();
4865 for(
int k=0;k<nbExclude;k++)
4869 if(min>mNbPointUsed)
break;
4870 if(max>mNbPointUsed)max=mNbPointUsed;
4872 for(;l<min;l++) b += *p1++ * *p2++ * *p3++;
4878 for(;l<mNbPointUsed;l++) b += *p1++ * *p2++ * *p3++;
4883 for(
int k=0;k<nbExclude;k++)
4887 if(min>mNbPointUsed)
break;
4888 if(max>mNbPointUsed)max=mNbPointUsed;
4890 for(;l<min;l++) b += (*p1++ - *p4++) * *p2++ * *p3++;
4896 for(;l<mNbPointUsed;l++) b += (*p1++ - *p4++) * *p2++ * *p3++;
4898 mFitScaleFactorB(i,0) =b;
4901 if(1==nbScale) mFitScaleFactorX=mFitScaleFactorB(0)/mFitScaleFactorM(0);
4903 mFitScaleFactorX=product(InvertMatrix(mFitScaleFactorM),mFitScaleFactorB);
4904 VFN_DEBUG_MESSAGE(
"B, M, X"<<endl<<mFitScaleFactorB<<endl<<mFitScaleFactorM<<endl<<mFitScaleFactorX,2)
4905 for(
int i=0;i<nbScale;i++)
4907 const REAL * p1=mPowderPatternComponentRegistry.GetObj(mScalableComponentIndex(i))
4908 .mPowderPatternCalc.data();
4910 const REAL s = mFitScaleFactorX(i)
4911 -mScaleFactor(mScalableComponentIndex(i));
4914 (*fpObjCrystInformUser)(
"Warning:FitScaleFactorForRw working around NaN scale factor...");
4917 for(
unsigned long j=0;j<mNbPointUsed;j++) *p0++ += s * *p1++;
4918 VFN_DEBUG_MESSAGE(
"-> Old:"<<mScaleFactor(mScalableComponentIndex(i)) <<
" Change:"<<mFitScaleFactorX(i),3);
4919 mScaleFactor(mScalableComponentIndex(i)) = mFitScaleFactorX(i);
4923 VFN_DEBUG_EXIT(
"PowderPattern::FitScaleFactorForRw():End",3);
4926 void PowderPattern::FitScaleFactorForIntegratedRw()
const
4935 VFN_DEBUG_ENTRY(
"PowderPattern::FitScaleFactorForIntegratedRw()",3);
4936 TAU_PROFILE(
"PowderPattern::FitScaleFactorForIntegratedRw()",
"void ()",TAU_DEFAULT);
4938 mScalableComponentIndex.resize(mPowderPatternComponentRegistry.GetNb());
4939 CrystVector_REAL vScalableVarianceIndex(mPowderPatternComponentRegistry.GetNb());
4942 for(
int i=0;i<mPowderPatternComponentRegistry.GetNb();i++)
4944 if(mPowderPatternComponentRegistry.GetObj(i).IsScalable())
4945 mScalableComponentIndex(nbScale++)=i;
4946 if(mPowderPatternComponentRegistry.GetObj(i).HasPowderPatternCalcVariance())
4948 if(0!=mPowderPatternComponentRegistry.GetObj(i)
4949 .GetPowderPatternIntegratedCalcVariance().first->numElements())
4950 vScalableVarianceIndex(nbVarCalc++)=i;
4953 VFN_DEBUG_MESSAGE(
"-> Number of Scale Factors:"<<nbScale<<
"("<<nbVarCalc
4954 <<
"with variance). Index:"<<endl<<mScalableComponentIndex,2);
4957 VFN_DEBUG_EXIT(
"PowderPattern::FitScaleFactorForIntegratedRw(): No scalable component!",3);
4961 unsigned int ctagain=0;
4962 RECALC_SCALE_FACTOR_VARIANCE_FitScaleFactorForIntegratedRw:
4966 double a=0.,b=0.,c=0.;
4969 const REAL * RESTRICT p1=mIntegratedObs.data();
4970 const REAL * RESTRICT p2=mPowderPatternComponentRegistry.GetObj(mScalableComponentIndex(0))
4971 .GetPowderPatternIntegratedCalc().first->data();
4972 const REAL * RESTRICT p1v=mIntegratedVarianceObs.data();
4973 const REAL * RESTRICT p2v=mPowderPatternComponentRegistry.GetObj(mScalableComponentIndex(0))
4974 .GetPowderPatternIntegratedCalcVariance().first->data();
4977 for(
unsigned long k=mNbIntegrationUsed;k>0;k--)
4979 a += *p2v * *p1 * *p2;
4980 b += *p2 * *p2 * *p1v - *p1 * *p1 * *p2v;
4981 c -= *p1 * *p2 * *p1v;
4982 p1++;p2++;p1v++;p2v++;
4988 for(
unsigned long k=mNbIntegrationUsed;k>0;k--)
4990 a += *p2v * (*p1 - *p3) * *p2;
4991 b += *p2 * *p2 * *p1v - (*p1 - *p3) * (*p1 - *p3) * *p2v;
4992 c -= (*p1 - *p3) * *p2 * *p1v;
4993 p1++;p2++;p1v++;p2v++;p3++;
4997 newscale=(REAL)((-b+sqrt(b*b-4*a*c))/(2*a));
5000 const REAL s = newscale-mScaleFactor(mScalableComponentIndex(0));
5001 const REAL s2 = newscale*newscale
5002 -mScaleFactor(mScalableComponentIndex(0))
5003 *mScaleFactor(mScalableComponentIndex(0));
5005 const REAL * RESTRICT p1=mPowderPatternComponentRegistry.GetObj(mScalableComponentIndex(0))
5006 .GetPowderPatternIntegratedCalc().first->data();
5008 const REAL * RESTRICT p1v=mPowderPatternComponentRegistry.GetObj(mScalableComponentIndex(0))
5009 .GetPowderPatternIntegratedCalcVariance().first->data();
5010 REAL * RESTRICT p0w = mIntegratedWeight.data();
5011 for(
unsigned long j=mNbIntegrationUsed;j>0;j--)
5014 *p0v += s2 * *p1v++;
5015 if(*p0v <=0) {*p0w++ =0;p0v++;}
else *p0w++ = 1. / *p0v++;
5018 mScaleFactor(mScalableComponentIndex(0)) = newscale;
5025 mScalableComponentIndex.resizeAndPreserve(nbScale);
5027 mFitScaleFactorM.resize(nbScale,nbScale);
5028 mFitScaleFactorB.resize(nbScale,1);
5029 mFitScaleFactorX.resize(nbScale,1);
5031 VFN_DEBUG_MESSAGE(
"PowderPattern::FitScaleFactorForIntegratedRw():1",2);
5032 vector<const CrystVector_REAL*> integratedCalc;
5033 for(
int i=0;i<nbScale;i++)
5035 integratedCalc.push_back(mPowderPatternComponentRegistry.
5036 GetObj(mScalableComponentIndex(i))
5037 .GetPowderPatternIntegratedCalc().first);
5039 VFN_DEBUG_MESSAGE(
"PowderPattern::FitScaleFactorForIntegratedRw():3",2);
5040 for(
int i=0;i<nbScale;i++)
5042 for(
int j=i;j<nbScale;j++)
5044 const REAL * RESTRICT p1=integratedCalc[i]->data();
5045 const REAL * RESTRICT p2=integratedCalc[j]->data();
5046 const REAL * RESTRICT p3;
5047 if(mIntegratedWeight.numElements()==0)
5049 p3=mIntegratedWeightObs.data();
5050 if(ctagain>5) VFN_DEBUG_MESSAGE(
"ctagain="<<ctagain<<
", using mIntegratedWeightObs",5);
5054 p3=mIntegratedWeight.data();
5055 if(ctagain>5) VFN_DEBUG_MESSAGE(
"ctagain="<<ctagain<<
", using mIntegratedWeight",5);
5059 for(
unsigned long k=mNbIntegrationUsed;k>0;k--)
5060 {m += *p1 * *p1 * *p3++; p1++;}
5062 for(
unsigned long k=mNbIntegrationUsed;k>0;k--)
5063 m += *p1++ * *p2++ * *p3++;
5064 mFitScaleFactorM(i,j)=m;
5065 mFitScaleFactorM(j,i)=m;
5068 VFN_DEBUG_MESSAGE(
"PowderPattern::FitScaleFactorForIntegratedRw():4",2);
5069 for(
int i=0;i<nbScale;i++)
5071 const REAL * RESTRICT p1=mIntegratedObs.data();
5072 const REAL * RESTRICT p2=integratedCalc[i]->data();
5073 const REAL * RESTRICT p4;
5074 if(mIntegratedWeight.numElements()==0) p4=mIntegratedWeightObs.data();
5075 else p4=mIntegratedWeight.data();
5079 for(
unsigned long k=mNbIntegrationUsed;k>0;k--)
5081 b += *p1++ * *p2++ * *p4++;
5088 for(
unsigned long k=mNbIntegrationUsed;k>0;k--)
5093 b += (*p1++ - *p3++) * *p2++ * *p4++;
5096 mFitScaleFactorB(i,0) =b;
5098 VFN_DEBUG_MESSAGE(
"PowderPattern::FitScaleFactorForIntegratedRw():5",2);
5100 if(1==nbScale) mFitScaleFactorX=mFitScaleFactorB(0)/mFitScaleFactorM(0);
5102 mFitScaleFactorX=product(InvertMatrix(mFitScaleFactorM),mFitScaleFactorB);
5103 VFN_DEBUG_MESSAGE(
"B, M, X"<<endl<<mFitScaleFactorB<<endl<<mFitScaleFactorM<<endl<<mFitScaleFactorX,3)
5108 for(
int i=0;i<nbScale;i++)
5110 if(mPowderPatternComponentRegistry.GetObj(mScalableComponentIndex(i)).HasPowderPatternCalcVariance())
5112 if(0!=mPowderPatternComponentRegistry.GetObj(mScalableComponentIndex(i))
5113 .GetPowderPatternIntegratedCalcVariance().first->numElements())
5115 const REAL * RESTRICT p1=mPowderPatternComponentRegistry.GetObj(mScalableComponentIndex(i))
5116 .GetPowderPatternIntegratedCalcVariance().first->data();
5120 const REAL s2 = mFitScaleFactorX(i)*mFitScaleFactorX(i)
5121 -mScaleFactor(mScalableComponentIndex(i))
5122 *mScaleFactor(mScalableComponentIndex(i));
5123 for(
unsigned long j=mNbIntegrationUsed;j>0;j--) *p0++ += s2 * *p1++;
5129 REAL * RESTRICT p0 = mIntegratedWeight.data();
5131 for(
unsigned long j=mNbIntegrationUsed;j>0;j--)
5132 if(*p1 <=0) {*p0++ =0;p1++;}
5133 else *p0++ = 1. / *p1++;
5137 for(
int i=0;i<nbScale;i++)
5139 const REAL * RESTRICT p1=mPowderPatternComponentRegistry.GetObj(mScalableComponentIndex(i))
5140 .GetPowderPatternIntegratedCalc().first->data();
5142 const REAL s = mFitScaleFactorX(i)
5143 -mScaleFactor(mScalableComponentIndex(i));
5146 (*fpObjCrystInformUser)(
"Warning:FitScaleFactorForIntegratedRw: working around NaN scale factor...");
5151 if(abs(s/mFitScaleFactorX(i))>0.001)
5157 if((!again)&&(ctagain>0))
5159 VFN_DEBUG_MESSAGE(
"log(scale) :"<<log(mScaleFactor(mScalableComponentIndex(i)))
5160 <<
"->"<<log(mFitScaleFactorX(i))<<
" Again="<<ctagain,5);
5163 for(
unsigned long j=mNbIntegrationUsed;j>0;j--) *p0++ += s * *p1++;
5166 VFN_DEBUG_MESSAGE(
"->log(scale) Old :"<<log(mScaleFactor(mScalableComponentIndex(i))) <<
" New:"<<log(mFitScaleFactorX(i)),10);
5168 mScaleFactor(mScalableComponentIndex(i)) = mFitScaleFactorX(i);
5174 if(again && (ctagain<20))
5179 VFN_DEBUG_MESSAGE(
"PowderPattern::FitScaleFactorForIntegratedRw(), scaling again #"<<ctagain,10)
5181 goto RECALC_SCALE_FACTOR_VARIANCE_FitScaleFactorForIntegratedRw;
5183 VFN_DEBUG_EXIT(
"PowderPattern::FitScaleFactorForIntegratedRw():End",3);
5188 VFN_DEBUG_MESSAGE(
"PowderPattern::SetSigmaToSqrtIobs()",5);
5194 VFN_DEBUG_MESSAGE(
"PowderPattern::SetWeightToInvSigmaSq()",5);
5208 VFN_DEBUG_MESSAGE(
"PowderPattern::SetWeightToSinTheta()",5);
5214 const REAL minRelatIobs)
5216 VFN_DEBUG_MESSAGE(
"PowderPattern::SetWeightPolynomial()",5);
5229 const bool enableRestraints)
5232 if(0 == mOptProfileIntegration.GetChoice()) this->FitScaleFactorForIntegratedRw();
5254 VFN_DEBUG_MESSAGE(
"PowderPattern::AddExcludedRegion()",5)
5270 CrystVector_long subs;
5272 CrystVector_REAL tmp1,tmp2;
5281 VFN_DEBUG_MESSAGE(
"PowderPattern::Add2ThetaExcludedRegion():End",5)
5287 if(mOptProfileIntegration.GetChoice()==0) tmp+=mIntegratedChi2LikeNorm;
5288 else tmp+=mChi2LikeNorm;
5294 const CrystVector_REAL&
5297 TAU_PROFILE(
"PowderPattern::GetLSQCalc()",
"void ()",TAU_DEFAULT);
5316 const CrystVector_REAL&
5319 TAU_PROFILE(
"PowderPattern::GetLSQObs()",
"void ()",TAU_DEFAULT);
5338 const CrystVector_REAL&
5341 TAU_PROFILE(
"PowderPattern::GetLSQWeight()",
"void ()",TAU_DEFAULT);
5349 if(mIntegratedWeight.numElements()==0)
5366 TAU_PROFILE(
"PowderPattern::GetLSQ_FullDeriv()",
"void ()",TAU_DEFAULT);
5370 this->CalcPowderPatternIntegrated_FullDeriv(vPar);
5372 std::map<RefinablePar*, CrystVector_REAL> fullderiv_old;
5373 std::vector<const CrystVector_REAL*> v;
5376 cout<<
"PowderPattern::GetLSQ_FullDeriv(integrated):parameters:"<<endl;
5377 for(std::set<RefinablePar*>::iterator par=vPar.begin();par!=vPar.end();++par)
5379 v.push_back(&(mPowderPatternIntegrated_FullDeriv[*par]));
5380 fullderiv_old[*par]=this->
GetLSQDeriv(idx,*(*par));
5381 v.push_back(&(fullderiv_old[*par]));
5382 cout<<(*par)->GetName()<<
":"<<mPowderPatternIntegrated_FullDeriv[*par].size()<<
","<<fullderiv_old[*par].size()<<endl;
5385 cout<<
"PowderPattern::GetLSQ_FullDeriv(integrated):"<<endl<<FormatVertVector<REAL>(v,12,1,20)<<endl;
5388 return mPowderPatternIntegrated_FullDeriv;
5390 mPowderPattern_FullDeriv=this->GetPowderPattern_FullDeriv(vPar);
5391 return mPowderPattern_FullDeriv;
5396 VFN_DEBUG_MESSAGE(
"PowderPattern::Prepare()",5);
5397 for(
int i=0;i<mPowderPatternComponentRegistry.GetNb();i++)
5400 mPowderPatternComponentRegistry.GetObj(i).Prepare();
5404 CrystVector_uint & groupIndex,
5405 unsigned int &first)
const
5408 unsigned int scaleIndex=0;
5409 unsigned int thetaIndex=0;
5410 VFN_DEBUG_MESSAGE(
"PowderPattern::GetGeneGroup()",4)
5412 for(
long j=0;j<this->
GetNbPar();j++)
5417 if(scaleIndex==0) scaleIndex=first++;
5418 groupIndex(i)=scaleIndex;
5422 if(thetaIndex==0) thetaIndex=first++;
5423 groupIndex(i)=thetaIndex;
5431 for(
int i=0;i<mPowderPatternComponentRegistry.GetNb();i++)
5442 return mIntegratedPatternMin;
5448 return mIntegratedPatternMax;
5453 return mClockIntegratedFactorsPrep;
5459 if(this->
GetRadiation().GetWavelengthType()==WAVELENGTH_TOF)
5461 if(stol>0) x = 1.0/(2*stol);
else return 0;
5462 x = mDIFC*x+mDIFA*x*x;
5463 VFN_DEBUG_MESSAGE(
"PowderPattern::STOL2X("<<stol<<
","<<1.0/(2*stol)<<
")="<<x,2)
5468 if(abs(x)<1.0) x=2*asin(x);
else x=2*M_PI;
5476 if(this->
GetRadiation().GetWavelengthType()==WAVELENGTH_TOF)
5478 if(abs(mDIFA)>abs(mDIFC*1e-6))
5480 const REAL delta=mDIFC*mDIFC+4.0*mDIFA*x;
5481 stol = (-mDIFC+sqrt(delta))/(2.0*mDIFA);
5482 stol = 1/(2.0*stol);
5484 else stol=mDIFC/(2.0*x);
5488 VFN_DEBUG_MESSAGE(
"PowderPattern::X2STOL("<<x<<
")="<<stol,2)
5502 if(this->
GetRadiation().GetWavelengthType()!=WAVELENGTH_TOF)
5505 for(finish=0;finish<nb;++finish)
5514 for(start=nb-1;start>=0;--start)
5521 unsigned int width_golay=10;
5523 CrystVector_REAL obs;
5524 const int nbwidth=9;
5525 CrystVector_long width(nbwidth);
5535 if(max>=obs.numElements()) max = obs.numElements();
5536 for(
long j=min;j<max;j++) obs(j) = 0;
5538 const long nb=obs.numElements();
5539 for(
int j=0;j<nbwidth;j++)
5541 const long imax=obs.imax(nb/10,nb-1);
5542 const REAL iobs_max=obs(imax);
5543 REAL thres=iobs_max;
5545 for(i=imax-100;i<(imax+100);++i)
5547 if(i<0){i=0;
continue;}
5549 if(obs(i)<thres) thres=obs(i);
5551 thres=(iobs_max+thres)/2;
5553 while(obs(i)>=thres)
5561 while(obs(i)>=thres)
5568 cout<<endl<<
" => "<<width(j)<<endl;
5569 for(i=imax-width(j)*2;i<(imax+width(j)*2);++i)
5576 cout<<
"Width of "<<nbwidth<<
" strongest peaks:"<<endl<<width;
5577 width_golay=width(SortSubs(width)(nbwidth/2));
5578 cout<<
"median width:"<<width_golay<<endl;
5579 if(width_golay<=4)width_golay=4;
5580 if(width_golay>=16)width_golay=16;
5584 CrystVector_REAL obsd2;
5593 if(max>=obsd2.numElements()) max = obsd2.numElements();
5594 for(
long j=min;j<max;j++) obsd2(j) = 0;
5596 const float norm=-obsd2.min();
5603 CrystVector_REAL tmp;
5605 tmp.resizeAndPreserve(tmp.numElements()/4);
5607 QuickSortSubs(tmp,sub,tmp.numElements()-1,0);
5608 min_iobs=5*(tmp(tmp.numElements()/2)-tmp(tmp.numElements()/4));
5615 const long imax=obsd2.imax(start,finish);
5616 REAL iobs=obsd2(imax);
5618 REAL xmax=
mX(imax)*iobs;
5621 REAL lastiobs=obsd2(i);
5622 const REAL iobs_max=lastiobs;
5627 if(obsd2(--i)>=lastiobs)
break;
5631 xmax+=
mX(i)*lastiobs;nbav++;
5632 if(lastiobs<=0)
break;
5635 float dleft=
mX(i+1);
5640 if(i>=(nb-2))
break;
5641 if(obsd2(++i)>=lastiobs)
break;
5645 xmax+=
mX(i)*lastiobs;nbav++;
5646 if(lastiobs<=0)
break;
5649 float dright=
mX(i-1);
5652 REAL dmax=this->
X2STOL(xmax)*2;
5653 dright=this->
X2STOL(dright)*2;
5654 dleft =this->
X2STOL(dleft)*2;
5659 min_iobs=iobs/nbav*maxratio;
5665 if(nbav_min<3)nbav_min=3;
5671 cout<<
"Peak #"<<pl.
GetPeakList().size()<<
"imax="<<imax<<
", x="<<xmax<<
",d="<<1/dmax<<
", d2iobs_max="<<iobs_max
5672 <<
", d2Iobs="<<iobs<<
", nbav="<<nbav<<
", min_iobs="<<min_iobs<<
",sigma="<<sigma<<endl;
5673 if( ((nbav>=nbav_min) &&(iobs_max>min_iobs)&&((iobs/nbav)>min_iobs))
5674 ||((nbav>=nbav_min) &&(iobs_max>min_iobs)&&((iobs/nbav)>min_iobs*.2)&&((iobs/nbav)>3*sigma))
5675 ||((nbav>=nbav_min/2)&&(iobs_max>min_iobs)&&((iobs/nbav)>min_iobs*2 )&&((iobs/nbav)>6*sigma)))
5677 pl.
AddPeak(dmax,iobs,abs(dright-dleft)*.25);
5678 if((pl.
GetPeakList().size()==1)&&(maxratio<0)&&(min_iobs<0.005*iobs/nbav)) min_iobs=0.005*iobs/nbav;
5692 exportAtom(
string n,REAL X, REAL Y, REAL Z,REAL b,REAL o,
const ScatteringPower *pow):
5693 name(n),x(X),y(Y),z(Z),biso(b),occ(o),occMult(1),mpScattPow(pow){}
5695 REAL x,y,z,biso,occ;
5698 const ScatteringPower *mpScattPow;
5703 exportBond(
const string &a1,
const string &a2, REAL d, REAL s):
5704 at1(a1),at2(a2),dist(d),sigma(s){}
5711 exportAngle(
const string &a1,
const string &a2,
const string &a3, REAL a, REAL s):
5712 at1(a1),at2(a2),at3(a3),ang(a),sigma(s){}
5721 vector<const PowderPatternDiffraction*> vDiff;
5729 if((pBackground==0)||vDiff.size()==0)
return;
5732 ofstream dat((prefix+
".dat").c_str());
5734 <<
"INTER 1.0 1.0 0"<<endl<<endl<<endl<<endl<<endl;
5736 CrystVector_REAL ttheta;
5738 if(this->
GetRadiation().GetWavelengthType()!=WAVELENGTH_TOF) ttheta *= RAD2DEG;
5743 ofstream pcr((prefix+
".pcr").c_str());
5747 pcr<<
"Fox/ObjCryst exported file:"<<this->
GetName()<<endl;
5749 pcr<<
"NPATT 1"<<endl;
5751 pcr<<
"W_PAT 1.0"<<endl;
5753 pcr<<
"! Nph Dum Ias Nre Cry Opt Aut"<<endl;
5754 pcr<<
" "<<vDiff.size()
5755 <<
" 0 0 0 0 1 1 "<<endl;
5762 pcr<<
"! Job Npr Nba Nex Nsc Nor Iwg Ilo Res Ste Uni Cor Anm"<<endl
5764 <<
" 5 "<<pBackground->GetInterpPoints().first->numElements()
5765 <<
" 0 0 1 0 0 0 1 0 0 0"<<endl;
5770 std::string::size_type idx =prefix.rfind(
"/");
5771 std::string::size_type idx2=prefix.rfind(
"\\");
5772 std::string::size_type idx3=prefix.rfind(
":");
5773 if(((
long)idx2!=(
long)string::npos)&&((
long)idx2>(
long)idx))idx=idx2;
5774 if(((
long)idx3!=(
long)string::npos)&&((
long)idx3>(
long)idx))idx=idx3;
5775 if(idx==string::npos)
5778 shortName=prefix.substr(idx+1);
5780 pcr<<
"! File names of data files"<<endl;
5781 pcr<<shortName<<
".dat"<<endl;
5783 pcr<<
"! Mat Pcr Syo Rpa Sym Sho"<<endl
5784 <<
" 1 1 0 -1 0 0 "<<endl;
5786 pcr<<
"! Ipr Ppl Ioc Ls1 Ls2 Ls3 Prf Ins Hkl Fou Ana"<<endl
5787 <<
" 0 0 0 0 0 0 1 10 0 0 1 "<<endl;
5791 pcr<<
"!lambda1 lambda2 Ratio Bkpos Wdt Cthm muR AsyLim Rpolarz -> Patt #1"<<endl
5794 <<
" 0 0 0 0.95"<<endl;
5796 pcr<<
"!NCY Eps R_at R_an R_pr R_gl"<<endl
5797 <<
" 5 0.2 1.0 1.0 1.0 1.0"<<endl;
5799 pcr<<
"! Thmin Step Thmax PSD Sent0 -> Patt #1"<<endl
5800 <<
" 0 0 0 0 0"<<endl;
5802 pcr<<
"!2Theta Background for Pattern #1"<<endl;
5803 for(
long i=0;i<pBackground->GetInterpPoints().first->numElements();i++)
5804 pcr<<(*(pBackground->GetInterpPoints().first))(i)*RAD2DEG<<
" "
5805 <<(*(pBackground->GetInterpPoints().second))(i)<<
" 0.0"<<endl;
5807 pcr<<
"!"<<endl<<
"!"<<endl<<
"1 !Number of refined parameters"<<endl;
5809 pcr<<
"! Zero Code Sycos Code Sysin Code Lambda Code More -> Patt #1"<<endl;
5810 pcr<<
" "<<
mXZero*RAD2DEG <<
" 0.0 "
5811 <<m2ThetaDisplacement*RAD2DEG <<
" 0.0 "
5812 <<m2ThetaTransparency*RAD2DEG <<
" 0.0 "
5813 <<
"0.000 0.0 0"<<endl;
5815 for(
unsigned int i=0;i<vDiff.size();++i)
5817 pcr<<
"!-------------------------------------------------------------------------------"<<endl
5818 <<
"! Data for PHASE number: "<<i<<
" ==> Current R_Bragg for Pattern# 1: 0.00 "<<endl
5819 <<
"!-------------------------------------------------------------------------------"<<endl;
5821 pcr<<vDiff[i]->GetCrystal().GetName()<<endl;
5824 map<int,exportAtom> vExportAtom;
5825 list<exportBond> vExportBond;
5826 list<exportAngle> vExportAngle;
5828 CrystMatrix_REAL minDistTable;
5829 minDistTable=vDiff[i]->GetCrystal().GetMinDistanceTable(-1.);
5834 for(
int s=0;s<vDiff[i]->GetCrystal().GetScattererRegistry().GetNb();s++)
5841 if(vDiff[i]->GetCrystal().GetScatt(s).GetClassName()==
"Molecule")
5842 pMol=
dynamic_cast<const Molecule*
>(&(vDiff[i]->GetCrystal().GetScatt(s)));
5843 map<const MolAtom*,string> vMolAtomName;
5847 if(0==list(j).mpScattPow)
continue;
5848 bool redundant=
false;
5849 for(
unsigned long l=0;l<k;++l)
5850 if(abs(minDistTable(l,k))<0.5)
5852 map<int,exportAtom>::iterator pos=vExportAtom.find(l);
5853 if(pos!=vExportAtom.end()) pos->second.occMult+=1;
5860 name<<list(j).mpScattPow->GetName()<<k+1;
5861 vExportAtom.insert(make_pair(k,exportAtom(name.str(),
5862 list(j).
mX,list(j).mY,list(j).mZ,
5863 list(j).mpScattPow->GetBiso(),
5864 list(j).mOccupancy*list0(k).mDynPopCorr,
5865 list(j).mpScattPow)));
5866 if(pMol!=0) vMolAtomName.insert(make_pair(pMol->GetAtomList()[j],name.str()));
5872 for(vector<MolBond*>::const_iterator pos=pMol->GetBondList().begin();
5873 pos!=pMol->GetBondList().end();++pos)
5875 map<const MolAtom*,string>::const_iterator p1,p2;
5876 p1=vMolAtomName.find(&((*pos)->GetAtom1()));
5877 p2=vMolAtomName.find(&((*pos)->GetAtom2()));
5878 if( (p1!=vMolAtomName.end()) && (p2!=vMolAtomName.end()))
5879 vExportBond.push_back(exportBond(p1->second, p2->second,
5880 (*pos)->GetLength0(),(*pos)->GetLengthSigma()));
5883 for(vector<MolBondAngle*>::const_iterator pos=pMol->GetBondAngleList().begin();
5884 pos!=pMol->GetBondAngleList().end();++pos)
5886 map<const MolAtom*,string>::const_iterator p1,p2,p3;
5887 p1=vMolAtomName.find(&((*pos)->GetAtom1()));
5888 p2=vMolAtomName.find(&((*pos)->GetAtom2()));
5889 p3=vMolAtomName.find(&((*pos)->GetAtom3()));
5890 if( (p1!=vMolAtomName.end()) && (p2!=vMolAtomName.end()) && (p3!=vMolAtomName.end()))
5891 vExportAngle.push_back(exportAngle(p1->second, p2->second,p3->second,
5892 (*pos)->GetAngle0(),(*pos)->GetAngleSigma()));
5903 pcr<<
"!Nat Dis Ang Jbt Isy Str Furth ATZ Nvk More"<<endl
5904 << vExportAtom.size() <<
" "<<vExportBond.size()<<
" "<<vExportAngle.size()
5905 <<
" 0 0 0 0 1.0 0 1"<<endl;
5906 pcr<<
"!Jvi Jdi Hel Sol Mom Ter N_Domains"<<endl
5907 <<
" 0 3 0 0 0 0 0"<<endl;
5909 pcr<<
"!Contributions (0/1) of this phase to the patterns"<<endl
5912 pcr<<
"!Irf Npr Jtyp Nsp_Ref Ph_Shift for Pattern#"<<i<<endl
5913 <<
" 0 0 0 0 0"<<endl;
5914 pcr<<
"! Pr1 Pr2 Pr3 Brind. Rmua Rmub Rmuc for Pattern#"<<i<<endl
5915 <<
" 1.0 1.0 1.0 1.0 0.0 0.0 0.0"<<endl;
5918 pcr<<
"! Max_dst(dist) (angles) Bond-Valence Calc."<<endl
5919 <<
" 2.7000 1.5000 0"<<endl;
5922 pcr<<vDiff[i]->GetCrystal().GetSpaceGroup().GetCCTbxSpg().match_tabulated_settings().hermann_mauguin()
5923 <<
" <- Space Group Symbol"<<endl;
5925 pcr<<
"!Atom Typ X Y Z Biso Occ In Fin N_t Spc / Codes"<<endl;
5926 for(map<int,exportAtom>::const_iterator pos=vExportAtom.begin();pos!=vExportAtom.end();++pos)
5928 pcr<<pos->second.name
5929 <<
" "<<pos->second.mpScattPow->GetSymbol()<<
" "
5930 <<pos->second.x<<
" "<<pos->second.y<<
" "<<pos->second.z<<
" "
5931 <<pos->second.biso<<
" "
5932 <<pos->second.occ*pos->second.occMult<<
" 0 0 0 0"<<endl
5933 <<
" 0 0 0 0 0"<<endl;
5936 REAL eta0=vDiff[0]->GetProfile().GetPar(
"Eta0").GetHumanValue();
5937 if(eta0<.01) eta0=.01;
5938 else if(eta0>.99) eta0=.99;
5939 pcr<<
"!Scale Shape1 Bov Str1 Str2 Str3 Strain-Model"<<endl
5941 <<
" 0.0 0.0 0.0 0.0 0"<<endl
5942 <<
" 1.0 0.0 0.0 0.0 0.0 0.0 0"<<endl;
5945 pcr<<
"! U V W X Y GauSiz LorSiz Size-Model"<<endl;
5948 pcr<<vDiff[i]->GetProfile().GetPar(
"U").GetHumanValue()<<
" ";
5949 pcr<<vDiff[i]->GetProfile().GetPar(
"V").GetHumanValue()<<
" ";
5950 pcr<<vDiff[i]->GetProfile().GetPar(
"W").GetHumanValue()<<
" ";
5951 pcr<<
" 0.0 0.0 0.0 0.0 "<<endl
5952 <<
" 0.0 0.0 0.0 0.0 0.0 0.0 0.0 "<<endl;
5954 pcr<<
"! a b c alpha beta gamma #Cell Info"<<endl
5955 <<vDiff[i]->GetCrystal().GetLatticePar(0)<<
" "
5956 <<vDiff[i]->GetCrystal().GetLatticePar(1)<<
" "
5957 <<vDiff[i]->GetCrystal().GetLatticePar(2)<<
" "
5958 <<vDiff[i]->GetCrystal().GetLatticePar(3)*RAD2DEG<<
" "
5959 <<vDiff[i]->GetCrystal().GetLatticePar(4)*RAD2DEG<<
" "
5960 <<vDiff[i]->GetCrystal().GetLatticePar(5)*RAD2DEG<<endl
5961 <<
" 0.0 0.0 0.0 0.0 0.0 0.0"<<endl;
5962 pcr<<
"! Pref1 Pref2 alpha0 beta0 alpha1 beta1 ?"<<endl
5963 <<
" 0.0 0.0 0.0 0.0 0.0 0.0"<<endl
5964 <<
" 0.0 0.0 0.0 0.0 0.0 0.0"<<endl;
5968 if(vExportBond.size()>0)
5970 pcr<<
"!Soft distance constraints"<<endl;
5971 for(list<exportBond>::const_iterator pos=vExportBond.begin();pos!=vExportBond.end();++pos)
5973 pcr<<pos->at1<<
" "<<pos->at2<<
" 1 0 0 0 "<<pos->dist<<
" "<<pos->sigma<<endl;
5976 if(vExportBond.size()>0)
5978 pcr<<
"!Soft angle constraints"<<endl;
5979 for(list<exportAngle>::const_iterator pos=vExportAngle.begin();pos!=vExportAngle.end();++pos)
5981 pcr<<pos->at1<<
" "<<pos->at2<<
" "<<pos->at3<<
" 1 1 0 0 0 0 0 0 "
5982 <<pos->ang*RAD2DEG<<
" "<<pos->sigma*RAD2DEG<<endl;
5998 TAU_PROFILE(
"PowderPattern::CalcPowderPattern()",
"void ()",TAU_DEFAULT);
5999 VFN_DEBUG_ENTRY(
"PowderPattern::CalcPowderPattern()",3);
6000 if(mPowderPatternComponentRegistry.GetNb()==0)
6011 for(
unsigned long j=0;j<
mNbPoint;j++)
6013 if(*p0 <=0) {*p1 =0.;}
6014 else *p1 = 1. / *p0;
6017 VFN_DEBUG_EXIT(
"PowderPattern::CalcPowderPattern():no components!",3);
6020 TAU_PROFILE_TIMER(timer1,
"PowderPattern::CalcPowderPattern1()Calc components",
"", TAU_FIELD);
6021 TAU_PROFILE_TIMER(timer2,
"PowderPattern::CalcPowderPattern2(Add spectrums-scaled)"\
6023 TAU_PROFILE_TIMER(timer3,
"PowderPattern::CalcPowderPattern2(Add spectrums-backgd1)"\
6025 TAU_PROFILE_TIMER(timer4,
"PowderPattern::CalcPowderPattern2(Add spectrums-backgd2)"\
6027 TAU_PROFILE_TIMER(timer5,
"PowderPattern::CalcPowderPattern3(Variance)",
"", TAU_FIELD);
6028 TAU_PROFILE_START(timer1);
6029 for(
int i=0;i<mPowderPatternComponentRegistry.GetNb();i++)
6030 mPowderPatternComponentRegistry.GetObj(i).CalcPowderPattern();
6031 TAU_PROFILE_STOP(timer1);
6032 VFN_DEBUG_MESSAGE(
"PowderPattern::CalcPowderPattern():Calculated components..",3);
6039 for(
int i=0;i<mPowderPatternComponentRegistry.GetNb();i++)
6041 mPowderPatternComponentRegistry.GetObj(i).GetClockPowderPatternCalc())
6049 VFN_DEBUG_EXIT(
"PowderPattern::CalcPowderPattern():no need to recalc",3);
6054 for(
int i=0;i<mPowderPatternComponentRegistry.GetNb();i++)
6056 VFN_DEBUG_MESSAGE(
"PowderPattern::CalcPowderPattern():Adding "<< mPowderPatternComponentRegistry.GetObj(i).GetName(),3);
6057 if(
true==mPowderPatternComponentRegistry.GetObj(i).IsScalable())
6059 TAU_PROFILE_START(timer2);
6062 const REAL * p1=mPowderPatternComponentRegistry.GetObj(i)
6063 .mPowderPatternCalc.data();
6065 const REAL s = mScaleFactor(i);
6066 for(
unsigned long j=0;j<mNbPointUsed;j++) *p0++ = s * *p1++;
6071 const REAL * p1=mPowderPatternComponentRegistry.GetObj(i)
6072 .mPowderPatternCalc.data();
6074 const REAL s = mScaleFactor(i);
6075 for(
unsigned long j=0;j<mNbPointUsed;j++) *p0++ += s * *p1++;
6077 TAU_PROFILE_STOP (timer2);
6081 TAU_PROFILE_START(timer3);
6084 const REAL * p1=mPowderPatternComponentRegistry.GetObj(i)
6085 .mPowderPatternCalc.data();
6087 for(
unsigned long j=0;j<mNbPointUsed;j++) *p0++ = *p1++;
6092 const REAL * p1=mPowderPatternComponentRegistry.GetObj(i)
6093 .mPowderPatternCalc.data();
6095 for(
unsigned long j=0;j<mNbPointUsed;j++) *p0++ += *p1++;
6098 TAU_PROFILE_STOP(timer3);
6099 TAU_PROFILE_START(timer4);
6105 const REAL *p1=mPowderPatternComponentRegistry.GetObj(i)
6106 .mPowderPatternCalc.data();
6107 for(
unsigned long j=0;j<mNbPointUsed;j++) *p0++ = *p1++;
6113 const REAL *p1=mPowderPatternComponentRegistry.GetObj(i)
6114 .mPowderPatternCalc.data();
6115 for(
unsigned long j=0;j<mNbPointUsed;j++) *p0++ += *p1++;
6118 TAU_PROFILE_STOP(timer4);
6123 TAU_PROFILE_START(timer5);
6126 VFN_DEBUG_MESSAGE(
"PowderPattern::CalcPowderPattern():variance",2);
6129 for(
int i=0;i<mPowderPatternComponentRegistry.GetNb();i++)
6131 if(mPowderPatternComponentRegistry.GetObj(i).HasPowderPatternCalcVariance())
6133 if(0==mPowderPatternComponentRegistry.GetObj(i).GetPowderPatternCalcVariance()
6134 .numElements())
break;
6137 const REAL *p1=mPowderPatternComponentRegistry.GetObj(i)
6138 .GetPowderPatternCalcVariance().data();
6140 if(
true==mPowderPatternComponentRegistry.GetObj(i).IsScalable())
6142 const REAL s2 = mScaleFactor(i) * mScaleFactor(i);
6143 for(
unsigned long j=0;j<mNbPointUsed;j++) *p0++ += s2 * *p1++;
6145 else for(
unsigned long j=0;j<mNbPointUsed;j++) *p0++ += *p1++;
6150 for(
unsigned long j=0;j<mNbPointUsed;j++)
6151 if(*p1 <=0) {*p0++ =0;p1++;}
6152 else *p0++ = 1. / *p1++;
6155 TAU_PROFILE_STOP(timer5);
6156 VFN_DEBUG_EXIT(
"PowderPattern::CalcPowderPattern():End",3);
6159 void PowderPattern::CalcPowderPattern_FullDeriv(std::set<RefinablePar*> &vPar)
6161 TAU_PROFILE(
"PowderPattern::CalcPowderPattern_FullDeriv()",
"void ()",TAU_DEFAULT);
6163 mPowderPattern_FullDeriv.clear();
6164 if(mPowderPatternComponentRegistry.GetNb()==0)
return;
6165 std::vector<std::map<RefinablePar*,CrystVector_REAL>*> comps;
6166 for(
int i=0;i<mPowderPatternComponentRegistry.GetNb();i++)
6167 comps.push_back(&(mPowderPatternComponentRegistry.GetObj(i).GetPowderPattern_FullDeriv(vPar)));
6169 for(std::set<RefinablePar *>::iterator par=vPar.begin();par!=vPar.end();++par)
6171 if(*par==0)
continue;
6172 for(
int i=0;i<mPowderPatternComponentRegistry.GetNb();i++)
6174 if((*par)->GetPointer()==mScaleFactor.data()+i)
6176 mPowderPattern_FullDeriv[*par]=mPowderPatternComponentRegistry.GetObj(i).GetPowderPatternCalc();
6181 if((*(comps[i]))[*par].size()==0)
continue;
6183 if(
true==mPowderPatternComponentRegistry.GetObj(i).IsScalable())
6185 if(mPowderPattern_FullDeriv[*par].size()==0)
6187 mPowderPattern_FullDeriv[*par].resize(
mNbPoint);
6188 const REAL * p1=(*comps[i])[*par].data();
6189 REAL * p0 = mPowderPattern_FullDeriv[*par].data();
6190 const REAL s = mScaleFactor(i);
6191 for(
unsigned long j=0;j<mNbPointUsed;j++) *p0++ = s * *p1++;
6192 for(
unsigned long j=mNbPointUsed;j<
mNbPoint;j++) *p0++ = 0;
6196 const REAL * p1=(*comps[i])[*par].data();
6197 REAL * p0 = mPowderPattern_FullDeriv[*par].data();
6198 const REAL s = mScaleFactor(i);
6199 for(
unsigned long j=0;j<mNbPointUsed;j++) *p0++ += s * *p1++;
6204 if(mPowderPattern_FullDeriv[*par].size()==0)
6206 mPowderPattern_FullDeriv[*par].resize(
mNbPoint);
6207 const REAL * p1=(*comps[i])[*par].data();
6208 REAL * p0 = mPowderPattern_FullDeriv[*par].data();
6209 for(
unsigned long j=0;j<mNbPointUsed;j++) *p0++ = *p1++;
6210 for(
unsigned long j=mNbPointUsed;j<
mNbPoint;j++) *p0++ = 0;
6214 const REAL * p1=(*comps[i])[*par].data();
6215 REAL * p0 = mPowderPattern_FullDeriv[*par].data();
6216 for(
unsigned long j=0;j<mNbPointUsed;j++) *p0++ += *p1++;
6223 std::map<RefinablePar*, CrystVector_REAL> oldDeriv;
6224 std::vector<const CrystVector_REAL*> v;
6226 for(std::map<RefinablePar*, CrystVector_REAL>::reverse_iterator pos=mPowderPattern_FullDeriv.rbegin();pos!=mPowderPattern_FullDeriv.rend();++pos)
6228 if(pos->first==0)
continue;
6229 if(pos->second.size()==0)
continue;
6231 const REAL step=pos->first->GetDerivStep();
6232 pos->first->Mutate(step);
6235 pos->first->Mutate(-2*step);
6238 oldDeriv[pos->first]/=2*step;
6239 pos->first->Mutate(step);
6241 v.push_back(&(pos->second));
6242 v.push_back(&(oldDeriv[pos->first]));
6243 cout<<pos->first->GetName()<<
":"<<pos->second.size()<<
","<<oldDeriv[pos->first].size()<<endl;
6246 cout<<FormatVertVector<REAL>(v,16)<<endl;
6257 TAU_PROFILE(
"PowderPattern::CalcPowderPatternIntegrated()",
"void ()",TAU_DEFAULT);
6258 VFN_DEBUG_ENTRY(
"PowderPattern::CalcPowderPatternIntegrated()",4);
6259 if(mPowderPatternComponentRegistry.GetNb()==0)
6263 VFN_DEBUG_EXIT(
"PowderPattern::CalcPowderPatternIntegrated():no components!",4);
6266 TAU_PROFILE_TIMER(timer1,
"PowderPattern::CalcPowderPatternIntegrated()1:Calc components",\
6268 TAU_PROFILE_TIMER(timer2,
"PowderPattern::CalcPowderPatternIntegrated()2:Add comps-scaled"\
6270 TAU_PROFILE_TIMER(timer3,
"PowderPattern::CalcPowderPatternIntegrated()2:Add backgd1"\
6272 TAU_PROFILE_TIMER(timer4,
"PowderPattern::CalcPowderPatternIntegrated()2:Add backgd2"\
6274 TAU_PROFILE_TIMER(timer5,
"PowderPattern::CalcPowderPatternIntegrated()3:Variance"
6276 TAU_PROFILE_START(timer1);
6277 vector< pair<const CrystVector_REAL*,const RefinableObjClock*> > comps;
6278 for(
int i=0;i<mPowderPatternComponentRegistry.GetNb();i++)
6280 comps.push_back(mPowderPatternComponentRegistry.GetObj(i).
6281 GetPowderPatternIntegratedCalc());
6283 TAU_PROFILE_STOP(timer1);
6290 for(vector< pair<const CrystVector_REAL*,const RefinableObjClock*> >::iterator
6291 pos=comps.begin();pos!=comps.end();++pos)
6300 VFN_DEBUG_EXIT(
"PowderPattern::CalcPowderPatternIntegrated():no need to recalc",4);
6303 VFN_DEBUG_MESSAGE(
"PowderPattern::CalcPowderPatternIntegrated():Recalc",3);
6306 for(
int i=0;i<mPowderPatternComponentRegistry.GetNb();i++)
6308 VFN_DEBUG_MESSAGE(
"PowderPattern::CalcPowderPatternIntegrated():Adding "\
6309 << mPowderPatternComponentRegistry.GetObj(i).GetName(),3);
6310 if(
true==mPowderPatternComponentRegistry.GetObj(i).IsScalable())
6312 TAU_PROFILE_START(timer2);
6315 const REAL * RESTRICT p1= comps[i].first->data();
6317 const REAL s = mScaleFactor(i);
6318 for(
unsigned long j=mNbIntegrationUsed;j>0;j--) *p0++ = s * *p1++;
6322 const REAL * RESTRICT p1= comps[i].first->data();
6324 const REAL s = mScaleFactor(i);
6325 for(
unsigned long j=mNbIntegrationUsed;j>0;j--) *p0++ += s * *p1++;
6327 TAU_PROFILE_STOP (timer2);
6331 TAU_PROFILE_START(timer3);
6334 const REAL * RESTRICT p1= comps[i].first->data();
6336 for(
unsigned long j=mNbIntegrationUsed;j>0;j--) *p0++ = *p1++;
6340 const REAL * RESTRICT p1= comps[i].first->data();
6342 for(
unsigned long j=mNbIntegrationUsed;j>0;j--) *p0++ += *p1++;
6345 TAU_PROFILE_STOP(timer3);
6346 TAU_PROFILE_START(timer4);
6351 const REAL * RESTRICT p1= comps[i].first->data();
6353 for(
unsigned long j=mNbIntegrationUsed;j>0;j--) *p0++ = *p1++;
6357 const REAL * RESTRICT p1= comps[i].first->data();
6359 for(
unsigned long j=mNbIntegrationUsed;j>0;j--) *p0++ += *p1++;
6362 TAU_PROFILE_STOP(timer4);
6365 TAU_PROFILE_START(timer5);
6368 bool useCalcVariance=
false;
6369 for(
int i=0;i<mPowderPatternComponentRegistry.GetNb();i++)
6370 if(mPowderPatternComponentRegistry.GetObj(i).HasPowderPatternCalcVariance())
6371 if(mPowderPatternComponentRegistry.GetObj(i)
6372 .GetPowderPatternIntegratedCalcVariance().first->numElements() !=0)
6373 useCalcVariance=
true;
6376 VFN_DEBUG_MESSAGE(
"PowderPattern::CalcPowderPatternIntegrated():variance",3);
6379 mIntegratedWeight.resize(mNbIntegrationUsed);
6380 const REAL * RESTRICT p1= mIntegratedVarianceObs.data();
6382 for(
unsigned long j=mNbIntegrationUsed;j>0;j--) *p0++ = *p1++;
6386 for(
int i=0;i<mPowderPatternComponentRegistry.GetNb();i++)
6388 if(mPowderPatternComponentRegistry.GetObj(i).HasPowderPatternCalcVariance())
6390 if(0==mPowderPatternComponentRegistry.GetObj(i)
6391 .GetPowderPatternIntegratedCalcVariance().first->numElements())
break;
6393 const REAL * RESTRICT p1= mPowderPatternComponentRegistry.GetObj(i)
6394 .GetPowderPatternIntegratedCalcVariance().first->data();
6399 if(
true==mPowderPatternComponentRegistry.GetObj(i).IsScalable())
6401 const REAL s2 = mScaleFactor(i) * mScaleFactor(i);
6402 for(
unsigned long j=mNbIntegrationUsed;j>0;j--) *p0++ += s2 * *p1++;
6404 else for(
unsigned long j=mNbIntegrationUsed;j>0;j--) *p0++ += *p1++;
6409 REAL *p0 = mIntegratedWeight.data();
6411 for(
unsigned long j=mNbIntegrationUsed;j>0;j--)
6412 if(*p1 <=0) {*p0++ =0;p1++;}
6413 else *p0++ = 1. / *p1++;
6415 else mIntegratedWeight.resize(0);
6417 TAU_PROFILE_STOP(timer5);
6443 VFN_DEBUG_EXIT(
"PowderPattern::CalcPowderPatternIntegrated():End",4);
6446 void PowderPattern::CalcPowderPatternIntegrated_FullDeriv(std::set<RefinablePar *> &vPar)
6448 TAU_PROFILE(
"PowderPattern::CalcPowderPatternIntegrated_FullDeriv()",
"void ()",TAU_DEFAULT);
6453 mPowderPatternUsed_FullDeriv.clear();
6454 if(mPowderPatternComponentRegistry.GetNb()==0)
return;
6455 std::vector<map<RefinablePar*,CrystVector_REAL>*> comps;
6456 for(
int i=0;i<mPowderPatternComponentRegistry.GetNb();i++)
6458 comps.push_back(&(mPowderPatternComponentRegistry.GetObj(i).
6459 GetPowderPatternIntegrated_FullDeriv(vPar)));
6463 mPowderPatternIntegrated_FullDeriv.clear();
6464 for(std::set<RefinablePar *>::iterator par=vPar.begin();par!=vPar.end();++par)
6466 if(*par==0)
continue;
6467 for(
int i=0;i<mPowderPatternComponentRegistry.GetNb();i++)
6469 if((*par)->GetPointer()==mScaleFactor.data()+i)
6472 mPowderPatternIntegrated_FullDeriv[*par]=*(mPowderPatternComponentRegistry.GetObj(i).GetPowderPatternIntegratedCalc().first);
6478 if((*(comps[i]))[*par].size()==0)
continue;
6480 if(
true==mPowderPatternComponentRegistry.GetObj(i).IsScalable())
6482 if(mPowderPatternIntegrated_FullDeriv[*par].size()==0)
6484 mPowderPatternIntegrated_FullDeriv[*par].resize(mNbIntegrationUsed);
6485 const REAL * RESTRICT p1= (*(comps[i]))[*par].data();
6486 REAL * RESTRICT p0 = mPowderPatternIntegrated_FullDeriv[*par].data();
6487 const REAL s = mScaleFactor(i);
6488 for(
unsigned long j=mNbIntegrationUsed;j>0;j--)
6494 if((j==mNbIntegrationUsed)&&((*par)->GetName()==
"Cimetidine_C11_x")) cout<<__FILE__<<
":"<<__LINE__<<
":"<<*p0<<
","<<s<<
"*"<<*p1<<endl;
6501 const REAL * RESTRICT p1= (*(comps[i]))[*par].data();
6502 REAL * RESTRICT p0 = mPowderPatternIntegrated_FullDeriv[*par].data();
6503 const REAL s = mScaleFactor(i);
6504 for(
unsigned long j=mNbIntegrationUsed;j>0;j--)
6510 if((j==mNbIntegrationUsed)&&((*par)->GetName()==
"Cimetidine_C11_x")) cout<<__FILE__<<
":"<<__LINE__<<
":"<<*p0<<
","<<s<<
"*"<<*p1<<endl;
6518 if(mPowderPatternIntegrated_FullDeriv[*par].size()==0)
6520 mPowderPatternIntegrated_FullDeriv[*par].resize(mNbIntegrationUsed);
6521 const REAL * RESTRICT p1= (*(comps[i]))[*par].data();
6522 REAL * RESTRICT p0 = mPowderPatternIntegrated_FullDeriv[*par].data();
6523 for(
unsigned long j=mNbIntegrationUsed;j>0;j--)
6529 if((j==mNbIntegrationUsed)&&((*par)->GetName()==
"Cimetidine_C11_x")) cout<<__FILE__<<
":"<<__LINE__<<
":"<<*p0<<
","<<*p1<<endl;
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 if((*par)->GetName()==
"Cimetidine_C11_x")
6552 cout<<
"PowderPattern::CalcPowderPatternIntegrated_FullDeriv():"
6553 <<(*par)->GetName()<<
":s="<<mScaleFactor(i)<<
", d[0]="<<(*(comps[i]))[*par](0)
6554 <<
", integ="<<mPowderPatternIntegrated_FullDeriv[*par](0)<<endl;
6568 REFPAR_DERIV_STEP_ABSOLUTE,
true,
true,
true,
false,RAD2DEG);
6575 REFPAR_DERIV_STEP_ABSOLUTE,
true,
true,
true,
false,RAD2DEG);
6582 REFPAR_DERIV_STEP_ABSOLUTE,
true,
true,
true,
false,RAD2DEG);
6589 REFPAR_DERIV_STEP_ABSOLUTE,
true,
true,
true,
false,1.0);
6596 REFPAR_DERIV_STEP_ABSOLUTE,
true,
true,
true,
false,1.0);
6603 REFPAR_DERIV_STEP_ABSOLUTE,
true,
true,
true,
false,1.0);
6611 bool needPrep=
false;
6612 for(
int i=0;i<mPowderPatternComponentRegistry.GetNb();i++)
6614 mPowderPatternComponentRegistry.GetObj(i).GetBraggLimits();
6615 if(mPowderPatternComponentRegistry.GetObj(i).GetClockBraggLimits()
6616 >mClockIntegratedFactorsPrep)
6625 if(mClockIntegratedFactorsPrep<mClockNbPointUsed) needPrep=
true;
6627 if(
false==needPrep)
return;
6628 VFN_DEBUG_ENTRY(
"PowderPattern::PrepareIntegratedRfactor()",3);
6629 TAU_PROFILE(
"PowderPattern::PrepareIntegratedRfactor()",
"void ()",TAU_DEFAULT);
6635 for(
int i=0;i<mPowderPatternComponentRegistry.GetNb();i++)
6637 const CrystVector_long vLim=mPowderPatternComponentRegistry.GetObj(i).GetBraggLimits();
6638 for(
int j=0;j<vLim.numElements();j++) vLimits.push_back(vLim(j));
6640 if(vLimits.size()<2)
6642 mIntegratedPatternMin.resize(0);
6643 mIntegratedPatternMax.resize(0);
6644 mNbIntegrationUsed=0;
6645 mClockIntegratedFactorsPrep.Click();
6647 VFN_DEBUG_EXIT(
"PowderPattern::PrepareIntegratedRfactor(): no intervals !",3);
6650 if(*(vLimits.begin())<0)
6652 vLimits.push_back(0);
6655 for(list<long>::iterator pos=vLimits.begin();pos!=vLimits.end();)
6657 if( (*pos<0) || (*pos>=
long(mNbPointUsed)) ) pos=vLimits.erase(pos);
6662 list<long> vLimits2;
6663 list<long>::iterator pos1=vLimits.begin();
6664 list<long>::iterator pos2=pos1;pos2++;
6665 for(;pos2!=vLimits.end();)
6667 const long pix1=*pos1;
6669 vLimits2.push_back(pix1);
6673 if(pos2==vLimits.end())
break;
6674 if(*pos2>(pix1+8))
break;
6677 vLimits2.push_back(*pos1);
6680 pos1=vLimits2.begin();
6682 for(;pos2!=vLimits2.end();)
6685 if( *pos2<=((*pos1)+2))
6688 pos2=vLimits2.erase(pos2);
6691 else {pos1++;pos2++;}
6695 list<pair<long,long> > vLimits3;
6696 pos1=vLimits2.begin();
6698 for(;pos2!=vLimits2.end();)
6700 if(*pos2!=(
long(mNbPointUsed)-1)) vLimits3.push_back(make_pair(*pos1++,*pos2++-1));
6701 else vLimits3.push_back(make_pair(*pos1++,*pos2++));
6705 mIntegratedPatternMin.resize(vLimits3.size());
6706 mIntegratedPatternMax.resize(vLimits3.size());
6708 for(list<pair<long,long> >::iterator pos=vLimits3.begin();pos!=vLimits3.end();++pos)
6710 mIntegratedPatternMin(i)=pos->first;
6711 mIntegratedPatternMax(i++)=pos->second;
6714 long numInterval=mIntegratedPatternMin.numElements();
6715 CrystVector_bool keep(numInterval);
6722 VFN_DEBUG_MESSAGE(
"PowderPattern::PrepareIntegratedRfactor():5:Excluded regions("<<nbExclude<<
")",3);
6724 long minExcl,maxExcl;
6727 for(
int i=0;i<nbExclude;i++)
6729 while(mIntegratedPatternMax(j)<minExcl)
6732 if(j>=numInterval)
break;
6734 if(j>=numInterval)
break;
6735 while(mIntegratedPatternMin(j)<maxExcl)
6737 if( (mIntegratedPatternMin(j)>minExcl) &&(mIntegratedPatternMax(j)<maxExcl))
6739 if( (mIntegratedPatternMin(j)<minExcl) &&(mIntegratedPatternMax(j)<maxExcl))
6740 mIntegratedPatternMax(j)=minExcl;
6741 if( (mIntegratedPatternMin(j)>minExcl) &&(mIntegratedPatternMax(j)>maxExcl))
6742 mIntegratedPatternMin(j)=maxExcl;
6743 if(j==(numInterval-1))
break;
6750 while(mIntegratedPatternMax(j)>=minExcl)
6758 VFN_DEBUG_MESSAGE(
"PowderPattern::PrepareIntegratedRfactor():6",3);
6760 for(
int i=0;i<numInterval;i++)
6764 mIntegratedPatternMin(j )=mIntegratedPatternMin(i);
6765 mIntegratedPatternMax(j++)=mIntegratedPatternMax(i);
6769 mIntegratedPatternMax.resizeAndPreserve(numInterval);
6770 mIntegratedPatternMin.resizeAndPreserve(numInterval);
6772 VFN_DEBUG_MESSAGE(
"PowderPattern::PrepareIntegratedRfactor():intervals"<<endl\
6775 mIntegratedObs.resize(numInterval);
6776 mIntegratedVarianceObs.resize(numInterval);
6777 mIntegratedVarianceObs=0;
6779 mIntegratedWeightObs.resize(numInterval);
6780 for(
int i=0;i<numInterval;i++)
6782 for(
int j=mIntegratedPatternMin(i);j<=mIntegratedPatternMax(i);j++)
6787 if(mIntegratedVarianceObs(i) <= 0) mIntegratedWeightObs(i)=0;
6788 else mIntegratedWeightObs(i)=1./mIntegratedVarianceObs(i);
6795 mNbIntegrationUsed=mIntegratedPatternMin.numElements();
6796 mClockIntegratedFactorsPrep.Click();
6797 VFN_DEBUG_EXIT(
"PowderPattern::PrepareIntegratedRfactor()",3);
6804 if(this->
GetRadiation().GetWavelengthType()==WAVELENGTH_TOF)
6811 if(1>fabs(sinth)) tmp=(
unsigned long)(this->
X2PixelCorr(2*asin(sinth)));
else tmp=
mNbPoint;
6814 if(tmp !=mNbPointUsed)
6817 mClockNbPointUsed.Click();
6818 VFN_DEBUG_MESSAGE(
"PowderPattern::CalcNbPointUsed():"<<mNbPointUsed<<
" max(sin(theta)/lambda)="<<
mMaxSinThetaOvLambda, 3)
6825 VFN_DEBUG_MESSAGE(
"PowderPattern::InitOptions()",5)
6826 static string OptProfileIntegrationName;
6827 static string OptProfileIntegrationChoices[2];
6829 static bool needInitNames=
true;
6830 if(
true==needInitNames)
6832 OptProfileIntegrationName=
"Use Integrated Profiles";
6833 OptProfileIntegrationChoices[0]=
"Yes (recommended)";
6834 OptProfileIntegrationChoices[1]=
"No";
6836 needInitNames=
false;
6838 mOptProfileIntegration.Init(2,&OptProfileIntegrationName,OptProfileIntegrationChoices);
6839 this->
AddOption(&mOptProfileIntegration);
6842 #ifdef __WX__CRYST__
6847 return mpWXCrystObj;
6860 SPGScore::SPGScore(
const string &s,
const REAL r,
const REAL g,
const unsigned int nbextinct,
const REAL ngof,
const unsigned int nbrefl):
6861 hm(s),rw(r),gof(g),ngof(ngof),nbextinct446(nbextinct),nbreflused(nbrefl)
6867 if(s1.ngof > 0.00001 && s2.ngof >0.00001)
return s1.ngof < s2.ngof;
6868 return s1.gof < s2.gof;
6876 std::vector<bool> fingerprint(5*5*7-1+6);
6884 for(
long h=0;h<5;++h)
6885 for(
long k=0;k<5;++k)
6886 for (
long l=0;l<7;++l)
6888 if((h+k+l)==0)
continue;
6889 cctbx::miller::index<long> hkl(scitbx::vec3<long>(h,k,l));
6890 if(i>=fingerprint.size()) cout<<
"WHOOOOOOOOOOOOOPS"<<endl;
6891 fingerprint[i++] =spg.is_sys_absent(hkl);
6899 SpaceGroupExplorer::SpaceGroupExplorer(PowderPatternDiffraction *pd):
6902 SPGScore SpaceGroupExplorer::Run(
const string &spgId,
const bool fitprofile,
6903 const bool verbose,
const bool restore_orig,
6904 const bool update_display,
const REAL relative_length_tolerance,
6905 const REAL absolute_angle_tolerance_degree)
6907 cctbx::sgtbx::space_group sg;
6911 cctbx::sgtbx::space_group_symbols sgs=cctbx::sgtbx::space_group_symbols(spgId);
6912 sg = cctbx::sgtbx::space_group(sgs);
6914 catch(exception ex1)
6919 sg = cctbx::sgtbx::space_group(spgId);
6921 catch(exception ex2)
6923 string emsg =
"Space group symbol '" + spgId +
"' not recognized";
6924 throw ObjCrystException(emsg);
6927 return this->Run(sg, fitprofile, verbose, restore_orig, update_display,
6928 relative_length_tolerance, absolute_angle_tolerance_degree);
6931 SPGScore SpaceGroupExplorer::Run(
const cctbx::sgtbx::space_group &spg,
const bool fitprofile,
const bool verbose,
6932 const bool restore_orig,
const bool update_display,
6933 const REAL relative_length_tolerance,
const REAL absolute_angle_tolerance_degree)
6935 TAU_PROFILE(
"SpaceGroupExplorer::Run()",
"void (wxCommandEvent &)",TAU_DEFAULT);
6936 TAU_PROFILE_TIMER(timer1,
"SpaceGroupExplorer::Run()LSQ-P1",
"", TAU_FIELD);
6937 TAU_PROFILE_TIMER(timer2,
"SpaceGroupExplorer::Run()LSQ1",
"", TAU_FIELD);
6938 TAU_PROFILE_TIMER(timer3,
"SpaceGroupExplorer::Run()LSQ2",
"", TAU_FIELD);
6939 Crystal *pCrystal=&(mpDiff->GetCrystal());
6941 const REAL a=pCrystal->GetLatticePar(0),
6942 b=pCrystal->GetLatticePar(1),
6943 c=pCrystal->GetLatticePar(2),
6944 d=pCrystal->GetLatticePar(3),
6945 e=pCrystal->GetLatticePar(4),
6946 f=pCrystal->GetLatticePar(5);
6947 const string spghm=pCrystal->GetSpaceGroup().GetCCTbxSpg().match_tabulated_settings().hermann_mauguin();
6948 const string name=pCrystal->GetName();
6950 const cctbx::sgtbx::space_group_symbols s = spg.match_tabulated_settings();
6951 const string hm=s.universal_hermann_mauguin();
6952 const cctbx::uctbx::unit_cell uc(scitbx::af::double6(a,b,c,d*RAD2DEG,e*RAD2DEG,f*RAD2DEG));
6953 if(!spg.is_compatible_unit_cell(uc,relative_length_tolerance,absolute_angle_tolerance_degree))
6955 throw ObjCrystException(
"Spacegroup is not compatible with unit cell.");
6957 mpDiff->GetCrystal().Init(a,b,c,d,e,f,hm,name);
6958 mpDiff->SetExtractionMode(
true,
true);
6959 unsigned int nbcycle=1;
6960 if(update_display) mpDiff->GetParentPowderPattern().UpdateDisplay();
6962 unsigned int nbfreepar=mpDiff->GetProfileFitNetNbObs();
6963 if(nbfreepar<1) nbfreepar=1;
6967 lsq.SetRefinedObj(mpDiff->GetParentPowderPattern(),0,
true,
true);
6968 lsq.PrepareRefParList(
true);
6969 const unsigned int saved_par = lsq.GetCompiledRefinedObj().CreateParamSet(
"SpaceGroupExplorer saved parameters");
6970 lsq.GetCompiledRefinedObj().SaveParamSet(saved_par);
6973 for(
unsigned int j=0;j<nbcycle;j++)
6976 mpDiff->SetExtractionMode(
true,
true);
6977 const float t0=chrono.seconds();
6978 if(verbose) cout<<
"Doing Le Bail, t="<<
FormatFloat(t0,6,2)<<
"s";
6979 mpDiff->ExtractLeBail(5);
6980 if(verbose) cout<<
", dt="<<
FormatFloat(chrono.seconds()-t0,6,2)<<
"s"<<endl;
6984 TAU_PROFILE_START(timer2);
6985 lsq.SetParIsFixed(gpRefParTypeObjCryst,
true);
6986 lsq.SetParIsFixed(gpRefParTypeScattDataScale,
false);
6987 std::list<RefinablePar*> vnewpar;
6988 std::list<const RefParType*> vnewpartype;
6990 if(s.number()==1) vnewpar.push_back(&lsq.GetCompiledRefinedObj().GetPar(
"Zero"));
6991 lsq.SetParIsFixed(gpRefParTypeUnitCell,
false);
6992 lsq.SafeRefine(vnewpar, vnewpartype,1.01,2,
true,
true);
6994 TAU_PROFILE_STOP(timer2);
6997 TAU_PROFILE_START(timer1);
6998 vnewpar.push_back(&lsq.GetCompiledRefinedObj().GetPar(
"2ThetaDispl"));
6999 vnewpar.push_back(&lsq.GetCompiledRefinedObj().GetPar(
"2ThetaTransp"));
7000 lsq.SafeRefine(vnewpar, vnewpartype,1.01,2,
true,
true);
7002 lsq.SetParIsFixed(gpRefParTypeScattDataBackground,
false);
7004 const unsigned int nbcomp= mpDiff->GetParentPowderPattern().GetNbPowderPatternComponent();
7005 for(
unsigned int i=0;i<nbcomp;++i)
7006 if(mpDiff->GetParentPowderPattern().GetPowderPatternComponent(i).GetClassName()==
"PowderPatternBackground")
7008 PowderPatternBackground *pback=
dynamic_cast<PowderPatternBackground *
>
7009 (&(mpDiff->GetParentPowderPattern().GetPowderPatternComponent(i)));
7010 pback->FixParametersBeyondMaxresolution(lsq.GetCompiledRefinedObj());
7012 for(
unsigned int i=0; i<lsq.GetCompiledRefinedObj().GetNbPar();i++)
7013 if( (lsq.GetCompiledRefinedObj().GetPar(i).IsFixed()==
false)
7015 vnewpar.push_back(&lsq.GetCompiledRefinedObj().GetPar(i));
7016 lsq.SafeRefine(vnewpar,vnewpartype,1.01,2,
true,
true);
7018 TAU_PROFILE_STOP(timer1);
7021 mpDiff->SetExtractionMode(
true,
true);
7022 mpDiff->ExtractLeBail(5);
7023 TAU_PROFILE_START(timer3);
7024 lsq.SafeRefine(vnewpar,vnewpartype,1.01,3,
true,
true);
7025 TAU_PROFILE_STOP(timer3);
7029 if(update_display) mpDiff->GetParentPowderPattern().UpdateDisplay();
7030 const REAL rw=mpDiff->GetParentPowderPattern().GetRw()*100;
7031 const REAL gof=mpDiff->GetParentPowderPattern().GetChi2()/mpDiff->GetParentPowderPattern().GetNbPointUsed();
7032 if(verbose) cout << boost::format(
" (cycle #%u)\n Rwp=%5.2f%%\n GoF=%9.2f") % j % rw % gof <<endl;
7034 const REAL rw=mpDiff->GetParentPowderPattern().GetRw()*100;
7035 const REAL gof=mpDiff->GetParentPowderPattern().GetChi2()/nbfreepar;
7036 const REAL ngof = this->GetP1IntegratedGoF();
7037 unsigned int nbextinct446=0;
7039 for(
unsigned int i=6;i<fgp.size();++i) nbextinct446+=(
unsigned int)(fgp[i]);
7040 const unsigned int nbrefl = mpDiff->GetNbReflBelowMaxSinThetaOvLambda();
7041 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;
7045 mpDiff->GetCrystal().Init(a,b,c,d,e,f,spghm,name);
7046 lsq.GetCompiledRefinedObj().RestoreParamSet(saved_par);
7049 mpDiff->GetCrystal().UpdateDisplay();
7050 this->mpDiff->GetParentPowderPattern().UpdateDisplay();
7054 return SPGScore(hm.c_str(),rw,gof,nbextinct446, ngof, nbrefl);
7057 void SpaceGroupExplorer::RunAll(
const bool fitprofile_all,
const bool verbose,
const bool keep_best,
7058 const bool update_display,
const bool fitprofile_p1,
7059 const REAL relative_length_tolerance,
const REAL absolute_angle_tolerance_degree)
7061 Crystal *pCrystal=&(mpDiff->GetCrystal());
7064 const REAL a=pCrystal->GetLatticePar(0),
7065 b=pCrystal->GetLatticePar(1),
7066 c=pCrystal->GetLatticePar(2),
7067 d=pCrystal->GetLatticePar(3),
7068 e=pCrystal->GetLatticePar(4),
7069 f=pCrystal->GetLatticePar(5);
7070 const cctbx::uctbx::unit_cell uc(scitbx::af::double6(a,b,c,d*RAD2DEG,e*RAD2DEG,f*RAD2DEG));
7071 const string spghm=pCrystal->GetSpaceGroup().GetCCTbxSpg().match_tabulated_settings().hermann_mauguin();
7072 const string name=pCrystal->GetName();
7074 cctbx::sgtbx::space_group_symbol_iterator it=cctbx::sgtbx::space_group_symbol_iterator();
7076 unsigned int nbspg=0;
7079 cctbx::sgtbx::space_group_symbols s=it.next();
7080 if(s.number()==0)
break;
7081 cctbx::sgtbx::space_group spg(s);
7082 if(spg.is_compatible_unit_cell(uc,relative_length_tolerance, absolute_angle_tolerance_degree)) nbspg++;
7085 if(verbose) cout << boost::format(
"Beginning spacegroup exploration... %u to go...\n") % nbspg;
7088 mvSPGExtinctionFingerprint.clear();
7091 unsigned int nb_refl_p1=1;
7093 it=cctbx::sgtbx::space_group_symbol_iterator();
7098 cctbx::sgtbx::space_group_symbols s=it.next();
7099 if(s.number()==0)
break;
7100 cctbx::sgtbx::space_group spg(s);
7101 bool compat=spg.is_compatible_unit_cell(uc,relative_length_tolerance,absolute_angle_tolerance_degree);
7105 const string hm=s.universal_hermann_mauguin();
7107 pCrystal->Init(a,b,c,d,e,f,hm,name);
7110 std::map<std::vector<bool>,SPGScore>::iterator posfgp=mvSPGExtinctionFingerprint.find(fgp);
7111 if(posfgp!=mvSPGExtinctionFingerprint.end())
7113 pCrystal->Init(a,b,c,d,e,f,hm,name);
7114 mpDiff->SetExtractionMode(
true,
true);
7115 unsigned int nbrefl = mpDiff->GetNbReflBelowMaxSinThetaOvLambda();
7116 REAL ngof = (posfgp->second.ngof * nbrefl) / posfgp->second.nbreflused;
7117 mvSPG.push_back(SPGScore(hm.c_str(),posfgp->second.rw,posfgp->second.gof,posfgp->second.nbextinct446, ngof, nbrefl));
7118 if(verbose) cout<<boost::format(
" (#%3d) %-14s: Rwp= %5.2f%% GoF=%9.2f nGoF=%9.2f (%3u reflections, %3u extinct)")
7119 % s.number() % hm.c_str() % mvSPG.back().rw % mvSPG.back().gof % mvSPG.back().ngof % mvSPG.back().nbreflused % mvSPG.back().nbextinct446
7120 <<
" [same extinctions as:"<<posfgp->second.hm<<
"]\n";
7124 if(((s.number()==1) && fitprofile_p1) || fitprofile_all) mvSPG.push_back(this->Run(spg,
true,
false,
false, update_display,
7125 relative_length_tolerance, absolute_angle_tolerance_degree));
7126 else mvSPG.push_back(this->Run(spg,
false,
false,
true, update_display,relative_length_tolerance,absolute_angle_tolerance_degree));
7127 if(s.number() == 1) nb_refl_p1 = mvSPG.back().nbreflused;
7128 mvSPG.back().ngof *= mpDiff->GetNbReflBelowMaxSinThetaOvLambda() / (float)nb_refl_p1;
7129 mvSPGExtinctionFingerprint.insert(make_pair(fgp, mvSPG.back()));
7131 if(verbose) cout<<boost::format(
" (#%3d) %-14s: Rwp= %5.2f%% GoF=%9.2f nGoF=%9.2f (%3u reflections, %3u extinct)\n")
7132 % s.number() % hm.c_str() % mvSPG.back().rw % mvSPG.back().gof % mvSPG.back().ngof % mvSPG.back().nbreflused % mvSPG.back().nbextinct446;
7136 mvSPG.sort(compareSPGScore);
7139 if(verbose) cout<<
"Restoring best spacegroup: "<<mvSPG.front().hm<<endl;
7140 pCrystal->ChangeSpaceGroup(mvSPG.front().hm);
7145 pCrystal->Init(a,b,c,d,e,f,spghm,name);
7147 mpDiff->SetExtractionMode(
true,
true);
7148 mpDiff->ExtractLeBail(5);
7151 pCrystal->UpdateDisplay();
7152 this->mpDiff->GetParentPowderPattern().UpdateDisplay();
7157 const list<SPGScore>& SpaceGroupExplorer::GetScores()
const
7162 REAL SpaceGroupExplorer::GetP1IntegratedGoF()
7164 if(mpDiff->GetCrystal().GetSpaceGroup().GetSpaceGroupNumber()==1)
7166 mpDiff->GetPowderPatternIntegratedCalc();
7167 mP1IntegratedProfileMin = mpDiff->GetParentPowderPattern().GetIntegratedProfileMin();
7168 mP1IntegratedProfileMax = mpDiff->GetParentPowderPattern().GetIntegratedProfileMax();
7172 else if (mP1IntegratedProfileMin.size()==0)
return 0;
7175 REAL integratedChi2=0.;
7176 REAL integratedChi2LikeNorm=0.;
7177 const REAL * RESTRICT p1, * RESTRICT p2, * RESTRICT p3;
7178 CrystVector_REAL
const* pcalc = &(mpDiff->GetParentPowderPattern().GetPowderPatternCalc());
7179 CrystVector_REAL
const* pobs = &(mpDiff->GetParentPowderPattern().GetPowderPatternObs());
7180 CrystVector_REAL
const* psigma = &(mpDiff->GetParentPowderPattern().GetPowderPatternObsSigma());
7181 const unsigned int jmax = mpDiff->GetParentPowderPattern().GetNbPointUsed();
7183 unsigned int nbpoint = 0;
7184 for(
unsigned long i=0;i<mP1IntegratedProfileMin.size();i++)
7186 if(mP1IntegratedProfileMin(i) > jmax)
break;
7187 if(mP1IntegratedProfileMax(i) < 0)
continue;
7189 for(
unsigned long j=mP1IntegratedProfileMin(i); j<=mP1IntegratedProfileMax(i); j++)
7192 if(j >= jmax)
break;
7196 v += (*psigma)(j)*(*psigma)(j);
7198 if(v>0) chi2 += (c-o)*(c-o)/v;
7200 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.
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.