FOX/ObjCryst++  1.10.X (development)
PowderPattern.cpp
1 /* ObjCryst++ Object-Oriented Crystallographic Library
2  (c) 2000-2002 Vincent Favre-Nicolin vincefn@users.sourceforge.net
3  2000-2001 University of Geneva (Switzerland)
4 
5  This program is free software; you can redistribute it and/or modify
6  it under the terms of the GNU General Public License as published by
7  the Free Software Foundation; either version 2 of the License, or
8  (at your option) any later version.
9 
10  This program is distributed in the hope that it will be useful,
11  but WITHOUT ANY WARRANTY; without even the implied warranty of
12  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13  GNU General Public License for more details.
14 
15  You should have received a copy of the GNU General Public License
16  along with this program; if not, write to the Free Software
17  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
18 */
19 /*
20 * source file for LibCryst++ PowderPattern class
21 *
22 */
23 
24 #include <cstdlib>
25 
26 #include <typeinfo>
27 #include <stdio.h> //for sprintf()
28 
29 #include "cctbx/sgtbx/space_group.h" // For fullprof export
30 
31 #include "ObjCryst/ObjCryst/PowderPattern.h"
32 #include "ObjCryst/ObjCryst/Molecule.h" // For fullprof export
33 #include "ObjCryst/ObjCryst/PowderPatternBackgroundBayesianMinimiser.h"
34 #include "ObjCryst/RefinableObj/Simplex.h"
35 #include "ObjCryst/Quirks/VFNDebug.h"
36 #include "ObjCryst/Quirks/VFNStreamFormat.h"
37 #include "ObjCryst/ObjCryst/CIF.h"
38 #ifdef __WX__CRYST__
39  #include "ObjCryst/wxCryst/wxPowderPattern.h"
40 #endif
41 
42 #include <fstream>
43 #include <iomanip>
44 #include <sstream>
45 
46 #ifdef _MSC_VER // MS VC++ predefined macros....
47 #undef min
48 #undef max
49 #endif
50 
51 //#define USE_BACKGROUND_MAXLIKE_ERROR
52 
53 namespace ObjCryst
54 {
56 //
57 // PowderPatternComponent
58 //
60 ObjRegistry<PowderPatternComponent>
61  gPowderPatternComponentRegistry("List of all PowderPattern Components");
62 
63 PowderPatternComponent::PowderPatternComponent():
64 mIsScalable(false),mpParentPowderPattern(0)
65 {
66  gPowderPatternComponentRegistry.Register(*this);
67  mClockMaster.AddChild(mClockBraggLimits);
68 }
69 
70 PowderPatternComponent::PowderPatternComponent(const PowderPatternComponent &old):
71 mIsScalable(old.mIsScalable),
72 mpParentPowderPattern(old.mpParentPowderPattern)
73 {
74  mClockMaster.AddChild(mClockBraggLimits);
75  if(mpParentPowderPattern!=0)
76  {
77  mClockMaster.AddChild(mpParentPowderPattern->GetClockPowderPatternPar());
78  mClockMaster.AddChild(mpParentPowderPattern->GetClockPowderPatternXCorr());
79  mClockMaster.AddChild(mpParentPowderPattern->GetClockPowderPatternRadiation());
80  }
81 }
82 
83 PowderPatternComponent::~PowderPatternComponent()
84 {
85  gPowderPatternComponentRegistry.DeRegister(*this);
86 }
87 
88 const string& PowderPatternComponent::GetClassName() const
89 {
90  const static string className="PowderPatternComponent";
91  return className;
92 }
93 
94 const PowderPattern& PowderPatternComponent::GetParentPowderPattern()const
95 {
96  return *mpParentPowderPattern;
97 }
98 
99 std::map<RefinablePar*,CrystVector_REAL>& PowderPatternComponent::GetPowderPattern_FullDeriv(std::set<RefinablePar *> &vPar)
100 {
101  this->CalcPowderPattern_FullDeriv(vPar);
102  return mPowderPattern_FullDeriv;
103 }
104 
105 std::map<RefinablePar*,CrystVector_REAL>& PowderPatternComponent::GetPowderPatternIntegrated_FullDeriv(std::set<RefinablePar *> &vPar)
106 {
107  this->CalcPowderPatternIntegrated_FullDeriv(vPar);
108  return mPowderPatternIntegrated_FullDeriv;
109 }
110 
111 PowderPattern& PowderPatternComponent::GetParentPowderPattern()
112 {
113  return *mpParentPowderPattern;
114 }
115 
116 bool PowderPatternComponent::IsScalable()const {return mIsScalable;}
117 
118 const RefinableObjClock& PowderPatternComponent::GetClockPowderPatternCalc()const
119 {
120  return mClockPowderPatternCalc;
121 }
122 const RefinableObjClock& PowderPatternComponent::GetClockBraggLimits()const
123 {
124  return mClockBraggLimits;
125 }
126 
127 const list<pair<const REAL, const string> >& PowderPatternComponent::GetPatternLabelList()const
128 {
129  return mvLabel;
130 }
131 
132 void PowderPatternComponent::CalcPowderPattern_FullDeriv(std::set<RefinablePar *> &vPar)
133 {
134  TAU_PROFILE("PowderPatternComponent::CalcPowderPattern_FullDeriv()","void ()",TAU_DEFAULT);
135  mPowderPattern_FullDeriv.clear();
136  for(std::set<RefinablePar *>::iterator par=vPar.begin();par!=vPar.end();par++)
137  {
138  if(*par==0)
139  {
140  mPowderPattern_FullDeriv[*par]=this->GetPowderPatternCalc();
141  continue;
142  }
143  const REAL step=(*par)->GetDerivStep();
144  (*par)->Mutate(step);
145  mPowderPattern_FullDeriv[*par] =this->GetPowderPatternCalc();
146  (*par)->Mutate(-2*step);
147  mPowderPattern_FullDeriv[*par]-=this->GetPowderPatternCalc();
148  (*par)->Mutate(step);
149  mPowderPattern_FullDeriv[*par]/= step*2;
150  if(MaxAbs(mPowderPattern_FullDeriv[*par])==0)
151  mPowderPattern_FullDeriv[*par].resize(0);
152  }
153 }
154 
155 void PowderPatternComponent::CalcPowderPatternIntegrated_FullDeriv(std::set<RefinablePar *> &vPar)
156 {
157  TAU_PROFILE("PowderPatternComponent::CalcPowderPatternIntegrated_FullDeriv()","void ()",TAU_DEFAULT);
158  mPowderPatternIntegrated_FullDeriv.clear();
159  for(std::set<RefinablePar *>::iterator par=vPar.begin();par!=vPar.end();par++)
160  {
161  if(*par==0)
162  {
163  mPowderPatternIntegrated_FullDeriv[*par]=*(this->GetPowderPatternIntegratedCalc().first);
164  continue;
165  }
166  const REAL step=(*par)->GetDerivStep();
167  (*par)->Mutate(step);
168  mPowderPatternIntegrated_FullDeriv[*par] =*(this->GetPowderPatternIntegratedCalc().first);
169  (*par)->Mutate(-2*step);
170  mPowderPatternIntegrated_FullDeriv[*par]-=*(this->GetPowderPatternIntegratedCalc().first);
171  (*par)->Mutate(step);
172  mPowderPatternIntegrated_FullDeriv[*par]/= step*2;
173  if(MaxAbs(mPowderPatternIntegrated_FullDeriv[*par])==0)
174  mPowderPatternIntegrated_FullDeriv[*par].resize(0);
175  }
176 }
177 
179 //
180 // PowderPatternBackground
181 //
183 PowderPatternBackground::PowderPatternBackground():
184 mBackgroundNbPoint(0),
185 mMaxSinThetaOvLambda(10),mModelVariance(0)
186 {
187  mClockMaster.AddChild(mClockBackgroundPoint);
188  this->InitOptions();
189 }
190 
191 PowderPatternBackground::PowderPatternBackground(const PowderPatternBackground &old):
192 mBackgroundNbPoint(old.mBackgroundNbPoint),
193 mBackgroundInterpPointX(old.mBackgroundInterpPointX),
194 mBackgroundInterpPointIntensity(old.mBackgroundInterpPointIntensity),
195 mMaxSinThetaOvLambda(10),mModelVariance(0)
196 {
197  mClockMaster.AddChild(mClockBackgroundPoint);
198  this->InitOptions();
199  mInterpolationModel.SetChoice(old.mInterpolationModel.GetChoice());
200 }
201 
202 PowderPatternBackground::~PowderPatternBackground(){}
204 {
205  const static string className="PowderPatternBackground";
206  return className;
207 }
208 
210 {
211  if(mpParentPowderPattern!=0)
218 }
219 const CrystVector_REAL& PowderPatternBackground::GetPowderPatternCalc()const
220 {
221  this->CalcPowderPattern();
222  return mPowderPatternCalc;
223 }
224 
225 pair<const CrystVector_REAL*,const RefinableObjClock*>
227 {
228  VFN_DEBUG_MESSAGE("PowderPatternBackground::GetPowderPatternIntegratedCalc()",3)
229  this->CalcPowderPatternIntegrated();
231 }
232 
234 {
235  VFN_DEBUG_MESSAGE("PowderPatternBackground::ImportUserBackground():"<<filename,5)
236  ifstream fin (filename.c_str());
237  if(!fin)
238  {
239  throw ObjCrystException("PowderPatternBackground::ImportUserBackground() : \
240 Error opening file for input:"+filename);
241  }
242  long max=100;
243  CrystVector_REAL bckgd2Theta(max),bckgd(max);
244 
245  long nbPoints=0;
246  do
247  {
248  fin >> bckgd2Theta(nbPoints);
249  fin >> bckgd(nbPoints);
250  if(!fin) break;
251  VFN_DEBUG_MESSAGE("Background=" << bckgd(nbPoints)\
252  <<" at 2theta="<<bckgd2Theta(nbPoints),3)
253  nbPoints++;
254  if(nbPoints==max)
255  {
256  max += 100;
257  bckgd2Theta.resizeAndPreserve(max);
258  bckgd.resizeAndPreserve(max);
259  }
260  } while(fin.eof() == false);
261  fin.close();
262  bckgd2Theta.resizeAndPreserve(nbPoints);
263  bckgd.resizeAndPreserve(nbPoints);
264  if(mpParentPowderPattern!=0)
265  { if((this->GetParentPowderPattern().GetRadiation().GetWavelengthType()==WAVELENGTH_MONOCHROMATIC)
266  ||(this->GetParentPowderPattern().GetRadiation().GetWavelengthType()==WAVELENGTH_ALPHA12))
267  bckgd2Theta*= DEG2RAD;
268  } else bckgd2Theta*= DEG2RAD;
269 
270  this->SetInterpPoints(bckgd2Theta,bckgd);
271  this->InitRefParList();
273  {
274  char buf [200];
275  sprintf(buf,"Imported %d background points",(int)nbPoints);
276  (*fpObjCrystInformUser)((string)buf);
277  }
278  this->UpdateDisplay();
279  VFN_DEBUG_MESSAGE("PowderPatternBackground::ImportUserBackground():finished",5)
280 }
281 void PowderPatternBackground::SetInterpPoints(const CrystVector_REAL tth,
282  const CrystVector_REAL backgd)
283 {
284  VFN_DEBUG_ENTRY("PowderPatternBackground::SetInterpPoints():",5)
285  if( (tth.numElements()!=backgd.numElements())
286  ||(tth.numElements()<2))
287  {
288  throw ObjCrystException("PowderPatternBackground::SetInterpPoints() : \
289 number of points differ or less than 2 points !");
290  }
291  mBackgroundNbPoint=tth.numElements();
294  // Sort in ascending order, disregarding radiation type.
295  CrystVector<long> subs;
296  subs=SortSubs(tth);
297 
298  for(long i=0;i<mBackgroundNbPoint;++i)
299  {
300  mBackgroundInterpPointX(i)=tth(subs(i));
301  mBackgroundInterpPointIntensity(i)=backgd(subs(i));
302  }
303  this->InitRefParList();
305  VFN_DEBUG_EXIT("PowderPatternBackground::SetInterpPoints()",5)
306 }
307 
308 const pair<const CrystVector_REAL*,const CrystVector_REAL*> PowderPatternBackground::GetInterpPoints()const
309 {
311 }
312 
314  CrystVector_uint & groupIndex,
315  unsigned int &first) const
316 {
317  // One group for all background points
318  unsigned int index=0;
319  VFN_DEBUG_MESSAGE("PowderPatternBackground::GetGeneGroup()",4)
320  for(long i=0;i<obj.GetNbPar();i++)
321  for(long j=0;j<this->GetNbPar();j++)
322  if(&(obj.GetPar(i)) == &(this->GetPar(j)))
323  {
324  if(index==0) index=first++;
325  groupIndex(i)=index;
326  }
327 }
328 
329 void PowderPatternBackground::BeginOptimization(const bool allowApproximations,
330  const bool enableRestraints)
331 {
332  this->RefinableObj::BeginOptimization(allowApproximations,enableRestraints);
333 }
334 
336 {
337  this->CalcPowderPattern();
339 }
340 
341 pair<const CrystVector_REAL*,const RefinableObjClock*>
343 {
344  this->CalcPowderPatternIntegrated();
345  return make_pair(&mPowderPatternIntegratedCalcVariance,
347 }
348 
350 {
351  #ifdef USE_BACKGROUND_MAXLIKE_ERROR
352  return true;
353  #else
354  return false;
355  #endif
356 }
357 
359 {
360 }
361 
363 {
364  VFN_DEBUG_ENTRY("PowderPatternBackground::OptimizeBayesianBackground()",5);
365  TAU_PROFILE("PowderPatternBackground::OptimizeBayesianBackground()","void ()",TAU_DEFAULT);
367  SimplexObj simplex("Simplex Test");
368  simplex.AddRefinableObj(min);
369 
370  long nbcycle;
371  REAL llk=simplex.GetLogLikelihood();
372  long ct=0;
373  cout<<"Initial Chi^2(BayesianBackground)="<<llk<<endl;
375  mBackgroundInterpPointIntensity.max()/1000.0);
376 
377  {
378  char buf [200];
379  sprintf(buf,"Optimizing Background, Cycle %d, Chi^2(Background)=%f",
380  (int)ct,(float)llk);
381  (*fpObjCrystInformUser)((string)buf);
382  }
383  nbcycle=500*mBackgroundNbPoint;
384  simplex.Optimize(nbcycle,false);
385  llk=simplex.GetLogLikelihood();
386  cout<<ct<<", Chi^2(BayesianBackground)="<<llk<<endl;
387 
389  {
390  char buf [200];
391  sprintf(buf,"Done Optimizing Bayesian Background, Chi^2(Background)=%f",(float)llk);
392  (*fpObjCrystInformUser)((string)buf);
393  }
394 
395  LSQNumObj lsq;
396  lsq.SetRefinedObj(min,0,true,true);
397  lsq.PrepareRefParList(true);
399  try {lsq.Refine(10,true,false);}
400  catch(const ObjCrystException &except){};
401 
403  VFN_DEBUG_EXIT("PowderPatternBackground::OptimizeBayesianBackground()",5);
404 }
405 
407 {
408  //Auto-fix points beyond used range
409  unsigned long nbpoint=this->GetParentPowderPattern().GetNbPointUsed();
410  for(long j=0;j<mBackgroundNbPoint;j++)
411  if(this->GetParentPowderPattern().X2Pixel(mBackgroundInterpPointX(j))>nbpoint)
412  {
413  obj.GetPar(&mBackgroundInterpPointIntensity(j)).Print();
414  obj.GetPar(&mBackgroundInterpPointIntensity(j)).SetIsFixed(true);
415  }
416 }
417 
419 {
421 
422  //:TODO: This needs serious optimization !
425  &&(mClockPowderPatternCalc>mInterpolationModel.GetClock())) return;
426  TAU_PROFILE("PowderPatternBackground::CalcPowderPattern()","void ()",TAU_DEFAULT);
427  VFN_DEBUG_MESSAGE("PowderPatternBackground::CalcPowderPattern()",3);
428 
429  const unsigned long nb=mpParentPowderPattern->GetNbPoint();
430  mPowderPatternCalc.resize(nb);
431  if(nb!=0)
432  switch(mInterpolationModel.GetChoice())
433  {
434  case POWDER_BACKGROUND_LINEAR:
435  {
436  VFN_DEBUG_MESSAGE("PowderPatternBackground::CalcPowderPattern()..Linear",2)
437  REAL p1,p2;
438  REAL b1,b2;
439  if(mBackgroundNbPoint==0)
440  {
442  break;
443  }
444  VFN_DEBUG_MESSAGE("PowderPatternBackground::CalcPowderPattern()"<<nb,2)
445  this->InitSpline();
446  REAL *b=mPowderPatternCalc.data();
451  long point=1;
452  for(unsigned long i=0;i<nb;i++)
453  {
454  if(i >= p2)
455  {
456  if(point < mBackgroundNbPoint-1)
457  {
458  b1=b2;
459  p1=p2;
462  point++ ;
463  }
464  }
465  *b = (b1*(p2-i)+b2*(i-p1))/(p2-p1) ;
466  b++;
467  }
468  break;
469  }
470  case POWDER_BACKGROUND_CUBIC_SPLINE:
471  {
472  if(mBackgroundNbPoint==0) mPowderPatternCalc=0;
473  else
474  {
475  this->InitSpline();
476  mPowderPatternCalc=mvSpline((REAL)0,(REAL)1,nb);
477  }
478  break;
479  }
480  }
481  VFN_DEBUG_MESSAGE("PowderPatternBackground::CalcPowderPattern()",3);
482  #ifdef USE_BACKGROUND_MAXLIKE_ERROR
483  {
484  mPowderPatternCalcVariance.resize(nb);
485  const REAL step=mModelVariance*mModelVariance/(REAL)nbPoint;
486  REAL var=0;
487  REAL *p=mPowderPatternCalcVariance.data();
488  for(long i=0;i<nb;i++) {*p++ = var;var +=step;}
489  }
491  #endif
493  VFN_DEBUG_MESSAGE("PowderPatternBackground::CalcPowderPattern():End",3);
494 }
495 
496 void PowderPatternBackground::CalcPowderPattern_FullDeriv(std::set<RefinablePar*> &vPar)
497 {
498  TAU_PROFILE("PowderPatternBackground::CalcPowderPattern_FullDeriv()","void ()",TAU_DEFAULT);
499  const unsigned long nb=mpParentPowderPattern->GetNbPoint();
500  mPowderPattern_FullDeriv.clear();
501  if((nb==0)||(mBackgroundNbPoint==0)) return;
502  for(std::set<RefinablePar*>::iterator par=vPar.begin();par!=vPar.end();++par)
503  {
504  if((*par)==0) mPowderPattern_FullDeriv[*par]=this->GetPowderPatternCalc();
505  else
506  for(int j = 0; j < mBackgroundNbPoint; j++)
507  {
508  if((*par)->GetPointer()!=mBackgroundInterpPointIntensity.data()+j) continue;
509  const REAL step=(*par)->GetDerivStep();
510  (*par)->Mutate(step);
511  mPowderPattern_FullDeriv[*par] =this->GetPowderPatternCalc();
512  (*par)->Mutate(-2*step);
513  mPowderPattern_FullDeriv[*par]-=this->GetPowderPatternCalc();
514  (*par)->Mutate(step);
515  mPowderPattern_FullDeriv[*par]/= step*2;
516  if(MaxAbs(mPowderPattern_FullDeriv[*par])==0)
517  mPowderPattern_FullDeriv[*par].resize(0);
518  }
519  }
520 }
521 
522 void PowderPatternBackground::CalcPowderPatternIntegrated() const
523 {
525 
526  this->CalcPowderPattern();// :TODO: Optimize
529  return;
530 
531  VFN_DEBUG_ENTRY("PowderPatternBackground::CalcPowderPatternIntegrated()",3)
532  TAU_PROFILE("PowderPatternBackground::CalcPowderPatternIntegrated()","void ()",TAU_DEFAULT);
533  const CrystVector_long *pMin=&(mpParentPowderPattern->GetIntegratedProfileMin());
534  const CrystVector_long *pMax=&(mpParentPowderPattern->GetIntegratedProfileMax());
535 
536  const long numInterval=pMin->numElements();
537  mPowderPatternIntegratedCalc.resize(numInterval);
538  REAL * RESTRICT p2=mPowderPatternIntegratedCalc.data();
539  for(int j=0;j<numInterval;j++)
540  {
541  const long max=(*pMax)(j);
542  const REAL * RESTRICT p1=mPowderPatternCalc.data()+(*pMin)(j);
543  *p2=0;
544  for(int k=(*pMin)(j);k<=max;k++) *p2 += *p1++;
545  p2++;
546  }
547  #ifdef USE_BACKGROUND_MAXLIKE_ERROR
548  mPowderPatternIntegratedCalcVariance.resize(numInterval);
550  for(int j=0;j<numInterval;j++)
551  {
552  const long max=(*pMax)(j);
553  const REAL *p1=mPowderPatternCalcVariance.data()+(*pMin)(j);
554  *p2=0;
555  for(int k=(*pMin)(j);k<=max;k++) *p2 += *p1++;
556  p2++;
557  }
559  #endif
561  VFN_DEBUG_EXIT("PowderPatternBackground::CalcPowderPatternIntegrated()",3)
562 }
563 
565 {
566  TAU_PROFILE("PowderPatternBackground::CalcPowderPatternIntegrated_FullDeriv()","void ()",TAU_DEFAULT);
567  //cout<<"PowderPatternBackground::CalcPowderPatternIntegrated_FullDeriv"<<endl;
568  const unsigned long nb=mpParentPowderPattern->GetNbPoint();
569  mPowderPatternIntegrated_FullDeriv.clear();
570  if((nb==0)||(mBackgroundNbPoint==0)) return;
571  for(std::set<RefinablePar*>::iterator par=vPar.begin();par!=vPar.end();++par)
572  {
573  if((*par)==0) mPowderPattern_FullDeriv[*par]=this->GetPowderPatternCalc();
574  else
575  for(int j = 0; j < mBackgroundNbPoint; j++)
576  {
577  if((*par)->GetPointer()!=mBackgroundInterpPointIntensity.data()+j) continue;
578  const REAL step=(*par)->GetDerivStep();
579  (*par)->Mutate(step);
580  mPowderPatternIntegrated_FullDeriv[*par] =*(this->GetPowderPatternIntegratedCalc().first);
581  (*par)->Mutate(-2*step);
582  mPowderPatternIntegrated_FullDeriv[*par]-=*(this->GetPowderPatternIntegratedCalc().first);
583  (*par)->Mutate(step);
584  mPowderPatternIntegrated_FullDeriv[*par]/= step*2;
585  if(MaxAbs(mPowderPatternIntegrated_FullDeriv[*par])==0)
586  mPowderPatternIntegrated_FullDeriv[*par].resize(0);
587  }
588  }
589  #if 0
590  std::map<RefinablePar*, CrystVector_REAL> newDeriv=mPowderPatternIntegrated_FullDeriv;
592  std::vector<const CrystVector_REAL*> v;
593  int n=0;
594  for(std::map<RefinablePar*, CrystVector_REAL>::reverse_iterator pos=mPowderPatternIntegrated_FullDeriv.rbegin();pos!=mPowderPatternIntegrated_FullDeriv.rend();++pos)
595  {
596  cout<<pos->first->GetName()<<":"<<pos->second.size()<<","<<newDeriv[pos->first].size()<<endl;
597  if(pos->second.size()==0) continue;
598  v.push_back(&(pos->second));
599  v.push_back(&(newDeriv[pos->first]));
600  if(++n>8) break;
601  }
602  if(v.size()>0) cout<<"PowderPatternBackground::CalcPowderPatternIntegrated_FullDeriv():"<<endl<<FormatVertVector<REAL>(v,12,1,20)<<endl;
603  //exit(0);
604  #endif
605 }
606 
608 {
609 }
610 const CrystVector_long& PowderPatternBackground::GetBraggLimits()const
611 {
612  // no integration interval for the background
613  mIntegratedReflLimits.resize(0);
614  return mIntegratedReflLimits;
615 }
616 
618 {
621 }
622 
623 void PowderPatternBackground::InitRefParList()
624 {
625  this->ResetParList();
626  REAL *p=mBackgroundInterpPointIntensity.data();
627  char buf [10];
628  string str="Background_Point_";
629  //for(int i=0;i<3;i++)
630  for(int i=0;i<mBackgroundNbPoint;i++)
631  {
632  sprintf(buf,"%d",i);
633  RefinablePar tmp(str+(string)buf,p++,
634  0.,1000.,gpRefParTypeScattDataBackground,REFPAR_DERIV_STEP_RELATIVE,
635  false,true,true,false,1.);
636  tmp.SetGlobalOptimStep(10.);
637  tmp.AssignClock(mClockBackgroundPoint);
638  tmp.SetDerivStep(1e-3);
639  this->AddPar(tmp);
640  }
641  #ifdef USE_BACKGROUND_MAXLIKE_ERROR
642  {
643  RefinablePar tmp("ML Model Error",&mModelVariance,
644  0.,100000.,gpRefParTypeObjCryst,REFPAR_DERIV_STEP_RELATIVE,
645  true,true,true,false,1.);
646  tmp.AssignClock(mClockBackgroundPoint);
647  tmp.SetDerivStep(1e-3);
648  //tmp.SetGlobalOptimStep(10.);
649  tmp.SetGlobalOptimStep(sqrt(mBackgroundInterpPointIntensity.sum()
650  /mBackgroundInterpPointIntensity.numElements()));
651  this->AddPar(tmp);
652  }
653  #endif
654 }
655 
656 void PowderPatternBackground::InitOptions()
657 {
658  VFN_DEBUG_MESSAGE("PowderPatternBackground::InitOptions()",5)
659  static string InterpolationModelName;
660  static string InterpolationModelChoices[2];
661 
662  static bool needInitNames=true;
663  if(true==needInitNames)
664  {
665  InterpolationModelName="Interpolation Model";
666  InterpolationModelChoices[0]="Linear";
667  InterpolationModelChoices[1]="Spline";
668  //InterpolationModelChoices[2]="Chebyshev";
669 
670  needInitNames=false;//Only once for the class
671  }
672  mInterpolationModel.Init(2,&InterpolationModelName,InterpolationModelChoices);
675  mInterpolationModel.SetChoice(1);
676 }
677 
678 void PowderPatternBackground::InitSpline()const
679 {
683 
684  mvSplinePixel.resize(mBackgroundNbPoint);
685 
686  // The points must be in ascending order
687  // Take care later of neutron TOF, as the powder apttern data may not have been initialized yet.
688  CrystVector<long> subs;
689  subs=SortSubs(mBackgroundInterpPointX);
690 
691  if(this->GetParentPowderPattern().GetRadiation().GetWavelengthType()==WAVELENGTH_TOF)
692  {
693  mPointOrder.resize(mBackgroundNbPoint);
694  for(long i=0;i<mBackgroundNbPoint;++i)
695  mPointOrder(i)=subs(mBackgroundNbPoint-1-i);
696  }
697  else mPointOrder=subs;
698 
699  CrystVector_REAL ipixel(mBackgroundNbPoint);
700 
701  for(long i=0;i<mBackgroundNbPoint;++i)
702  {
703  mvSplinePixel(i)=
706  }
707 
708  mvSpline.Init(mvSplinePixel,ipixel);
710 }
711 
712 #ifdef __WX__CRYST__
713 WXCrystObjBasic* PowderPatternBackground::WXCreate(wxWindow* parent)
714 {
715  //:TODO: Check mpWXCrystObj==0
716  mpWXCrystObj=new WXPowderPatternBackground(parent,this);
717  return mpWXCrystObj;
718 }
719 #endif
720 //
722 // PowderPatternDiffraction
723 //
725 PowderPatternDiffraction::PowderPatternDiffraction():
726 mpReflectionProfile(0),
727 mCorrLorentz(*this),mCorrPolar(*this),mCorrSlitAperture(*this),
728 mCorrTextureMarchDollase(*this),mCorrTextureEllipsoid(*this),mCorrTOF(*this),mExtractionMode(false),
729 mpLeBailData(0),mFrozenLatticePar(6),mFreezeLatticePar(false),mFrozenBMatrix(3,3)
730 {
731  VFN_DEBUG_MESSAGE("PowderPatternDiffraction::PowderPatternDiffraction()",10)
732  mIsScalable=true;
733  this->InitOptions();
734  this->SetProfile(new ReflectionProfilePseudoVoigt);
735  this->SetIsIgnoringImagScattFact(true);
736  this->AddSubRefObj(mCorrTextureMarchDollase);
737  this->AddSubRefObj(mCorrTextureEllipsoid);
738  mClockMaster.AddChild(mClockProfilePar);
739  mClockMaster.AddChild(mClockLorentzPolarSlitCorrPar);
740  mClockMaster.AddChild(mpReflectionProfile->GetClockMaster());
741  for(unsigned int i=0;i<3;++i) mFrozenLatticePar(i)=5;
742  for(unsigned int i=3;i<6;++i) mFrozenLatticePar(i)=M_PI/2;
743 }
744 
745 PowderPatternDiffraction::PowderPatternDiffraction(const PowderPatternDiffraction &old):
746 mpReflectionProfile(0),
747 mCorrLorentz(*this),mCorrPolar(*this),mCorrSlitAperture(*this),
748 mCorrTextureMarchDollase(*this),mCorrTextureEllipsoid(*this),mCorrTOF(*this),mExtractionMode(false),
749 mpLeBailData(0),mFrozenLatticePar(6),mFreezeLatticePar(old.FreezeLatticePar()),mFrozenBMatrix(3,3)
750 {
751  this->AddSubRefObj(mCorrTextureMarchDollase);
752  this->AddSubRefObj(mCorrTextureEllipsoid);
753  this->SetIsIgnoringImagScattFact(true);
754  this->SetProfile(old.mpReflectionProfile->CreateCopy());
755  #if 0 //:TODO:
756  if(old.mpLeBailData!=0)
757  {
758  mpLeBailData=new DiffractionDataSingleCrystal(false);
759  *mpLeBailData = *(old.mpLeBailData);
760  }
761  #endif
762  mClockMaster.AddChild(mClockProfilePar);
763  mClockMaster.AddChild(mClockLorentzPolarSlitCorrPar);
764  mClockMaster.AddChild(mpReflectionProfile->GetClockMaster());
765  for(unsigned int i=0;i<6;++i) mFrozenLatticePar(i)=old.GetFrozenLatticePar(i);
766  mFrozenBMatrix=old.GetBMatrix();
767 }
768 
769 PowderPatternDiffraction::~PowderPatternDiffraction()
770 {
771  if(mpReflectionProfile!=0)
772  {
773  this->RemoveSubRefObj(*mpReflectionProfile);
774  delete mpReflectionProfile;
775  }
776 }
778 {
779  const static string className="PowderPatternDiffraction";
780  return className;
781 }
782 
784 {
785  return new PowderPatternDiffraction(*this);
786 }
787 
789 {
790  if(mpParentPowderPattern!=0)
797 }
798 
800 {
801  this->CalcPowderPattern();
802  return mPowderPatternCalc;
803 }
804 
805 pair<const CrystVector_REAL*,const RefinableObjClock*>
807 {
808  this->CalcPowderPatternIntegrated();
810 }
811 
813  const REAL w, const REAL u, const REAL v,
814  const REAL eta0, const REAL eta1)
815 {
816  VFN_DEBUG_MESSAGE("PowderPatternDiffraction::SetReflectionProfilePar()",5)
817  ReflectionProfilePseudoVoigt* p=new ReflectionProfilePseudoVoigt();
818  p->SetProfilePar(w,u,v,eta0,eta1);
819  this->SetProfile(p);
820 }
821 
823 {
824  if(p==mpReflectionProfile) return;
825  if(mpReflectionProfile!=0)
826  {
827  this->RemoveSubRefObj(*mpReflectionProfile);
828  delete mpReflectionProfile;
829  }
830  mpReflectionProfile= p;
831  this->AddSubRefObj(*mpReflectionProfile);
832  mClockMaster.AddChild(mpReflectionProfile->GetClockMaster());
833 }
834 
836 {
837  return *mpReflectionProfile;
838 }
839 
841 {
842  return *mpReflectionProfile;
843 }
844 
845 // Disable the base-class function.
846 void PowderPatternDiffraction::GenHKLFullSpace(
847  const REAL maxTheta, const bool useMultiplicity)
848 {
849  // This should be never called.
850  abort();
851 }
852 
853 void PowderPatternDiffraction::GenHKLFullSpace()
854 {
855  VFN_DEBUG_ENTRY("PowderPatternDiffraction::GenHKLFullSpace():",5)
856  float stol;
857  if(this->GetRadiation().GetWavelengthType()==WAVELENGTH_TOF)
858  stol=mpParentPowderPattern->X2STOL(mpParentPowderPattern->GetPowderPatternXMin());
859  else
860  stol=mpParentPowderPattern->X2STOL(mpParentPowderPattern->GetPowderPatternXMax());
861  if(stol>1) stol=1; // Do not go beyond 0.5 A resolution (mostly for TOF data)
862  this->ScatteringData::GenHKLFullSpace2(stol,true);
863  if((mExtractionMode) && (mFhklObsSq.numElements()!=this->GetNbRefl()))
864  {// Reflections changed, so ScatteringData::PrepareHKLarrays() probably reseted mFhklObsSq
865  VFN_DEBUG_ENTRY("PowderPatternDiffraction::GenHKLFullSpace(): need to reset observed intensities",7)
866  mFhklObsSq.resize(this->GetNbRefl());
867  mFhklObsSq=100;
868  }
869  mCorrTextureEllipsoid.InitRefParList();// #TODO: SHould this be here ?
870  VFN_DEBUG_EXIT("PowderPatternDiffraction::GenHKLFullSpace():"<<this->GetNbRefl(),5)
871 }
872 void PowderPatternDiffraction::BeginOptimization(const bool allowApproximations,
873  const bool enableRestraints)
874 {
875  if(mUseFastLessPreciseFunc!=allowApproximations)
876  {
877  mClockProfileCalc.Reset();
878  mClockGeomStructFact.Reset();
879  mClockStructFactor.Reset();
881  }
882  mUseFastLessPreciseFunc=allowApproximations;
883  this->GetNbReflBelowMaxSinThetaOvLambda();
884  this->RefinableObj::BeginOptimization(allowApproximations,enableRestraints);
885 }
887 {
888  if(mOptimizationDepth==1)
889  {
890  if(mUseFastLessPreciseFunc==true)
891  {
892  mClockProfileCalc.Reset();
893  mClockGeomStructFact.Reset();
894  mClockStructFactor.Reset();
896  }
897  mUseFastLessPreciseFunc=false;
898  this->GetNbReflBelowMaxSinThetaOvLambda();
899  }
901 }
902 
904 {
905  if(mUseFastLessPreciseFunc!=allow)
906  {
907  mClockProfileCalc.Reset();
908  mClockGeomStructFact.Reset();
909  mClockStructFactor.Reset();
911  }
912  mUseFastLessPreciseFunc=allow;
913  this->GetNbReflBelowMaxSinThetaOvLambda();
915 }
916 
918  CrystVector_uint & groupIndex,
919  unsigned int &first) const
920 {
921  // One group for all profile parameters
922  unsigned int index=0;
923  VFN_DEBUG_MESSAGE("PowderPatternDiffraction::GetGeneGroup()",4)
924  for(long i=0;i<obj.GetNbPar();i++)
925  for(long j=0;j<this->GetNbPar();j++)
926  if(&(obj.GetPar(i)) == &(this->GetPar(j)))
927  {
928  //if(this->GetPar(j).GetType()->IsDescendantFromOrSameAs())
929  //{
930  if(index==0) index=first++;
931  groupIndex(i)=index;
932  //}
933  //else //no parameters other than unit cell
934  }
935 }
936 
938 {
939 
941 }
942 
943 pair<const CrystVector_REAL*,const RefinableObjClock*>
945 {
946  this->CalcPowderPatternIntegrated();
947  return make_pair(&mPowderPatternIntegratedCalcVariance,
949 }
950 
952 {
953  return true;
954 }
955 
957 {
958  bool reprep=(mpCrystal!=0);
959  this->ScatteringData::SetCrystal(crystal);
960  // Check if we use DE-PV
961  if(mpReflectionProfile!=0)
962  if(mpReflectionProfile->GetClassName()=="ReflectionProfileDoubleExponentialPseudoVoigt")
963  {
965  =dynamic_cast<ReflectionProfileDoubleExponentialPseudoVoigt*>(mpReflectionProfile);
966  p->SetUnitCell((UnitCell)crystal);
967  }
968  mClockHKL.Reset();
969  if(reprep) this->Prepare();
970 }
971 
974 
975 void PowderPatternDiffraction::SetExtractionMode(const bool extract,const bool init)
976 {
977  VFN_DEBUG_ENTRY("PowderPatternDiffraction::SetExtractionMode(),ExtractionMode="<<mExtractionMode<<", nbrefl="<<this->GetNbRefl(),7)
978  mExtractionMode=extract;
979  bool needInit=false;
980  if(extract)
981  {
982  this->FreezeLatticePar(false);
983  this->Prepare();
984  mFhklObsSq.resizeAndPreserve(this->GetNbRefl());
985  }
986  if(extract && (!init) && (mpLeBailData!=0))
987  {
988  // Re-use existing Le Bail data, if list of hkl's is consistent
989  const long nbrefl=this->GetNbReflBelowMaxSinThetaOvLambda();
990  if(nbrefl==mpLeBailData->GetNbRefl())
991  {
992  for(int i=0;i<nbrefl;++i)
993  {
994  if( (mpLeBailData->GetH()(i)==this->GetH()(i))
995  &&(mpLeBailData->GetK()(i)==this->GetK()(i))
996  &&(mpLeBailData->GetL()(i)==this->GetL()(i)))
997  {
998  mFhklObsSq(i)=mpLeBailData->GetFhklObsSq()(i);
999  }
1000  else
1001  {
1002  needInit=true;
1003  VFN_DEBUG_MESSAGE("PowderPatternDiffraction::SetExtractionMode():: Forcing initialize, cannot re-use: hkl list differs",10);
1004  break;
1005  }
1006 
1007  }
1008  }
1009  else
1010  {
1011  needInit=true;
1012  VFN_DEBUG_MESSAGE("PowderPatternDiffraction::SetExtractionMode():: Forcing initialize, cannot re-use: different number of reflections",10);
1013  }
1014  }
1015  if((extract && init) || needInit)
1016  {
1017  VFN_DEBUG_MESSAGE("PowderPatternDiffraction::SetExtractionMode():: Initializing intensities to 100",10);
1018  mFhklObsSq=100;
1019  }
1020  if((mExtractionMode==false)&&(mFhklObsSq.numElements()>0))
1021  {// Leaving extraction mode, so update extracted single crystal data
1022  VFN_DEBUG_ENTRY("PowderPatternDiffraction::SetExtractionMode(),LEAVING Le Bail Mode",7)
1023  if(mpLeBailData==0) mpLeBailData=new DiffractionDataSingleCrystal(this->GetCrystal(),false);
1024  // Update wavelength & name
1025  mpLeBailData->SetWavelength(this->GetRadiation().GetWavelength()(0));
1026  mpLeBailData->SetRadiationType(this->GetRadiation().GetRadiationType());
1027  char buf[200];
1028  sprintf(buf,"LeBail (d=%4.2fA):",1/(2*abs(mMaxSinThetaOvLambda)+1e-6));
1029  mpLeBailData->SetName(string(buf)+this->GetCrystal().GetName());
1030 
1031  const unsigned long nbrefl=this->GetNbReflBelowMaxSinThetaOvLambda();
1032  CrystVector_REAL iobs(nbrefl),sigma(nbrefl);
1033  CrystVector_long h(nbrefl),k(nbrefl),l(nbrefl);
1034  sigma=1;
1035  for(unsigned long i=0;i<nbrefl;++i)
1036  {
1037  h(i)=mIntH(i);
1038  k(i)=mIntK(i);
1039  l(i)=mIntL(i);
1040  iobs(i)=mFhklObsSq(i);
1041  }
1042  mpLeBailData->SetHklIobs(h,k,l,iobs,sigma);
1043  // Erase mFhklObsSq - only used during extraction mode.
1044  mFhklObsSq.resize(0);
1045  }
1046  mClockIhklCalc.Reset();mClockMaster.Reset();
1047  mClockFhklObsSq.Click();
1048  VFN_DEBUG_EXIT("PowderPatternDiffraction::SetExtractionMode(),ExtractionMode="<<mExtractionMode<<", nbrefl="<<this->GetNbRefl(),7)
1049 }
1050 
1051 bool PowderPatternDiffraction::GetExtractionMode()const{return mExtractionMode;}
1052 
1054 {
1055  VFN_DEBUG_ENTRY("PowderPatternDiffraction::ExtractLeBail()",7)
1056  TAU_PROFILE("PowderPatternDiffraction::ExtractLeBail()","void (int)",TAU_DEFAULT);
1057 
1058  if(mExtractionMode==false) this->SetExtractionMode(true,true);// Should not have to do this here !
1059  if(mFhklObsSq.numElements()!=this->GetNbRefl())
1060  {//Something went wrong !
1061  VFN_DEBUG_ENTRY("PowderPatternDiffraction::ExtractLeBail() mFhklObsSq.size() != NbRefl !!!!!!",7)
1062  mFhklObsSq.resize(this->GetNbRefl());
1063  mFhklObsSq=100;
1064  }
1065  // First get the observed powder pattern, minus the contribution of all other phases.
1066  CrystVector_REAL obs,iextract,calc;
1067  iextract=mFhklObsSq;
1068  mFhklObsSq=0;
1069  mClockFhklObsSq.Click();
1070  // Get the observed and calculated powder pattern (excluding this diffraction phase)
1073  mFhklObsSq=iextract;
1074  mClockFhklObsSq.Click();
1075  // We take here the reflections which are centered below the max(sin(theta)/lambda)
1076  // actually more reflections are calculated, but the pattern is only calculated up to
1077  // max(sin(theta)/lambda).
1078  const unsigned long nbrefl=this->ScatteringData::GetNbReflBelowMaxSinThetaOvLambda();
1079  iextract=0;
1080  for(;nbcycle>0;nbcycle--)
1081  {
1082  //cout<<"PowderPatternDiffraction::ExtractLeBail(): cycle #"<<nbcycle<<endl;
1083  calc=this->GetPowderPatternCalc();
1084  for(unsigned int k0=0;k0<nbrefl;++k0)
1085  {
1086  REAL s1=0;
1087  //cout<<mH(k0)<<" "<<mK(k0)<<" "<<mL(k0)<<" , Iobs=??"<<endl;
1088  long last=mvReflProfile[k0].last,first;
1090  if(mvReflProfile[k0].first<0)first=0;
1091  else first=(mvReflProfile[k0].first);
1092  const REAL *p1=mvReflProfile[k0].profile.data()+(first-mvReflProfile[k0].first);
1093  const REAL *p2=calc.data()+first;
1094  const REAL *pobs=obs.data()+first;
1095  for(unsigned int i=first;i<=last;++i)
1096  {
1097  const REAL s2=*p2++;
1098  const REAL tmp=*pobs++ * *p1++;
1099  if( (s2<1e-8) ) // || (tmp<=0)
1100  {// Avoid <0 intensities (should not happen, it means profile is <0)
1101  //cout<<"S2? "<< int(mH(k0))<<" "<<int(mK(k0))<<" "<<int(mL(k0)) <<" calc(i="<<i<<")"<<calc(i)<<" obs(i="<<i<<")="<<obs(i)<<", tmp="<<tmp<<" profile(i)="<<mvReflProfile[k0].profile(i-mvReflProfile[k0].first)<<" "<<mFhklObsSq(k0)<<endl;
1102  continue ;
1103  }
1104  s1 += tmp /s2;
1105  //cout<<" "<<s2<<" "<<obs(i)<<" "<<mvReflProfile[k0].profile(i-mvReflProfile[k0].first)<<" "<<mFhklObsSq(k0)<<endl;
1106  }
1107  if((s1>1e-8)&&(!ISNAN_OR_INF(s1))) iextract(k0)=s1*mFhklObsSq(k0);
1108  else iextract(k0)=1e-8;//:KLUDGE: should <0 intensities be allowed ?
1109  //if(nbcycle==1) cout<<" Le Bail "<<int(mH(k0))<<" "<<int(mK(k0))<<" "<<int(mL(k0))<<" , Iobs="<<iextract(k0)<<endl;
1110  }
1111  mFhklObsSq=iextract;
1112  if(this->GetCrystal().GetScatteringComponentList().GetNbComponent()>0)
1113  {// Change scale factor if we have some atoms in the structure
1114  const REAL* p1=this->GetFhklCalcSq() .data();
1115  const REAL* p2=mFhklObsSq.data();
1116  REAL tmp1=0,tmp2=0;
1117  for(long i=nbrefl;i>0;i--)
1118  {
1119  tmp1 += (*p1) * (*p2++);
1120  tmp2 += (*p1) * (*p1);
1121  p1++;
1122  }
1123  //cout<<"SCALING: tmp2="<<tmp2<<",tmp1="<<tmp1<<endl;
1124  mFhklObsSq*=tmp2/tmp1;
1125  }
1126  mClockFhklObsSq.Click();
1127  //cout<<"PowderPatternDiffraction::ExtractLeBail():results (scale factor="<<mpParentPowderPattern->GetScaleFactor(*this)*1e6<<")" <<endl<< FormatVertVectorHKLFloats<REAL>(mH,mK,mL,this->GetFhklCalcSq(),mFhklObsSq,10,4,nbrefl)<<endl;
1128  }
1129  // Store extracted data in a single crystal data object
1130  if(mpLeBailData==0) mpLeBailData=new DiffractionDataSingleCrystal(*mpCrystal,false);
1131  {
1132  VFN_DEBUG_MESSAGE("PowderPatternDiffraction::ExtractLeBail(): creating single crystal extracted data",7)
1133  CrystVector_REAL iobs(nbrefl),sigma(nbrefl);
1134  CrystVector_long h(nbrefl),k(nbrefl),l(nbrefl);
1135  sigma=1;
1136  for(unsigned long i=0;i<nbrefl;++i)
1137  {
1138  h(i)=mIntH(i);
1139  k(i)=mIntK(i);
1140  l(i)=mIntL(i);
1141  iobs(i)=mFhklObsSq(i);
1142  }
1143  mpLeBailData->SetHklIobs(h,k,l,iobs,sigma);
1144  }
1145  VFN_DEBUG_EXIT("PowderPatternDiffraction::ExtractLeBail()mFhklObsSq.size()=="<<mFhklObsSq.numElements(),7)
1146 }
1148 {
1149  if(this->IsBeingRefined()) return mNbReflUsed;
1150  VFN_DEBUG_MESSAGE("PowderPattern::GetNbReflBelowMaxSinThetaOvLambda()",4)
1151  this->CalcPowderReflProfile();
1152  const long nbpoint=mpParentPowderPattern->GetNbPointUsed();
1153  if((mNbReflUsed>0)&&(mNbReflUsed<mNbRefl))
1154  {
1155  if( (mvReflProfile[mNbReflUsed ].first>nbpoint)
1156  &&(mvReflProfile[mNbReflUsed-1].first<=nbpoint)) return mNbReflUsed;
1157  }
1158 
1159  if((mNbReflUsed==mNbRefl) && (mvReflProfile[mNbReflUsed-1].profile.numElements()>0))
1160  if(mvReflProfile[mNbReflUsed-1].first<=nbpoint)return mNbReflUsed;
1161 
1162 
1163  long i;
1164  for(i=0;i<mNbRefl;i++)
1165  {
1166  if(mvReflProfile[i].first>nbpoint) break;
1167  }
1168  if(i!=mNbReflUsed)
1169  {
1170  mNbReflUsed=i;
1171  mClockNbReflUsed.Click();
1172  VFN_DEBUG_MESSAGE("->Changed Max sin(theta)/lambda="<<mMaxSinThetaOvLambda\
1173  <<" nb refl="<<mNbReflUsed,4)
1174  }
1175  return mNbReflUsed;
1176 }
1177 
1178 void PowderPatternDiffraction::SetFrozenLatticePar(const unsigned int i, REAL v)
1179 {
1180  const REAL old=mFrozenLatticePar(i);
1181  cout<<"PowderPatternDiffraction::SetFrozenLatticePar("<<i<<":"<<v<<")"<<endl;
1182  if(old==v) return;
1183  mFrozenLatticePar(i)=v;
1184  this->CalcFrozenBMatrix();
1185 }
1186 
1187 REAL PowderPatternDiffraction::GetFrozenLatticePar(const unsigned int i) const {return mFrozenLatticePar(i);}
1188 
1190 {
1191  cout<<"PowderPatternDiffraction::FreezeLatticePar("<<use<<")"<<endl;
1192  VFN_DEBUG_MESSAGE("PowderPatternDiffraction::FreezeLatticePar("<<use<<")", 10)
1193  if(use==mFreezeLatticePar) return;
1194  mFreezeLatticePar=use;
1195  mFrozenLatticePar=this->GetCrystal().GetLatticePar();
1196  if(use) this->CalcFrozenBMatrix();
1197  mClockTheta.Reset();
1198  this->UpdateDisplay();
1199 }
1200 
1201 bool PowderPatternDiffraction::FreezeLatticePar() const {return mFreezeLatticePar;}
1202 
1204 {
1205  this->GetNbReflBelowMaxSinThetaOvLambda();
1207  TAU_PROFILE("PowderPatternDiffraction::CalcPowderPattern()-Apply profiles","void (bool)",TAU_DEFAULT);
1208 
1209  VFN_DEBUG_ENTRY("PowderPatternDiffraction::CalcPowderPattern():",3)
1210 
1211  // :TODO: Can't do this as this is non-const
1212  //if(this->GetCrystal().GetSpaceGroup().GetClockSpaceGroup()>mClockHKL)
1213  // this->GenHKLFullSpace();
1214  //
1215  // The workaround is to call Prepare() (non-const) before every calculation
1216  // when a modifictaion may have occured.
1217 
1218  this->CalcIhkl();
1219  this->CalcPowderReflProfile();
1220 
1221  if( (mClockPowderPatternCalc>mClockIhklCalc)
1222  &&(mClockPowderPatternCalc>mClockProfileCalc)) return;
1223 
1224  if(true) //:TODO: false == mUseFastLessPreciseFunc
1225  {
1226  const long nbRefl=this->GetNbRefl();
1227  VFN_DEBUG_MESSAGE("PowderPatternDiffraction::CalcPowderPattern\
1228 Applying profiles for "<<nbRefl<<" reflections",2)
1229  long step; // number of reflections at the same place and with the same (assumed) profile
1230  const long specNbPoints=mpParentPowderPattern->GetNbPoint();
1231  mPowderPatternCalc.resize(specNbPoints);
1233  const bool useML= (mIhklCalcVariance.numElements() != 0);
1234  if(useML)
1235  {
1236  mPowderPatternCalcVariance.resize(specNbPoints);
1238  }
1239  else mPowderPatternCalcVariance.resize(0);
1240  VFN_DEBUG_MESSAGE("PowderPatternDiffraction::CalcPowderPattern() Has variance:"<<useML,2)
1241 
1242  for(long i=0;i<mNbRefl;i += step)
1243  {
1244  if(mvReflProfile[i].profile.numElements()==0)
1245  {
1246  step=1;
1247  if(i>=mNbReflUsed) break;// After sin(theta)/lambda limit
1248  else continue; // before beginning of pattern ?
1249  }
1250 
1251  VFN_DEBUG_MESSAGE("PowderPatternDiffraction::CalcPowderPattern()#"<<i,2)
1252  REAL intensity=0.;
1253  REAL var=0.;
1254  //check if the next reflection is at the same theta. If this is true,
1255  //Then assume that the profile is exactly the same, unless it is anisotropic
1256  for(step=0; ;)
1257  {
1258  intensity += mIhklCalc(i + step);
1259  if(useML) var += mIhklCalcVariance(i + step);
1260  step++;
1261  if(mpReflectionProfile->IsAnisotropic()) break;// Anisotropic profiles
1262  if( (i+step) >= nbRefl) break;
1263  if(mSinThetaLambda(i+step) > (mSinThetaLambda(i)+1e-5) ) break;
1264  }
1265  VFN_DEBUG_MESSAGE("Apply profile(Monochromatic)Refl("<<i<<")"\
1266  <<mIntH(i)<<" "<<mIntK(i)<<" "<<mIntL(i)<<" "\
1267  <<" I="<<intensity<<" stol="<<mSinThetaLambda(i)\
1268  <<",pixel #"<<mvReflProfile[i].first<<"->"<<mvReflProfile[i].last,2)
1269  {
1270  const unsigned long first=mvReflProfile[i].first,last=mvReflProfile[i].last;
1271  const REAL *p2 = mvReflProfile[i].profile.data();
1272  REAL *p3 = mPowderPatternCalc.data()+first;
1273  for(unsigned long j=first;j<=last;j++) *p3++ += *p2++ * intensity;
1274  if(useML)
1275  {
1276  const REAL *p2 = mvReflProfile[i].profile.data();
1277  REAL *p3 = mPowderPatternCalcVariance.data()+first;
1278  for(unsigned long j=first;j<=last;j++) *p3++ += *p2++ * var;
1279  }
1280  }
1281  }
1282  }
1283  else
1284  {
1285  //:TODO:
1286  throw ObjCrystException("PowderPatternDiffraction::CalcPowderPattern() : \
1287  FAST option not yet implemented !");
1288  }
1290  VFN_DEBUG_EXIT("PowderPatternDiffraction::CalcPowderPattern: End.",3)
1291 }
1292 
1293 void PowderPatternDiffraction::CalcPowderPattern_FullDeriv(std::set<RefinablePar*> &vPar)
1294 {
1295  TAU_PROFILE("PowderPatternDiffraction::CalcPowderPattern_FullDeriv()","void ()",TAU_DEFAULT);
1296  //cout<<"PowderPatternDiffraction::CalcPowderPattern_FullDeriv"<<endl;
1297  this->CalcPowderPattern();
1298  bool notYetDerivProfiles=true;
1299  mIhkl_FullDeriv.clear();
1300  mvReflProfile_FullDeriv.clear();
1301  mPowderPattern_FullDeriv.clear();
1302  for(std::set<RefinablePar*>::iterator par=vPar.begin();par!=vPar.end();++par)
1303  {
1304  if(*par==0) continue;
1305  if((*par)->IsFixed()) continue;
1306  if((*par)->IsUsed()==false) continue;
1307  if( (*par)->GetType()->IsDescendantFromOrSameAs(gpRefParTypeScatt)
1308  ||(*par)->GetType()->IsDescendantFromOrSameAs(gpRefParTypeScattPow))
1309  {
1310  this->CalcIhkl_FullDeriv(vPar);
1311  }
1312  if(notYetDerivProfiles)
1313  {
1314  if( (*par)->GetType()->IsDescendantFromOrSameAs(gpRefParTypeRadiation)
1315  ||(*par)->GetType()->IsDescendantFromOrSameAs(gpRefParTypeUnitCell)
1316  ||(*par)->GetType()->IsDescendantFromOrSameAs(gpRefParTypeScattDataCorrPos)
1317  ||(*par)->GetType()->IsDescendantFromOrSameAs(gpRefParTypeScattDataProfile))
1318  {
1319  this->CalcPowderReflProfile_FullDeriv(vPar);
1320  notYetDerivProfiles=false;
1321  }
1322  }
1323  }
1324 
1325  //this->CalcPowderReflProfile();
1326  for(std::set<RefinablePar*>::iterator par=vPar.begin();par!=vPar.end();++par)
1327  {
1328  if(*par==0) mPowderPattern_FullDeriv[*par]=this->GetPowderPatternCalc();
1329  else
1330  {
1331  if((*par)->IsFixed()) continue;
1332  if((*par)->IsUsed()==false) continue;
1333  if(mIhkl_FullDeriv[*par].size()!=0)
1334  {
1335  const long nbRefl=this->GetNbRefl();
1336  long step; // number of reflections at the same place and with the same (assumed) profile
1337  const long specNbPoints=mpParentPowderPattern->GetNbPoint();
1338  mPowderPattern_FullDeriv[*par].resize(specNbPoints);
1339  mPowderPattern_FullDeriv[*par]=0;
1340 
1341  for(long i=0;i<mNbReflUsed;i += step)
1342  {
1343  if(mvReflProfile[i].profile.numElements()==0)
1344  {
1345  step=1;
1346  if(i>=mNbReflUsed) break;
1347  else continue;
1348  }
1349 
1350  REAL intensity=0.;
1351  //check if the next reflection is at the same theta. If this is true,
1352  //Then assume that the profile is exactly the same, unless it is anisotropic
1353  for(step=0; ;)
1354  {
1355  intensity += mIhkl_FullDeriv[*par](i + step);
1356  step++;
1357  if(mpReflectionProfile->IsAnisotropic()) break;// Anisotropic profiles
1358  if( (i+step) >= nbRefl) break;
1359  if(mSinThetaLambda(i+step) > (mSinThetaLambda(i)+1e-5) ) break;
1360  }
1361  {
1362  const unsigned long first=mvReflProfile[i].first,last=mvReflProfile[i].last;
1363  const REAL *p2 = mvReflProfile[i].profile.data();
1364  REAL *p3 = mPowderPattern_FullDeriv[*par].data()+first;
1365  for(unsigned long j=first;j<=last;j++) *p3++ += *p2++ * intensity;
1366  }
1367  }
1368  }
1369  if(mvReflProfile_FullDeriv[*par].size()!=0)
1370  {
1371  const long nbRefl=this->GetNbRefl();
1372  long step; // number of reflections at the same place and with the same (assumed) profile
1373  const long specNbPoints=mpParentPowderPattern->GetNbPoint();
1374  mPowderPattern_FullDeriv[*par].resize(specNbPoints);
1375  mPowderPattern_FullDeriv[*par]=0;// :TODO: use only the number of points actually used
1376  cout<<__FILE__<<":"<<__LINE__<<":PowderPatternDiffraction::CalcPowderPattern_FullDeriv():par="<<(*par)->GetName()<<endl;
1377  for(long i=0;i<mNbReflUsed;i += step)
1378  {
1379  if(mvReflProfile[i].profile.numElements()==0)
1380  {
1381  step=1;
1382  if(i>=mNbReflUsed) break;
1383  else continue;
1384  }
1385 
1386  REAL intensity=0.;
1387  //check if the next reflection is at the same theta. If this is true,
1388  //Then assume that the profile is exactly the same, unless it is anisotropic
1389  for(step=0; ;)
1390  {
1391  intensity += mIhklCalc(i + step);
1392  step++;
1393  if(mpReflectionProfile->IsAnisotropic()) break;// Anisotropic profiles
1394  if( (i+step) >= nbRefl) break;
1395  if(mSinThetaLambda(i+step) > (mSinThetaLambda(i)+1e-5) ) break;
1396  }
1397  if(mvReflProfile_FullDeriv[*par][i].size()>0)// Some profiles may be unaffected by a given parameter
1398  {
1399  const unsigned long first=mvReflProfile[i].first,last=mvReflProfile[i].last;
1400  const REAL *p2 = mvReflProfile_FullDeriv[*par][i].data();
1401  REAL *p3 = mPowderPattern_FullDeriv[*par].data()+first;
1402  for(unsigned long j=first;j<=last;j++) *p3++ += *p2++ * intensity;
1403  }
1404  }
1405  }
1406  }
1407  }
1408  #if 0
1409  std::map<RefinablePar*, CrystVector_REAL> newDeriv=mPowderPattern_FullDeriv;
1410  this->PowderPatternComponent::CalcPowderPattern_FullDeriv(vPar);
1411  std::vector<const CrystVector_REAL*> v;
1412  int n=0;
1413  for(std::map<RefinablePar*, CrystVector_REAL>::reverse_iterator pos=mPowderPattern_FullDeriv.rbegin();pos!=mPowderPattern_FullDeriv.rend();++pos)
1414  {
1415  if(pos->first==0) continue;
1416  if(pos->second.size()==0) continue;
1417  v.push_back(&(newDeriv[pos->first]));
1418  v.push_back(&(pos->second));
1419  cout<<pos->first->GetName()<<":"<<pos->second.size()<<","<<newDeriv[pos->first].size()<<endl;
1420  if(++n>8) break;
1421  }
1422  cout<<FormatVertVector<REAL>(v,16,4,1000)<<endl;
1423  //exit(0);
1424  #endif
1425 }
1426 
1427 void PowderPatternDiffraction::CalcPowderPatternIntegrated() const
1428 {
1429  this->GetNbReflBelowMaxSinThetaOvLambda();
1431  TAU_PROFILE("PowderPatternDiffraction::CalcPowderPatternIntegrated()","void (bool)",TAU_DEFAULT);
1432  TAU_PROFILE_TIMER(timer1,"PowderPatternDiffraction::CalcPowderPatternIntegrated()1","", TAU_FIELD);
1433  TAU_PROFILE_TIMER(timer2,"PowderPatternDiffraction::CalcPowderPatternIntegrated()2","", TAU_FIELD);
1434 
1435  this->CalcIhkl();
1436  TAU_PROFILE_START(timer1);
1437  this->PrepareIntegratedProfile();
1438  TAU_PROFILE_STOP(timer1);
1439 
1440  if( (mClockPowderPatternIntegratedCalc>mClockIhklCalc)
1441  &&(mClockPowderPatternIntegratedCalc>mClockIntegratedProfileFactor)
1443  return;
1444  VFN_DEBUG_ENTRY("PowderPatternDiffraction::CalcPowderPatternIntegrated()",3)
1445  const long nbRefl=this->GetNbRefl();
1446 
1447  const long nb=mpParentPowderPattern->GetIntegratedProfileMin().numElements();
1448  mPowderPatternIntegratedCalc.resize(nb);
1450  const bool useML= (mIhklCalcVariance.numElements() != 0);
1451  if(useML)
1452  {
1455  }
1456  else mPowderPatternIntegratedCalcVariance.resize(0);
1457  const REAL * RESTRICT psith=mSinThetaLambda.data();
1458  const REAL * RESTRICT pI=mIhklCalc.data();
1459  const REAL * RESTRICT pIvar=mIhklCalcVariance.data();
1460  vector< pair<unsigned long, CrystVector_REAL> >::const_iterator pos;
1461  pos=mIntegratedProfileFactor.begin();
1462  TAU_PROFILE_START(timer2);
1463  for(long i=0;i<mNbReflUsed;)
1464  {
1465  VFN_DEBUG_MESSAGE("PowderPatternDiffraction::CalcPowderPatternIntegrated():"<<i,2)
1466  REAL intensity=0.;
1467  REAL var=0.;
1468  const REAL thmax=*psith+1e-5;
1469  //check if the next reflection is at the same theta. If this is true,
1470  //Then assume that the profile is exactly the same, unless profiles are anisotropic.
1471  for(;;)
1472  {
1473  intensity += *pI++;
1474  if(useML) var += *pIvar++;
1475  if( ++i >= nbRefl) break;
1476  if( *(++psith) > thmax ) break;
1477  if(mpReflectionProfile->IsAnisotropic()) break;// Anisotropic profile
1478  ++pos;
1479  }
1480  VFN_DEBUG_MESSAGE("PowderPatternDiffraction::CalcPowderPatternIntegrated():"<<i,2)
1481  REAL * RESTRICT pData=mPowderPatternIntegratedCalc.data()+pos->first;
1482  const REAL * RESTRICT pFact=pos->second.data();
1483  const unsigned long nb=pos->second.numElements();
1484  //cout <<i<<" - "<< intensity<<"*:";
1485  for(unsigned long j=nb;j>0;j--)
1486  {
1487  //cout <<pos->first+j<<"("<<*pFact<<","<<*pData<<") ";
1488  *pData++ += intensity * *pFact++ ;
1489  }
1490  //cout<<endl;
1491 
1492  if(useML)
1493  {
1494  const REAL * RESTRICT pFact=pos->second.data();
1495  REAL * RESTRICT pVar=mPowderPatternIntegratedCalcVariance.data()+pos->first;
1496  for(unsigned long j=nb;j>0;j--) *pVar++ += var * *pFact++ ;
1497  }
1498  ++pos;
1499  }
1500  TAU_PROFILE_STOP(timer2);
1501  #ifdef __DEBUG__
1502  if(gVFNDebugMessageLevel<3)
1503  {
1504  this->CalcPowderPattern();
1505  CrystVector_REAL integr(nb),min(nb),max(nb),diff(nb),index(nb);
1506  integr=0;
1507  for(long i=0;i<nb;i++)
1508  {
1509  index(i)=i;
1512  integr(i)=0;
1515  {
1516  integr(i) += mPowderPatternCalc(j);
1517  }
1518  diff(i)=1.-mPowderPatternIntegratedCalc(i)/integr(i);
1519  }
1520  cout << "Integrated intensities, Component"<<endl
1522  (index,min,max,integr,mPowderPatternIntegratedCalc,diff,20,6)<<endl;
1523  }
1524  #endif
1526  VFN_DEBUG_EXIT("PowderPatternDiffraction::CalcPowderPatternIntegrated",3)
1527 }
1528 
1530 {
1531  TAU_PROFILE("PowderPatternDiffraction::CalcPowderPatternIntegrated_FullDeriv()","void ()",TAU_DEFAULT);
1532  //cout<<"PowderPatternDiffraction::CalcPowderPatternIntegrated_FullDeriv"<<endl;
1533  //this->PowderPatternComponent::CalcPowderPatternIntegrated_FullDeriv(vPar);
1534  //return;
1535 
1536  this->CalcPowderPatternIntegrated();
1537  this->CalcIhkl_FullDeriv(vPar);
1538  const long nbRefl=this->GetNbRefl();
1539  //#define PowderPatternDiffraction_CalcPowderPatternIntegrated_FullDerivDEBUG
1540  #ifdef PowderPatternDiffraction_CalcPowderPatternIntegrated_FullDerivDEBUG
1542  std::map<RefinablePar*, CrystVector_REAL> oldDeriv=mPowderPatternIntegrated_FullDeriv;
1543  #endif
1544  const long nbprof=mpParentPowderPattern->GetIntegratedProfileMin().size();
1545  long ctpar=0;
1546  mPowderPatternIntegrated_FullDeriv.clear();
1547  for(std::set<RefinablePar*>::iterator par=vPar.begin();par!=vPar.end();++par)
1548  {
1549  if(*par==0) mPowderPatternIntegrated_FullDeriv[*par]=mPowderPatternIntegratedCalc;
1550  else
1551  {
1552  if(mIhkl_FullDeriv[*par].size()==0) continue;
1553  if(mPowderPatternIntegrated_FullDeriv[*par].size()==0)
1554  {
1555  mPowderPatternIntegrated_FullDeriv[*par].resize(nbprof);
1556  mPowderPatternIntegrated_FullDeriv[*par]=0;
1557  }
1558  const REAL * RESTRICT psith=mSinThetaLambda.data();
1559  const REAL * RESTRICT pI=mIhkl_FullDeriv[*par].data();
1560  // :TODO: also handle derivatives of profile parameters ? Though one should not
1561  // refine profiles using integrated profiles !
1562  vector< pair<unsigned long, CrystVector_REAL> >::const_iterator pos=mIntegratedProfileFactor.begin();
1563 
1564  for(long i=0;i<mNbReflUsed;)
1565  {
1566  REAL intensity=0.;
1567  const REAL thmax=*psith+1e-5;
1568  for(;;)
1569  {
1570  intensity += *pI++;
1571  if( ++i >= nbRefl) break;
1572  if( *(++psith) > thmax ) break;
1573  if(mpReflectionProfile->IsAnisotropic()) break;// Anisotropic profile
1574  ++pos;
1575  }
1576  REAL * RESTRICT pData=mPowderPatternIntegrated_FullDeriv[*par].data()+pos->first;
1577  const REAL * RESTRICT pFact=pos->second.data();
1578  const unsigned long nb=pos->second.numElements();
1579  #ifdef PowderPatternDiffraction_CalcPowderPatternIntegrated_FullDerivDEBUG
1580  if((i<5)&&(ctpar<8)) cout<<__FILE__<<":"<<__LINE__<<":"<<(*par)->GetName()<<"i="<<setw(16)<<i<<":I="<<setw(16)<<intensity<<endl;
1581  #endif
1582  for(unsigned long j=nb;j>0;j--)
1583  {
1584  #ifdef PowderPatternDiffraction_CalcPowderPatternIntegrated_FullDerivDEBUG
1585  *pData += intensity * *pFact ;
1586  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;
1587  pData++;pFact++ ;
1588  #else
1589  *pData++ += intensity * *pFact++ ;
1590  #endif
1591  }
1592  ++pos;
1593  ctpar++;
1594  }
1595  }
1596  }
1597  #ifdef PowderPatternDiffraction_CalcPowderPatternIntegrated_FullDerivDEBUG
1598  std::vector<const CrystVector_REAL*> v;
1599  int n=0;
1600  cout<<"PowderPatternDiffraction::CalcPowderPatternIntegrated_FullDeriv():parameters:"<<endl;
1601  for(std::set<RefinablePar*>::iterator par=vPar.begin();par!=vPar.end();++par)
1602  {
1603  if(mPowderPatternIntegrated_FullDeriv[*par].size()==0) continue;
1604  v.push_back(&(mPowderPatternIntegrated_FullDeriv[*par]));
1605  v.push_back(&(oldDeriv[*par]));
1606  cout<<(*par)->GetName()<<":"<<mPowderPatternIntegrated_FullDeriv[*par].size()<<","<<oldDeriv[*par].size()<<endl;
1607  if(++n>6) break;
1608  }
1609  cout<<"PowderPatternDiffraction::CalcPowderPatternIntegrated_FullDeriv():"<<endl<<FormatVertVector<REAL>(v,12,1,20)<<endl;
1610  //exit(0);
1611  #endif
1612 }
1613 
1615 {
1616  this->CalcSinThetaLambda();
1618  if( (mClockProfileCalc>mClockProfilePar)
1619  &&(mClockProfileCalc>mpReflectionProfile->GetClockMaster())
1620  &&(mClockProfileCalc>mClockTheta)
1621  &&(mClockProfileCalc>this->GetRadiation().GetClockWavelength())
1622  &&(mClockProfileCalc>mpParentPowderPattern->GetClockPowderPatternXCorr())
1623  &&(mClockProfileCalc>mClockHKL)
1624  &&(mClockProfileCalc>mpParentPowderPattern->GetClockNbPointUsed())) return;
1625 
1626  TAU_PROFILE("PowderPatternDiffraction::CalcPowderReflProfile()","void (bool)",TAU_DEFAULT);
1627  VFN_DEBUG_ENTRY("PowderPatternDiffraction::CalcPowderReflProfile()",5)
1628 
1629  //Calc all profiles
1630  mvLabel.clear();
1631  stringstream label;
1632 
1633  unsigned int nbLine=1;
1634  CrystVector_REAL spectrumDeltaLambdaOvLambda;
1635  CrystVector_REAL spectrumFactor;//relative weigths of different lines of X-Ray tube
1636  switch(this->GetRadiation().GetWavelengthType())
1637  {
1638  case WAVELENGTH_MONOCHROMATIC:
1639  {
1640  spectrumDeltaLambdaOvLambda.resize(1);spectrumDeltaLambdaOvLambda=0.0;
1641  spectrumFactor.resize(1);spectrumFactor=1.0;
1642  break;
1643  }
1644  case WAVELENGTH_ALPHA12:
1645  {
1646  nbLine=2;
1647  spectrumDeltaLambdaOvLambda.resize(2);
1648  spectrumDeltaLambdaOvLambda(0)
1649  =-this->GetRadiation().GetXRayTubeDeltaLambda()
1650  *this->GetRadiation().GetXRayTubeAlpha2Alpha1Ratio()
1651  /(1+this->GetRadiation().GetXRayTubeAlpha2Alpha1Ratio())
1652  /this->GetRadiation().GetWavelength()(0);
1653  spectrumDeltaLambdaOvLambda(1)
1654  = this->GetRadiation().GetXRayTubeDeltaLambda()
1655  /(1+this->GetRadiation().GetXRayTubeAlpha2Alpha1Ratio())
1656  /this->GetRadiation().GetWavelength()(0);
1657 
1658  spectrumFactor.resize(2);
1659  spectrumFactor(0)=1./(1.+this->GetRadiation().GetXRayTubeAlpha2Alpha1Ratio());
1660  spectrumFactor(1)=this->GetRadiation().GetXRayTubeAlpha2Alpha1Ratio()
1661  /(1.+this->GetRadiation().GetXRayTubeAlpha2Alpha1Ratio());
1662  break;
1663  }
1664  case WAVELENGTH_TOF:
1665  {
1666  spectrumDeltaLambdaOvLambda.resize(1);spectrumDeltaLambdaOvLambda=0.0;
1667  spectrumFactor.resize(1);spectrumFactor=1.0;
1668  break;
1669  }
1670  default: throw ObjCrystException("PowderPatternDiffraction::PrepareIntegratedProfile():\
1671 Radiation must be either monochromatic, from an X-Ray Tube, or neutron TOF !!");
1672  }
1673 
1674 
1675  VFN_DEBUG_MESSAGE("PowderPatternDiffraction::CalcPowderReflProfile():\
1676 Computing all Profiles",5)
1677  REAL center,// center of current reflection (depends on line if several)
1678  x0; // theoretical (uncorrected for zero's, etc..) position of center of line
1679  long first,last;// first & last point of the stored profile
1680  CrystVector_REAL vx,reflProfile,tmpV;
1681  mvReflProfile.resize(this->GetNbRefl());
1682  for(unsigned int i=0;i<this->GetNbRefl();i++)
1683  {
1684  mvReflProfile[i].first=0;
1685  mvReflProfile[i].last=0;
1686  mvReflProfile[i].profile.resize(0);
1687  }
1688  VFN_DEBUG_MESSAGE("PowderPatternDiffraction::CalcPowderReflProfile()",5)
1689 
1690  for(unsigned int line=0;line<nbLine;line++)
1691  {
1692  for(long i=0;i<this->GetNbRefl();i++)
1693  {// Only the reflections contributing below the max(sin(theta)/lambda) will be computed
1694  VFN_DEBUG_ENTRY("PowderPatternDiffraction::CalcPowderReflProfile()#"<<i,5)
1695  x0=mpParentPowderPattern->STOL2X(mSinThetaLambda(i));
1696 
1697  VFN_DEBUG_MESSAGE("PowderPatternDiffraction::CalcPowderReflProfile()#"<<i,5)
1698  if(nbLine>1)
1699  {// we have several lines, not centered on the profile range
1700  center = mpParentPowderPattern->X2XCorr(
1701  x0+2*tan(x0/2.0)*spectrumDeltaLambdaOvLambda(line));
1702  }
1703  else center=mpParentPowderPattern->X2XCorr(x0);
1704  REAL fact=1.0;
1705  if(!mUseFastLessPreciseFunc) fact=5.0;
1706  const REAL halfwidth=mpReflectionProfile->GetFullProfileWidth(0.04,center,mH(i),mK(i),mL(i))*fact;
1707  if(line==0)
1708  {
1709  // For an X-Ray tube, label on first (strongest) of reflections lines (Kalpha1)
1710  label.str("");
1711  label<<mIntH(i)<<" "<<mIntK(i)<<" "<<mIntL(i);
1712  mvLabel.push_back(make_pair(center,label.str()));
1713  REAL spectrumwidth=0.0;
1714  if(this->GetRadiation().GetWavelengthType()==WAVELENGTH_ALPHA12)
1715  {// We need to shift the last point to include 2 lines in the profile
1716  spectrumwidth=2*this->GetRadiation().GetXRayTubeDeltaLambda()
1717  /this->GetRadiation().GetWavelength()(0)*tan(x0/2.0);
1718  }
1719  first=(long)(mpParentPowderPattern->X2Pixel(center-halfwidth));
1720  last =(long)(mpParentPowderPattern->X2Pixel(center+halfwidth+spectrumwidth));
1721  if(this->GetRadiation().GetWavelengthType()==WAVELENGTH_TOF)
1722  {
1723  const long f=first;
1724  first=last;
1725  last=f;
1726  }
1727  if(first>last)
1728  { // Whoops - should not happen !! Unless there is a strange (dis)order for the x coordinates...
1729  cout<<"PowderPatternDiffraction::CalcPowderReflProfile(), line"<<__LINE__<<"first>last !! :"<<first<<","<<last<<endl;
1730  first=(first+last)/2;
1731  last=first;
1732  }
1733  first -=1;
1734  last+=1;
1735  VFN_DEBUG_MESSAGE("PowderPatternDiffraction::CalcPowderReflProfile():"<<first<<","<<last<<","<<center,3)
1736  if(first>last)
1737  {
1738  cout<<__FILE__<<__LINE__<<endl;
1739  exit(0);
1740  }
1741  if((last>=0)&&(first<(long)(mpParentPowderPattern->GetNbPoint())))
1742  {
1743  if(first<0) first=0;
1744  if(last>=(long)(mpParentPowderPattern->GetNbPoint()))
1746  vx.resize(last-first+1);
1747  }
1748  else vx.resize(0); // store no profile if reflection out of pattern
1749  mvReflProfile[i].first=first;
1750  mvReflProfile[i].last=last;
1751  }
1752  else
1753  {
1754  first=mvReflProfile[i].first;
1755  last=mvReflProfile[i].last;
1756  if((last>=0)&&(first<(long)(mpParentPowderPattern->GetNbPoint())))
1757  vx.resize(last-first+1);
1758  else vx.resize(0);
1759  vx.resize(last-first+1);
1760  }
1761  if((last>=0)&&(first<(long)(mpParentPowderPattern->GetNbPoint())))
1762  {
1763  {
1764  const REAL *p0=mpParentPowderPattern->GetPowderPatternX().data()+first;
1765  REAL *p1=vx.data();
1766  for(long i=first;i<=last;i++) *p1++ = *p0++;
1767  }
1768 
1769  VFN_DEBUG_MESSAGE("PowderPatternDiffraction::CalcPowderReflProfile():"<<first<<","<<last<<","<<center,3)
1770  reflProfile=mpReflectionProfile->GetProfile(vx,center,mH(i),mK(i),mL(i));
1771  VFN_DEBUG_MESSAGE("PowderPatternDiffraction::CalcPowderReflProfile()",2)
1772  if(nbLine>1) reflProfile *=spectrumFactor(line);
1773  if(line==0) mvReflProfile[i].profile = reflProfile;
1774  else mvReflProfile[i].profile += reflProfile;
1775  }
1776  else
1777  { // reflection is out of pattern, so store no profile
1778  mvReflProfile[i].profile.resize(0);
1779  }
1780  VFN_DEBUG_EXIT("PowderPatternDiffraction::CalcPowderReflProfile():\
1781 Computing all Profiles: Reflection #"<<i,5)
1782  if(first>(long)(mpParentPowderPattern->GetNbPointUsed())) break;
1783  }
1784  }
1785  mClockProfileCalc.Click();
1786  VFN_DEBUG_EXIT("PowderPatternDiffraction::CalcPowderReflProfile()",5)
1787 }
1788 
1790 {
1791  TAU_PROFILE("PowderPatternDiffraction::CalcPowderReflProfile_FullDeriv()","void (bool)",TAU_DEFAULT);
1792  cout<<__FILE__<<":"<<__LINE__<<":PowderPatternDiffraction::CalcPowderReflProfile_FullDeriv()"<<endl;
1793  this->CalcPowderReflProfile();
1794  unsigned int nbLine=1;
1795  CrystVector_REAL spectrumDeltaLambdaOvLambda;
1796  CrystVector_REAL spectrumFactor;//relative weigths of different lines of X-Ray tube
1797  switch(this->GetRadiation().GetWavelengthType())
1798  {
1799  case WAVELENGTH_MONOCHROMATIC:
1800  {
1801  spectrumDeltaLambdaOvLambda.resize(1);spectrumDeltaLambdaOvLambda=0.0;
1802  spectrumFactor.resize(1);spectrumFactor=1.0;
1803  break;
1804  }
1805  case WAVELENGTH_ALPHA12:
1806  {
1807  nbLine=2;
1808  spectrumDeltaLambdaOvLambda.resize(2);
1809  spectrumDeltaLambdaOvLambda(0)
1810  =-this->GetRadiation().GetXRayTubeDeltaLambda()
1811  *this->GetRadiation().GetXRayTubeAlpha2Alpha1Ratio()
1812  /(1+this->GetRadiation().GetXRayTubeAlpha2Alpha1Ratio())
1813  /this->GetRadiation().GetWavelength()(0);
1814  spectrumDeltaLambdaOvLambda(1)
1815  = this->GetRadiation().GetXRayTubeDeltaLambda()
1816  /(1+this->GetRadiation().GetXRayTubeAlpha2Alpha1Ratio())
1817  /this->GetRadiation().GetWavelength()(0);
1818 
1819  spectrumFactor.resize(2);
1820  spectrumFactor(0)=1./(1.+this->GetRadiation().GetXRayTubeAlpha2Alpha1Ratio());
1821  spectrumFactor(1)=this->GetRadiation().GetXRayTubeAlpha2Alpha1Ratio()
1822  /(1.+this->GetRadiation().GetXRayTubeAlpha2Alpha1Ratio());
1823  break;
1824  }
1825  case WAVELENGTH_TOF:
1826  {
1827  spectrumDeltaLambdaOvLambda.resize(1);spectrumDeltaLambdaOvLambda=0.0;
1828  spectrumFactor.resize(1);spectrumFactor=1.0;
1829  break;
1830  }
1831  default: throw ObjCrystException("PowderPatternDiffraction::CalcPowderReflProfile_FullDeriv():\
1832 Radiation must be either monochromatic, from an X-Ray Tube, or neutron TOF !!");
1833  }
1834  REAL center,// center of current reflection (depends on line if several)
1835  x0; // theoretical (uncorrected for zero's, etc..) position of center of line
1836  long first,last;// first & last point of the stored profile
1837  CrystVector_REAL vx,reflProfile,tmpV;
1838 
1839  // Derivative vs the shift of the reflection center
1840  vector<CrystVector_REAL> vReflProfile_DerivCenter(mNbReflUsed);
1841 
1842  mvReflProfile_FullDeriv.clear();
1843  for(std::set<RefinablePar*>::iterator par=vPar.begin();par!=vPar.end();++par)
1844  {
1845  if(*par==0) continue;
1846  if( (*par)->GetType()->IsDescendantFromOrSameAs(gpRefParTypeRadiation)
1847  ||(*par)->GetType()->IsDescendantFromOrSameAs(gpRefParTypeUnitCell)
1848  ||(*par)->GetType()->IsDescendantFromOrSameAs(gpRefParTypeScattDataCorrPos)
1849  ||(*par)->GetType()->IsDescendantFromOrSameAs(gpRefParTypeScattDataProfile))
1850  {
1851  mvReflProfile_FullDeriv[*par].resize(mNbReflUsed);
1852 
1853  for(unsigned int line=0;line<nbLine;line++)
1854  {
1855  for(long i=0;i<mNbReflUsed;i++)
1856  {
1857  x0=mpParentPowderPattern->STOL2X(mSinThetaLambda(i));
1858 
1859  if(nbLine>1)
1860  {// we have several lines, not centered on the profile range
1861  center = mpParentPowderPattern->X2XCorr(
1862  x0+2*tan(x0/2.0)*spectrumDeltaLambdaOvLambda(line));
1863  }
1864  else center=mpParentPowderPattern->X2XCorr(x0);
1865  REAL fact=1.0;
1866  if(!mUseFastLessPreciseFunc) fact=5.0;
1867 
1868  first=mvReflProfile[i].first;
1869  last=mvReflProfile[i].last;
1870  if((last>=0)&&(first<(long)(mpParentPowderPattern->GetNbPoint())))
1871  vx.resize(last-first+1);
1872  else vx.resize(0);
1873  vx.resize(last-first+1);
1874  if((last>=0)&&(first<(long)(mpParentPowderPattern->GetNbPoint())))
1875  {
1876  {
1877  const REAL *p0=mpParentPowderPattern->GetPowderPatternX().data()+first;
1878  REAL *p1=vx.data();
1879  for(long i=first;i<=last;i++) *p1++ = *p0++;
1880  }
1881 
1882  if((*par)->GetType()->IsDescendantFromOrSameAs(gpRefParTypeScattDataProfile))
1883  {// Parameter only affects profile
1884  //if(i==0) cout<<"PowderPatternDiffraction::CalcPowderReflProfile_FullDeriv()par="<<(*par)->GetName()<<":refl #"<<i<<endl;
1885  //:TODO: analytical derivatives
1886  const REAL step=(*par)->GetDerivStep();
1887  (*par)->Mutate(step);
1888  reflProfile=mpReflectionProfile->GetProfile(vx,center,mH(i),mK(i),mL(i));
1889  (*par)->Mutate(-2*step);
1890  reflProfile-=mpReflectionProfile->GetProfile(vx,center,mH(i),mK(i),mL(i));
1891  (*par)->Mutate(step);
1892  reflProfile/=2*step;
1893  }
1894  else
1895  {// Parameter affects reflection center
1896  REAL dcenter=0;
1897  {
1898  //:TODO: analytical derivatives
1899  const REAL step=(*par)->GetDerivStep();
1900  (*par)->Mutate(step);
1901  REAL x1=mpParentPowderPattern->STOL2X(this->CalcSinThetaLambda(mH(i),mK(i),mL(i)));
1902  if(nbLine>1) dcenter = mpParentPowderPattern->X2XCorr(x1+2*tan(x1/2.0)*spectrumDeltaLambdaOvLambda(line));
1903  else dcenter = mpParentPowderPattern->X2XCorr(x1);
1904  (*par)->Mutate(-2*step);
1905  x1=mpParentPowderPattern->STOL2X(this->CalcSinThetaLambda(mH(i),mK(i),mL(i)));
1906  if(nbLine>1) dcenter-= mpParentPowderPattern->X2XCorr(x1+2*tan(x1/2.0)*spectrumDeltaLambdaOvLambda(line));
1907  else dcenter-= mpParentPowderPattern->X2XCorr(x1);
1908  (*par)->Mutate(step);
1909  dcenter/=2*step;
1910  }
1911 
1912  if(dcenter!=0)
1913  {
1914  //if(i==0) cout<<"PowderPatternDiffraction::CalcPowderReflProfile_FullDeriv()par="<<(*par)->GetName()<<":refl #"<<i<<", dcenter="<<setw(8)<<dcenter<<endl;
1915  if(vReflProfile_DerivCenter[i].size()==0)
1916  {
1917  const REAL step=1e-4;//:TODO: adapt for TOF
1918  vReflProfile_DerivCenter[i] =mpReflectionProfile->GetProfile(vx,center+step,mH(i),mK(i),mL(i));
1919  vReflProfile_DerivCenter[i]-=mpReflectionProfile->GetProfile(vx,center-step,mH(i),mK(i),mL(i));
1920  vReflProfile_DerivCenter[i]/=2*step;
1921  }
1922  reflProfile=vReflProfile_DerivCenter[i];
1923  reflProfile*=dcenter;
1924  }
1925  else
1926  {
1927  //if(i==0) cout<<"PowderPatternDiffraction::CalcPowderReflProfile_FullDeriv()par="<<(*par)->GetName()<<":refl #"<<i<<" => Parameter affects nothing ?"<<endl;
1928  reflProfile.resize(0);
1929  }
1930  }
1931  if(reflProfile.size()>0)
1932  {
1933  if(nbLine>1) reflProfile *=spectrumFactor(line);
1934  if(line==0) mvReflProfile_FullDeriv[*par][i] = reflProfile;
1935  else mvReflProfile_FullDeriv[*par][i] += reflProfile;
1936  }
1937  }
1938  }
1939  }
1940  }
1941  }
1942 }
1943 
1945 {
1946  bool needRecalc=false;
1947 
1948  this->CalcSinThetaLambda();
1949  if((mClockIntensityCorr<mClockTheta)||(mClockIntensityCorr<this->GetClockNbReflBelowMaxSinThetaOvLambda())) needRecalc=true;
1950 
1951  const CrystVector_REAL *mpCorr[5] = {0, 0, 0, 0, 0};
1952 
1953  if(this->GetRadiation().GetWavelengthType()==WAVELENGTH_TOF)
1954  {
1955  mpCorr[0]=&(mCorrTOF.GetCorr());
1956  if(mClockIntensityCorr<mCorrTOF.GetClockCorr()) needRecalc=true;
1957  }
1958  else
1959  {
1960  mpCorr[0]=&(mCorrLorentz.GetCorr());
1961  if(mClockIntensityCorr<mCorrLorentz.GetClockCorr()) needRecalc=true;
1962 
1963  if(this->GetRadiation().GetRadiationType()==RAD_XRAY)
1964  {
1965  mpCorr[1]=&(mCorrPolar.GetCorr());
1966  if(mClockIntensityCorr<mCorrPolar.GetClockCorr()) needRecalc=true;
1967  }
1968 
1969  mpCorr[2]=&(mCorrSlitAperture.GetCorr());
1970  if(mClockIntensityCorr<mCorrSlitAperture.GetClockCorr()) needRecalc=true;
1971  }
1972 
1973  if(mCorrTextureMarchDollase.GetNbPhase()>0)
1974  {
1975  mpCorr[3]=&(mCorrTextureMarchDollase.GetCorr());
1976  if(mClockIntensityCorr<mCorrTextureMarchDollase.GetClockCorr()) needRecalc=true;
1977  }
1978  mpCorr[4]=&(mCorrTextureEllipsoid.GetCorr());
1979  if(mClockIntensityCorr<mCorrTextureEllipsoid.GetClockCorr()) needRecalc=true;
1980 
1981 
1982  if(needRecalc==false) return;
1983 
1984  TAU_PROFILE("PowderPatternDiffraction::CalcIntensityCorr()","void ()",TAU_DEFAULT);
1985  VFN_DEBUG_MESSAGE("PowderPatternDiffraction::CalcIntensityCorr()",2)
1986  mIntensityCorr.resize(mNbRefl);
1987  REAL *pCorr=mIntensityCorr.data();
1988  const REAL *p=mpCorr[0]->data();
1989  for(long i=mNbReflUsed;i>0;i--) *pCorr++ = *p++;
1990  if(this->GetRadiation().GetWavelengthType()!=WAVELENGTH_TOF)
1991  {
1992  if(this->GetRadiation().GetRadiationType()==RAD_XRAY)
1993  {
1994  pCorr=mIntensityCorr.data();
1995  p=mpCorr[1]->data();
1996  const REAL* p2=mpCorr[2]->data();
1997  for(long i=mNbReflUsed;i>0;i--) *pCorr++ *= *p++ * *p2++;
1998  }
1999  else
2000  {
2001  pCorr=mIntensityCorr.data();
2002  p=mpCorr[2]->data();
2003  for(long i=mNbReflUsed;i>0;i--) *pCorr++ *= *p++;
2004  }
2005  }
2006  if(mCorrTextureMarchDollase.GetNbPhase()>0)
2007  {
2008  pCorr=mIntensityCorr.data();
2009  p=mpCorr[3]->data();
2010  for(long i=mNbReflUsed;i>0;i--) *pCorr++ *= *p++;
2011  }
2012  if(mpCorr[4]->numElements()>0)
2013  {
2014  pCorr=mIntensityCorr.data();
2015  p=mpCorr[4]->data();
2016  for(long i=mNbReflUsed;i>0;i--) *pCorr++ *= *p++;
2017  }
2018  mClockIntensityCorr.Click();
2019  VFN_DEBUG_MESSAGE("PowderPatternDiffraction::CalcIntensityCorr():finished",2)
2020 }
2021 
2023 {
2024  this->CalcIntensityCorr();
2025  if(mExtractionMode==true)
2026  {
2027  VFN_DEBUG_MESSAGE("PowderPatternDiffraction::CalcIhkl():"<<mFhklObsSq.numElements()<<","<<mIntensityCorr.numElements()<<","<<mMultiplicity.numElements(),7);
2028  mIhklCalc=mFhklObsSq;
2029  mIhklCalc*=mIntensityCorr;
2030  mIhklCalc*=mMultiplicity;
2031  mClockIhklCalc.Click();
2032  return;
2033  }
2034  this->CalcStructFactor();
2035  if( (mClockIhklCalc>mClockIntensityCorr)
2036  &&(mClockIhklCalc>mClockStructFactor)
2037  &&(mClockIhklCalc>mClockNbReflUsed)) return;
2038 
2039  VFN_DEBUG_MESSAGE("PowderPatternDiffraction::CalcIhkl()",3)
2040  TAU_PROFILE("PowderPatternDiffraction::CalcIhkl()","void ()",TAU_DEFAULT);
2041  const REAL * RESTRICT pr,* RESTRICT pi,* RESTRICT pcorr;
2042  const int * RESTRICT mult;
2043  REAL * RESTRICT p;
2044 
2045  pr=mFhklCalcReal.data();
2046  pi=mFhklCalcImag.data();
2047  pcorr=mIntensityCorr.data();
2048 
2049  mult=mMultiplicity.data();
2050  mIhklCalc.resize(mNbRefl);
2051  p=mIhklCalc.data();
2052  if(mFhklCalcVariance.numElements()>0)
2053  {
2054  const REAL * RESTRICT pv=mFhklCalcVariance.data();
2055  for(long i=mNbReflUsed;i>0;i--)
2056  {
2057  *p++ = *mult++ * (*pr * *pr + *pi * *pi + 2 * *pv++) * *pcorr++;
2058  pr++;
2059  pi++;
2060  }
2061  }
2062  else
2063  {
2064  for(long i=mNbReflUsed;i>0;i--)
2065  {
2066  *p++ = *mult++ * (*pr * *pr + *pi * *pi) * *pcorr++;
2067  pr++;
2068  pi++;
2069  }
2070  }
2071  if(mFhklCalcVariance.numElements()==0)
2072  {
2073  VFN_DEBUG_MESSAGE("PowderPatternDiffraction::CalcIhkl(): No Calc Variance",2)
2074  mIhklCalcVariance.resize(0);
2075  VFN_DEBUG_MESSAGE(endl<<
2076  FormatVertVectorHKLFloats<REAL>(mH,mK,mL,mSinThetaLambda,
2077  mFhklCalcReal,
2078  mFhklCalcImag,
2079  mIhklCalc,
2080  mIntensityCorr
2081  ),2)
2082  }
2083  else
2084  {
2085  VFN_DEBUG_MESSAGE("PowderPatternDiffraction::CalcIhkl(): Calc Variance",2)
2086  mIhklCalcVariance.resize(mNbRefl);
2087  REAL * RESTRICT pVar2=mIhklCalcVariance.data();
2088 
2089  const REAL * RESTRICT pInt=mIhklCalc.data();
2090  const REAL * RESTRICT pVar=mFhklCalcVariance.data();
2091  pcorr=mIntensityCorr.data();
2092  mult=mMultiplicity.data();
2093 
2094  for(long j=mNbReflUsed;j>0;j--)
2095  {
2096  *pVar2++ = (4* *mult) * *pcorr * *pVar *(*pInt++ - (*mult * *pcorr) * *pVar);
2097  pVar++;mult++;pcorr++;
2098  }
2099  VFN_DEBUG_MESSAGE(endl<<
2100  FormatVertVectorHKLFloats<REAL>(mH,mK,mL,mSinThetaLambda,
2101  mFhklCalcReal,
2102  mIhklCalc,
2103  mExpectedIntensityFactor,
2104  mIntensityCorr,
2105  mMultiplicity,
2106  mvLuzzatiFactor[&(mpCrystal->GetScatteringPowerRegistry().GetObj(0))],
2107  mFhklCalcVariance,
2108  mIhklCalcVariance),2);
2109  VFN_DEBUG_MESSAGE(mNbRefl<<" "<<mNbReflUsed,2)
2110  }
2111 
2112  //cout <<FormatVertVector<REAL>(mTheta,mIhklCalc,mMultiplicity,
2113  // mFhklCalcReal,mFhklCalcImag,mIntensityCorr);
2114  mClockIhklCalc.Click();
2115  VFN_DEBUG_MESSAGE("PowderPatternDiffraction::CalcIhkl():End",3)
2116 }
2117 
2118 void PowderPatternDiffraction::CalcIhkl_FullDeriv(std::set<RefinablePar*> &vPar)
2119 {
2120  TAU_PROFILE("PowderPatternDiffraction::CalcIhkl_FullDeriv()","void ()",TAU_DEFAULT);
2121  //cout<<"PowderPatternDiffraction::CalcIhkl_FullDeriv()"<<endl;
2122  this->CalcIntensityCorr();//:TODO: derivatives of intensity corrections (Texture, displacement parameters,...)
2123  mIhkl_FullDeriv.clear();
2124  if(mExtractionMode==true)
2125  {
2126  //:TODO: handle Pawley refinements of I(hkl)
2127  return;
2128  }
2129  this->CalcStructFactor_FullDeriv(vPar);
2130 
2131  for(std::set<RefinablePar*>::iterator par=vPar.begin();par!=vPar.end();++par)
2132  {
2133  if(*par==0) mIhkl_FullDeriv[*par]=mIhklCalc;
2134  else
2135  {
2136  if(mFhklCalcReal_FullDeriv[*par].size()==0) continue;
2137  const REAL * RESTRICT pr,* RESTRICT pi,* RESTRICT prd,* RESTRICT pid,* RESTRICT pcorr;
2138  const int * RESTRICT mult;
2139  REAL * RESTRICT p;
2140 
2141  pr=mFhklCalcReal.data();
2142  pi=mFhklCalcImag.data();
2143  prd=mFhklCalcReal_FullDeriv[*par].data();
2144  pid=mFhklCalcImag_FullDeriv[*par].data();
2145  pcorr=mIntensityCorr.data();//:TODO: derivatives of intensity corrections (Texture, displacement parameters,...)
2146 
2147  mult=mMultiplicity.data();
2148  mIhkl_FullDeriv[*par].resize(mNbRefl);
2149  p=mIhkl_FullDeriv[*par].data();
2150  if(mFhklCalcImag_FullDeriv[*par].size()==0)
2151  for(long i=mNbReflUsed;i>0;i--) *p++ = *mult++ * 2 * *pr++ * *prd++ * *pcorr++;
2152  else
2153  for(long i=mNbReflUsed;i>0;i--) *p++ = *mult++ * 2 *(*pr++ * *prd++ + *pi++ * *pid++) * *pcorr++;
2154  }
2155  }
2156  #if 0
2157  std::map<RefinablePar*, CrystVector_REAL> oldDeriv;
2158  std::vector<const CrystVector_REAL*> v;
2159  v.push_back(&mH);
2160  v.push_back(&mK);
2161  v.push_back(&mL);
2162  CrystVector_REAL m;
2163  m=mMultiplicity;
2164  v.push_back(&m);
2165  v.push_back(&mIntensityCorr);
2166  int n=0;
2167  for(std::set<RefinablePar*>::iterator par=vPar.begin();par!=vPar.end();++par)
2168  {
2169  if((*par)==0) continue;
2170  if(mIhkl_FullDeriv[*par].size()==0) continue;
2171 
2172  const REAL step=(*par)->GetDerivStep();
2173  (*par)->Mutate(step);
2174  this->CalcIhkl();
2175  oldDeriv[*par]=mIhklCalc;
2176  (*par)->Mutate(-2*step);
2177  this->CalcIhkl();
2178  oldDeriv[*par]-=mIhklCalc;
2179  oldDeriv[*par]/=2*step;
2180  (*par)->Mutate(step);
2181 
2182  v.push_back(&(mIhkl_FullDeriv[*par]));
2183  v.push_back(&(oldDeriv[*par]));
2184  cout<<(*par)->GetName()<<":"<<mIhkl_FullDeriv[*par].size()<<","<<oldDeriv[*par].size()<<", step="<<setw(16)<<step<<endl;
2185  if(++n>5) break;
2186  }
2187  cout<<"PowderPatternDiffraction::CalcIhkl_FullDeriv():"<<endl<<FormatVertVectorHKLFloats<REAL>(v,12,1,20)<<endl;
2188  //exit(0);
2189  #endif
2190 }
2191 
2193 {
2194  if( (this->GetCrystal().GetSpaceGroup().GetClockSpaceGroup()>mClockHKL)
2195  ||(this->GetCrystal().GetClockLatticePar()>mClockHKL)
2196  ||(this->GetRadiation().GetClockWavelength()>mClockHKL)
2198  this->GenHKLFullSpace();
2199  //if(0==this->GetNbRefl()) this->GenHKLFullSpace();
2200 }
2201 void PowderPatternDiffraction::InitOptions()
2202 {
2203  VFN_DEBUG_MESSAGE("PowderPatternDiffraction::InitOptions()",5)
2204  #if 0
2205  static string ReflectionProfileTypeName;
2206  static string ReflectionProfileTypeChoices[3];
2207 
2208  static bool needInitNames=true;
2209  if(true==needInitNames)
2210  {
2211  ReflectionProfileTypeName="Profile Type";
2212  ReflectionProfileTypeChoices[0]="Gaussian";
2213  ReflectionProfileTypeChoices[1]="Lorentzian";
2214  ReflectionProfileTypeChoices[2]="Pseudo-Voigt";
2215 
2216  needInitNames=false;//Only once for the class
2217  }
2218  mReflectionProfileType.Init(3,&ReflectionProfileTypeName,ReflectionProfileTypeChoices);
2219  this->AddOption(&mReflectionProfileType);
2220  #endif
2221 }
2222 const CrystVector_long& PowderPatternDiffraction::GetBraggLimits()const
2223 {
2224  this->CalcPowderReflProfile();
2225  if((mClockProfileCalc>mClockBraggLimits)&&(this->GetNbReflBelowMaxSinThetaOvLambda()>0))
2226  {
2227  VFN_DEBUG_ENTRY("PowderPatternDiffraction::GetBraggLimits(*min,*max)",3)
2228  TAU_PROFILE("PowderPatternDiffraction::GetBraggLimits()","void ()",TAU_DEFAULT);
2229  mIntegratedReflLimits.resize(this->GetNbReflBelowMaxSinThetaOvLambda());
2230  long i = 0;
2231  mIntegratedReflLimits(i)=mvReflProfile[0].first;
2232  for(;i<(this->GetNbReflBelowMaxSinThetaOvLambda()-1);++i)
2233  mIntegratedReflLimits(i+1)=(mvReflProfile[i].first+mvReflProfile[i].last+mvReflProfile[i+1].first+mvReflProfile[i+1].last)/4;
2234  mIntegratedReflLimits(i)=mvReflProfile[i].last;
2236  VFN_DEBUG_EXIT("PowderPatternDiffraction::GetBraggLimits(*min,*max)",3)
2237  }
2238  return mIntegratedReflLimits;
2239 }
2240 
2243 
2244 const CrystMatrix_REAL& PowderPatternDiffraction::GetBMatrix()const
2245 {
2246  if(mFreezeLatticePar) return mFrozenBMatrix;
2247  return this->ScatteringData::GetBMatrix();
2248 }
2249 
2251 {
2252  VFN_DEBUG_MESSAGE("PowderPatternDiffraction::CalcFrozenBMatrix()", 10)
2253  REAL a,b,c,alpha,beta,gamma;//direct space parameters
2254  REAL aa,bb,cc,alphaa,betaa,gammaa;//reciprocal space parameters
2255  REAL v;//volume of the unit cell
2256  a=mFrozenLatticePar(0);
2257  b=mFrozenLatticePar(1);
2258  c=mFrozenLatticePar(2);
2259  alpha=mFrozenLatticePar(3);
2260  beta=mFrozenLatticePar(4);
2261  gamma=mFrozenLatticePar(5);
2262 
2263  v=sqrt(1-cos(alpha)*cos(alpha)-cos(beta)*cos(beta)-cos(gamma)*cos(gamma)
2264  +2*cos(alpha)*cos(beta)*cos(gamma));
2265 
2266  aa=sin(alpha)/a/v;
2267  bb=sin(beta )/b/v;
2268  cc=sin(gamma)/c/v;
2269 
2270  alphaa=acos( (cos(beta )*cos(gamma)-cos(alpha))/sin(beta )/sin(gamma) );
2271  betaa =acos( (cos(alpha)*cos(gamma)-cos(beta ))/sin(alpha)/sin(gamma) );
2272  gammaa=acos( (cos(alpha)*cos(beta )-cos(gamma))/sin(alpha)/sin(beta ) );
2273 
2274  mFrozenBMatrix = aa , bb*cos(gammaa) , cc*cos(betaa) ,
2275  0 , bb*sin(gammaa) ,-cc*sin(betaa)*cos(alpha),
2276  0 , 0 ,1/c;
2277 }
2278 
2279 void PowderPatternDiffraction::PrepareIntegratedProfile()const
2280 {
2281  this->CalcPowderReflProfile();
2282 
2283  if( (mClockIntegratedProfileFactor>mClockProfileCalc)
2284  &&(mClockIntegratedProfileFactor>mpParentPowderPattern->GetIntegratedProfileLimitsClock())
2285  &&(mClockIntegratedProfileFactor>mClockNbReflUsed))
2286  return;
2287  VFN_DEBUG_ENTRY("PowderPatternDiffraction::PrepareIntegratedProfile()",7)
2288  TAU_PROFILE("PowderPatternDiffraction::PrepareIntegratedProfile()","void ()",TAU_DEFAULT);
2289  const CrystVector_long *pMin=&(mpParentPowderPattern->GetIntegratedProfileMin());
2290  const CrystVector_long *pMax=&(mpParentPowderPattern->GetIntegratedProfileMax());
2291 
2292  const long numInterval=pMin->numElements();
2293 
2294  vector< map<long, REAL> > vIntegratedProfileFactor;
2295  vIntegratedProfileFactor.resize(mNbReflUsed);
2296  vector< map<long, REAL> >::iterator pos1;
2297  pos1=vIntegratedProfileFactor.begin();
2298 
2299  mIntegratedProfileFactor.resize(mNbReflUsed);
2300  vector< pair<unsigned long, CrystVector_REAL> >::iterator pos2;
2301  pos2=mIntegratedProfileFactor.begin();
2302  for(long i=0;i<mNbReflUsed;i++)
2303  {
2304  pos1->clear();
2305  long firstInterval=numInterval;
2306  for(long j=0;j<numInterval;j++)
2307  {
2308  const long first0 = mvReflProfile[i].first;
2309  const long last0 = mvReflProfile[i].last ;
2310  const long first= first0>(*pMin)(j) ? first0:(*pMin)(j);
2311  const long last = last0 <(*pMax)(j) ? last0 :(*pMax)(j);
2312  if((first<=last) && (mvReflProfile[i].profile.size()>0))
2313  {
2314  if(firstInterval>j) firstInterval=j;
2315  if(pos1->find(j) == pos1->end()) (*pos1)[j]=0.;
2316  REAL *fact = &((*pos1)[j]);//this creates the 'j' entry if necessary
2317  const REAL *p2 = mvReflProfile[i].profile.data()+(first-first0);
2318  //cout << i<<","<<j<<","<<first<<","<<last<<":"<<*fact<<"/"<<mNbReflUsed<<","<<mNbRefl<<endl;
2319  for(int k=first;k<=last;k++) *fact += *p2++;
2320  }
2321  }
2322  pos2->first=firstInterval;
2323  pos2->second.resize(pos1->size());
2324  REAL *pFact=pos2->second.data();
2325  for(map<long, REAL>::const_iterator pos=pos1->begin();pos!=pos1->end();++pos)
2326  *pFact++ = pos->second;
2327  pos1++;
2328  pos2++;
2329  }
2330  mClockIntegratedProfileFactor.Click();
2331  #ifdef __DEBUG__
2332  if(gVFNDebugMessageLevel<3)
2333  {
2334  unsigned long i=0;
2335  for(vector< pair<unsigned long, CrystVector_REAL> >::const_iterator
2336  pos=mIntegratedProfileFactor.begin();
2337  pos!=mIntegratedProfileFactor.end();++pos)
2338  {
2339  cout <<"Integrated profile factors for reflection #"<<i++<<" ";
2340  for(int j=0;j<pos->second.numElements();++j)
2341  cout << j+pos->first<<"("<<pos->second(j)<<") ";
2342  cout<<endl;
2343  }
2344  }
2345  #endif
2346 
2347  VFN_DEBUG_EXIT("PowderPatternDiffraction::PrepareIntegratedProfile()",7)
2348 }
2349 
2350 #ifdef __WX__CRYST__
2351 WXCrystObjBasic* PowderPatternDiffraction::WXCreate(wxWindow* parent)
2352 {
2353  //:TODO: Check mpWXCrystObj==0
2354  mpWXCrystObj=new WXPowderPatternDiffraction(parent,this);
2355  return mpWXCrystObj;
2356 }
2357 #endif
2358 
2360 //
2361 // PowderPattern
2362 //
2364 ObjRegistry<PowderPattern>
2365  gPowderPatternRegistry("List of all PowderPattern objects");
2366 
2367 PowderPattern::PowderPattern():
2368 mIsXAscending(true),mNbPoint(0),
2369 mXZero(0.),m2ThetaDisplacement(0.),m2ThetaTransparency(0.),
2370 mDIFC(48277.14),mDIFA(-6.7),
2371 mScaleFactor(20),mUseFastLessPreciseFunc(false),
2372 mStatisticsExcludeBackground(false),mMaxSinThetaOvLambda(10),mNbPointUsed(0)
2373 {
2374  mScaleFactor=1;
2375  mSubObjRegistry.SetName("SubObjRegistry for a PowderPattern object");
2376  mPowderPatternComponentRegistry.SetName("Powder Pattern Components");
2377  this->AddSubRefObj(mRadiation);
2378  this->Init();
2379  gPowderPatternRegistry.Register(*this);
2380  gTopRefinableObjRegistry.Register(*this);
2381  mClockMaster.AddChild(mClockPowderPatternPar);
2382  mClockMaster.AddChild(mClockNbPointUsed);
2383  mClockMaster.AddChild(mClockPowderPatternXCorr);
2384  mClockMaster.AddChild(mClockScaleFactor);
2385  mClockMaster.AddChild(mClockPowderPatternRadiation);
2386 }
2387 
2388 PowderPattern::PowderPattern(const PowderPattern &old):
2389 mIsXAscending(old.mIsXAscending),mNbPoint(old.mNbPoint),
2390 mRadiation(old.mRadiation),
2391 mXZero(old.mXZero),m2ThetaDisplacement(old.m2ThetaDisplacement),
2392 m2ThetaTransparency(old.m2ThetaTransparency),
2393 mDIFC(old.mDIFC),mDIFA(old.mDIFA),
2394 mPowderPatternComponentRegistry(old.mPowderPatternComponentRegistry),
2395 mScaleFactor(old.mScaleFactor),
2396 mUseFastLessPreciseFunc(old.mUseFastLessPreciseFunc),
2397 mStatisticsExcludeBackground(old.mStatisticsExcludeBackground),
2398 mMaxSinThetaOvLambda(old.mMaxSinThetaOvLambda),mNbPointUsed(old.mNbPointUsed)
2399 {
2400  mX=old.mX;
2401  this->Init();
2402  mSubObjRegistry.SetName("SubObjRegistry for a PowderPattern :"+mName);
2403  gPowderPatternRegistry.Register(*this);
2404  gTopRefinableObjRegistry.Register(*this);
2405  this->AddSubRefObj(mRadiation);
2406  mClockMaster.AddChild(mClockPowderPatternPar);
2407  mClockMaster.AddChild(mClockNbPointUsed);
2408  mClockMaster.AddChild(mClockPowderPatternXCorr);
2409  mClockMaster.AddChild(mClockScaleFactor);
2410  mClockMaster.AddChild(mClockPowderPatternRadiation);
2411 }
2412 
2413 PowderPattern::~PowderPattern()
2414 {
2415  for(int i=0;i<mPowderPatternComponentRegistry.GetNb();i++)
2416  {
2417  mPowderPatternComponentRegistry.GetObj(i).DeRegisterClient(*this);
2418  this->RemoveSubRefObj(mPowderPatternComponentRegistry.GetObj(i));
2419  delete &(mPowderPatternComponentRegistry.GetObj(i));
2420  }
2421  gPowderPatternRegistry.DeRegister(*this);
2422  gTopRefinableObjRegistry.DeRegister(*this);
2423 }
2424 const string& PowderPattern::GetClassName() const
2425 {
2426  const static string className="PowderPattern";
2427  return className;
2428 }
2429 
2431 {
2432  VFN_DEBUG_ENTRY("PowderPattern::AddPowderPatternComponent():"<<comp.GetName(),5)
2433  comp.SetParentPowderPattern(*this);
2434  this->AddSubRefObj(comp);
2435  comp.RegisterClient(*this);
2437  mClockIntegratedFactorsPrep.Reset();
2438  mPowderPatternComponentRegistry.Register(comp);
2439  //:TODO: check if there are enough scale factors
2440  //mScaleFactor.resizeAndPreserve(mPowderPatternComponentRegistry.GetNb());
2441  mScaleFactor(mPowderPatternComponentRegistry.GetNb()-1)=1.;
2443  if(comp.IsScalable())
2444  {//Init refinable parameter
2445  RefinablePar tmp("Scale_"+comp.GetName(),mScaleFactor.data()+mPowderPatternComponentRegistry.GetNb()-1,
2446  1e-10,1e10,gpRefParTypeScattDataScale,REFPAR_DERIV_STEP_RELATIVE,
2447  false,true,true,false,1.);
2448  tmp.SetGlobalOptimStep(0.);
2450  tmp.SetDerivStep(1e-4);
2451  this->AddPar(tmp);
2452  }
2453  //this->UpdateDisplay();
2454  VFN_DEBUG_EXIT("PowderPattern::AddPowderPatternComponent():"<<comp.GetName(),5)
2455 }
2456 
2458 {
2459  return mPowderPatternComponentRegistry.GetNb();
2460 }
2461 
2463  (const string &name)const
2464 {
2465  return mPowderPatternComponentRegistry.GetObj(name);
2466 }
2467 
2469  (const int i) const
2470 {
2471  return mPowderPatternComponentRegistry.GetObj(i);
2472 }
2473 
2475  (const string &name)
2476 {
2477  return mPowderPatternComponentRegistry.GetObj(name);
2478 }
2479 
2481  (const int i)
2482 {
2483  return mPowderPatternComponentRegistry.GetObj(i);
2484 }
2485 
2486 REAL PowderPattern::GetScaleFactor(const int i)const{return mScaleFactor(i);}
2487 
2489 {
2490  unsigned int i=0;
2491  for(;i<mPowderPatternComponentRegistry.GetNb();++i)
2492  {
2493  if(&(mPowderPatternComponentRegistry.GetObj(i))==&comp) break;
2494  }
2495  if(i==mPowderPatternComponentRegistry.GetNb())
2496  throw ObjCrystException("PowderPattern::GetScaleFactor(comp) : no such component");
2497  return mScaleFactor(i);
2498 }
2499 void PowderPattern::SetScaleFactor(const int i, REAL s){ mScaleFactor(i)=s;}
2500 
2502 {
2503  unsigned int i=0;
2504  for(;i<mPowderPatternComponentRegistry.GetNb();++i)
2505  {
2506  if(&(mPowderPatternComponentRegistry.GetObj(i))==&comp) break;
2507  }
2508  if(i==mPowderPatternComponentRegistry.GetNb())
2509  throw ObjCrystException("PowderPattern::GetScaleFactor(comp) : no such component");
2510  mScaleFactor(i)=s;
2511 }
2512 
2514  const REAL step,
2515  unsigned long nbPoint)
2516 {
2517  VFN_DEBUG_MESSAGE("PowderPattern::SetPowderPatternPar():"<<min<<","<<step<<","<<nbPoint,3)
2518  mNbPoint=nbPoint;
2519  mX.resize(mNbPoint);
2520  for(unsigned long i=0;i<mNbPoint;i++) mX(i)=min+step*i;
2521  mPowderPatternObs.resizeAndPreserve(mNbPoint);
2522  mPowderPatternObsSigma.resizeAndPreserve(mNbPoint);
2523  mPowderPatternWeight.resizeAndPreserve(mNbPoint);
2525 }
2526 void PowderPattern::SetPowderPatternX(const CrystVector_REAL &x)
2527 {
2528  mNbPoint=x.numElements();
2529  if(&x != &mX) mX=x;
2530  mPowderPatternObs.resizeAndPreserve(mNbPoint);
2531  mPowderPatternObsSigma.resizeAndPreserve(mNbPoint);
2532  mPowderPatternWeight.resizeAndPreserve(mNbPoint);
2534  if(mX(mNbPoint-1)>mX(0))mIsXAscending=true;
2535  else mIsXAscending=false;
2536  VFN_DEBUG_MESSAGE("PowderPattern::SetPowderPatternX() is ascending="<<mIsXAscending,5)
2537 }
2538 
2539 unsigned long PowderPattern::GetNbPoint()const {return mNbPoint;}
2540 
2541 unsigned long PowderPattern::GetNbPointUsed()const
2542 {
2543  if(!this->IsBeingRefined()) this->CalcNbPointUsed();
2544  return mNbPointUsed;
2545 }
2546 
2547 const RefinableObjClock& PowderPattern::GetClockNbPointUsed()const{return mClockNbPointUsed;}
2548 
2550 {
2551  mRadiation=radiation;
2553 }
2555 
2557 
2559 {
2561 }
2562 
2564 void PowderPattern::SetWavelength(const REAL lambda)
2565 {
2566  VFN_DEBUG_MESSAGE("PowderPattern::SetWavelength(lambda)",3)
2567  mRadiation.SetWavelength(lambda);
2568 }
2569 
2570 void PowderPattern::SetWavelength(const string &XRayTubeElementName,const REAL alpha12ratio)
2571 {
2572  VFN_DEBUG_MESSAGE("PowderPattern::SetWavelength(wavelength)",3)
2573  mRadiation.SetWavelength(XRayTubeElementName,alpha12ratio);
2574 }
2575 
2577 
2578 const CrystVector_REAL& PowderPattern::GetPowderPatternCalc()const
2579 {
2580  this->CalcPowderPattern();
2581  return mPowderPatternCalc;
2582 }
2583 
2584 std::map<RefinablePar*,CrystVector_REAL>& PowderPattern::GetPowderPattern_FullDeriv(std::set<RefinablePar *> &vPar)
2585 {
2586  this->CalcPowderPattern_FullDeriv(vPar);
2587  return mPowderPattern_FullDeriv;
2588 }
2589 
2590 const CrystVector_REAL& PowderPattern::GetPowderPatternObs()const
2591 {
2592  return mPowderPatternObs;
2593 }
2594 
2595 const CrystVector_REAL& PowderPattern::GetPowderPatternObsSigma()const
2596 {
2597  return mPowderPatternObsSigma;
2598 }
2599 
2600 const CrystVector_REAL& PowderPattern::GetPowderPatternVariance()const
2601 {
2602  return mPowderPatternVariance;
2603 }
2604 
2605 const CrystVector_REAL& PowderPattern::GetPowderPatternWeight()const
2606 {
2607  return mPowderPatternWeight;
2608 }
2609 
2611 {
2612  if(mNbPoint==0) return 0;//:KLUDGE: ?
2613  if(true==mIsXAscending) return mX(0);
2614  return mX(mNbPoint-1);
2615 }
2616 
2618 {
2619  if(mNbPoint==0) return 0;//:KLUDGE: ?
2620  return abs((-mX(0)+mX(mNbPoint-1))/(mNbPoint-1));
2621 }
2622 
2624 {
2625  if(mNbPoint==0) return 0;//:KLUDGE: ?
2626  if(true==mIsXAscending) return mX(mNbPoint-1);
2627  return mX(0);
2628 }
2629 const CrystVector_REAL& PowderPattern::GetPowderPatternX()const
2630 {
2631  return mX;
2632 }
2633 
2634 const CrystVector_REAL& PowderPattern::GetChi2Cumul()const
2635 {
2636  VFN_DEBUG_ENTRY("PowderPattern::GetChi2Cumul()",3)
2637  mChi2Cumul.resize(mNbPoint);
2638  if(0 == mOptProfileIntegration.GetChoice())
2639  {
2641  if(mNbIntegrationUsed==0)
2642  mChi2Cumul=0;
2643  else
2644  {
2645  const REAL *pObs=mIntegratedObs.data();
2646  const REAL *pCalc=mPowderPatternIntegratedCalc.data();
2647  const REAL *pWeight;
2648  if(mIntegratedWeight.numElements()==0) pWeight=mIntegratedWeightObs.data();
2649  else pWeight=mIntegratedWeight.data();
2650 
2651  REAL *pC2Cu=mChi2Cumul.data();
2652  for(int i=0;i<mIntegratedPatternMin(0);i++) *pC2Cu++ = 0;
2653  REAL chi2cumul=0,tmp;
2654  for(unsigned long j=1;j<mNbIntegrationUsed;j++)
2655  {
2656  tmp=(*pObs++ - *pCalc++) ;
2657  chi2cumul += *pWeight++ * tmp*tmp;
2658  for(int i=mIntegratedPatternMin(j-1);i<mIntegratedPatternMin(j);i++) *pC2Cu++ =chi2cumul;
2659 
2660  if(mIntegratedPatternMin(j)>(int)mNbPointUsed)
2661  {
2662  for(unsigned int i=mIntegratedPatternMin(j);i<mNbPoint;i++) *pC2Cu++ =chi2cumul;
2663  break;
2664  }
2665  }
2666  pC2Cu=mChi2Cumul.data()+mIntegratedPatternMin(mNbIntegrationUsed-1);
2667  for(unsigned int i=mIntegratedPatternMin(mNbIntegrationUsed-1);i<mNbPoint;i++) *pC2Cu++ =chi2cumul;
2668  }
2669  }
2670  else
2671  {
2672  this->CalcPowderPattern();
2673  const REAL *pObs=mPowderPatternObs.data();
2674  const REAL *pCalc=mPowderPatternCalc.data();
2675  const REAL *pWeight=mPowderPatternWeight.data();
2676  REAL *pC2Cu=mChi2Cumul.data();
2677  REAL chi2cumul=0,tmp;
2678  for(unsigned int i=0;i<mNbPointUsed;i++)
2679  {
2680  VFN_DEBUG_MESSAGE("PowderPattern::GetChi2Cumul():"<<mIntegratedPatternMin(i)<<"->"<<mIntegratedPatternMax(i)<<":obs-calc="<<*pObs - *pCalc<<", weight="<<*pWeight,5);
2681  tmp = (*pObs++ - *pCalc++) ;
2682  chi2cumul += *pWeight++ * tmp*tmp;
2683  *pC2Cu++ = chi2cumul;
2684  }
2685  }
2686  VFN_DEBUG_EXIT("PowderPattern::GetChi2Cumul()",3)
2687  return mChi2Cumul;
2688 }
2689 
2691 { return mClockPowderPatternCalc;}
2692 
2694 { return mClockPowderPatternPar;}
2695 
2698 
2700 { return mClockPowderPatternXCorr;}
2701 
2702 void PowderPattern::SetXZero(const REAL newZero)
2703 {
2704  mXZero=newZero;
2706 }
2707 
2708 void PowderPattern::Set2ThetaDisplacement(const REAL displacement)
2709 {
2710  m2ThetaDisplacement=displacement;
2712 }
2713 
2714 void PowderPattern::Set2ThetaTransparency(const REAL transparency)
2715 {
2716  m2ThetaTransparency=transparency;
2718 }
2719 
2720 REAL PowderPattern::X2XCorr(const REAL x0)const
2721 {
2722  REAL x=x0;
2723  if( (mRadiation.GetWavelengthType()==WAVELENGTH_MONOCHROMATIC)
2724  ||(mRadiation.GetWavelengthType()==WAVELENGTH_ALPHA12))
2725  x += m2ThetaDisplacement*cos(x/2) +m2ThetaTransparency*sin(x);
2726 
2727  return x+mXZero;
2728 }
2729 
2730 REAL PowderPattern::X2PixelCorr(const REAL x0)const
2731 {
2732  return this->X2Pixel(this->X2XCorr(x0));
2733 }
2734 
2735 REAL PowderPattern::X2Pixel(const REAL x)const
2736 {
2737  //:TODO: faster if the step is actually constant.
2738  // Step may not be constant, so we guess twice before step-search
2739  REAL pixx;
2740  if(mIsXAscending==false)
2741  {
2742  VFN_DEBUG_MESSAGE("PowderPattern::X2Pixel()",1)
2743  long pix=(long)(mNbPoint-1-(x-this->GetPowderPatternXMin())/this->GetPowderPatternXStep());
2744  if((pix>0)&&(pix<((long)mNbPoint-1)))
2745  {
2746  // Why floor() and ceil() don't return a bloody integer is beyond me
2747  const REAL localStep=mX(pix)-mX(pix+1);
2748  if(localStep>0) pix -= (long)((x-mX(pix))/localStep);
2749  }
2750  VFN_DEBUG_MESSAGE("PowderPattern::X2Pixel():"<<x<<","<<pix,1)
2751  if(pix<1) pix=1;
2752  if(pix>((long)mNbPoint-2))pix=(long)mNbPoint-2;
2753  VFN_DEBUG_MESSAGE("PowderPattern::X2Pixel():"<<x<<","<<pix<<","<<mX(pix),1)
2754  if(mX(pix)<x)
2755  {
2756  for(;;pix--)
2757  {
2758  VFN_DEBUG_MESSAGE("PowderPattern::X2Pixel():"<<x<<","<<pix<<","<<mX(pix),1)
2759  if(mX(pix)>=x) break;
2760  if(pix==0) break;
2761  }
2762  }
2763  else
2764  {
2765  for(;;pix++)
2766  {
2767  if(mX(pix)<=x) {pix--;break;}
2768  if(pix==((long)mNbPoint-2)) break;
2769  }
2770  }
2771  // This assumes step is at least localy constant...
2772  VFN_DEBUG_MESSAGE("PowderPattern::X2Pixel():"<<x<<","<<pix<<","<<mX(pix),1)
2773  const REAL localStep=mX(pix)-mX(pix+1);
2774  pixx = (REAL)pix-(x-mX(pix))/localStep;
2775  VFN_DEBUG_MESSAGE("PowderPattern::X2Pixel():"<<x<<","<<pix<<","<<mX(pix),1)
2776  }
2777  else
2778  {
2779  VFN_DEBUG_MESSAGE("PowderPattern::X2Pixel():"<<x<<","<<this->GetPowderPatternXMin()<<","<<this->GetPowderPatternXMax(),1)
2780  long pix=(long)((x-this->GetPowderPatternXMin())/this->GetPowderPatternXStep());
2781  if((pix>0)&&(pix<((long)mNbPoint-1)))
2782  {
2783  // Why floor() and ceil() don't return a bloody integer is beyond me
2784  const REAL localStep=mX(pix+1)-mX(pix);
2785  if(localStep>0) pix += (long)((x-mX(pix))/localStep);
2786  }
2787  VFN_DEBUG_MESSAGE("PowderPattern::X2Pixel():"<<x<<","<<pix,1)
2788  if(pix<1) pix=1;
2789  if(pix>((long)mNbPoint-2))pix=(long)mNbPoint-2;
2790  VFN_DEBUG_MESSAGE("PowderPattern::X2Pixel():"<<x<<","<<pix<<","<<mX(pix),1)
2791  if(x<mX(pix))
2792  {
2793  for(;;pix--)
2794  {
2795  VFN_DEBUG_MESSAGE("PowderPattern::X2Pixel():"<<x<<","<<pix<<","<<mX(pix),1)
2796  if(mX(pix)<=x) break;
2797  if(pix==0) break;
2798  }
2799  }
2800  else
2801  {
2802  for(;;pix++)
2803  {
2804  VFN_DEBUG_MESSAGE("PowderPattern::X2Pixel():"<<x<<","<<pix<<","<<mX(pix),1)
2805  if(mX(pix)>=x) {pix-- ;break;}
2806  if(pix==((long)mNbPoint-2)) break;
2807  }
2808  }
2809  VFN_DEBUG_MESSAGE("PowderPattern::X2Pixel():"<<x<<","<<pix<<","<<mX(pix),1)
2810  if(pix>((long)mNbPoint-2))pix=(long)mNbPoint-2;
2811  // This assumes step is at least localy constant...
2812  const REAL localStep=mX(pix+1)-mX(pix);
2813  VFN_DEBUG_MESSAGE("PowderPattern::X2Pixel():"<<x<<","<<pix<<","<<mX(pix)<<","<<localStep,1)
2814  pixx = (REAL)pix+(x-mX(pix))/localStep;
2815  }
2816  VFN_DEBUG_MESSAGE("PowderPattern::X2Pixel():"<<x<<","<<pixx,1)
2817  return pixx;
2818 }
2819 
2821 {
2822  //15.000 0.030 70.000 LANI4FE#1 REC 800 4JRS
2823  //2447. 2418. 2384. 2457. 2398. 2374. 2378. 2383.
2824  //...
2825  VFN_DEBUG_MESSAGE("PowderPattern::ImportPowderPatternFullprof() : \
2826 from file : "+filename,5)
2827  ifstream fin(filename.c_str());
2828  if(!fin)
2829  {
2830  throw ObjCrystException("PowderPattern::ImportPowderPatternFullprof() : \
2831 Error opening file for input:"+filename);
2832  }
2833  REAL min,max,step;
2834  fin >> min >> step >> max;
2835  min *= DEG2RAD;
2836  max *= DEG2RAD;
2837  step *= DEG2RAD;
2838  this->SetPowderPatternPar(min,step,(long)((max-min)/step+1.001));
2839  VFN_DEBUG_MESSAGE("PowderPattern::ImportPowderPatternFullprof() :"\
2840  << " 2Theta min=" << min*RAD2DEG << " 2Theta max=" << max*RAD2DEG \
2841  << " NbPoints=" << mNbPoint,5)
2842  mPowderPatternObs.resize (mNbPoint);
2845 
2846  char tmpComment[200];
2847  fin.getline(tmpComment,100);
2848  //if(""==mName) mName.append(tmpComment);
2849 
2850  for(unsigned long i=0;i<mNbPoint;i++) fin >> mPowderPatternObs(i);
2851  fin.close();
2852  this->SetSigmaToSqrtIobs();
2853  this->SetWeightToInvSigmaSq();
2855  this->UpdateDisplay();
2856  {
2857  char buf [200];
2858  sprintf(buf,"Imported powder pattern: %d points, 2theta=%7.3f -> %7.3f, step=%6.3f",
2859  (int)mNbPoint,min*RAD2DEG,max*RAD2DEG,step*RAD2DEG);
2860  (*fpObjCrystInformUser)((string)buf);
2861  }
2862  VFN_DEBUG_MESSAGE("PowderPattern::ImportFullProfPattern():finished:"<<mNbPoint<<" points",5)
2863 }
2864 
2865 void PowderPattern::ImportPowderPatternPSI_DMC(const string &filename)
2866 {
2867  VFN_DEBUG_MESSAGE("PowderPattern::ImportPowderPatternPSI_DMC() : \
2868 from file : "+filename,5)
2869  ifstream fin(filename.c_str());
2870  if(!fin)
2871  {
2872  throw ObjCrystException("PowderPattern::ImportPowderPatternPSI_DMC() : \
2873 Error opening file for input:"+filename);
2874  }
2875  //Skip the first two lines
2876  char tmpComment[200];
2877  fin.getline(tmpComment,190);
2878  fin.getline(tmpComment,190);
2879  REAL min,max,step;
2880  fin >> min >> step >> max;
2881  min *= DEG2RAD;
2882  max *= DEG2RAD;
2883  step *= DEG2RAD;
2884  this->SetPowderPatternPar(min,step,(unsigned long)((max-min)/step+1.001));
2885  VFN_DEBUG_MESSAGE("PowderPattern::ImportPowderPatternPSI_DMC() :"\
2886  << " 2Theta min=" << min*RAD2DEG << " 2Theta max=" << max*RAD2DEG \
2887  << " NbPoints=" << mNbPoint,5)
2888  mPowderPatternObs.resize (mNbPoint);
2891 
2892  fin.getline(tmpComment,100);
2893  //if(""==mName) mName.append(tmpComment);
2894 
2895  for(unsigned long i=0;i<mNbPoint;i++) fin >> mPowderPatternObs(i);
2896  for(unsigned long i=0;i<mNbPoint;i++) fin >> mPowderPatternObsSigma(i);
2897  fin.close();
2898  this->SetWeightToInvSigmaSq();
2900  this->UpdateDisplay();
2901  {
2902  char buf [200];
2903  sprintf(buf,"Imported powder pattern: %d points, 2theta=%7.3f -> %7.3f, step=%6.3f",
2904  (int)mNbPoint,min*RAD2DEG,max*RAD2DEG,step*RAD2DEG);
2905  (*fpObjCrystInformUser)((string)buf);
2906  }
2907  VFN_DEBUG_MESSAGE("PowderPattern::ImportPowderPatternPSI_DMC():finished",5)
2908 }
2909 
2911 {
2912  VFN_DEBUG_MESSAGE("PowderPattern::ImportPowderPatternILL_D1AD2B() : \
2913 from file : "+filename,5)
2914  ifstream fin(filename.c_str());
2915  if(!fin)
2916  {
2917  throw ObjCrystException("PowderPattern::ImportPowderPatternILL_D1AD2B() : \
2918 Error opening file for input:"+filename);
2919  }
2920  //Skip the first three lines
2921  char tmpComment[200];
2922  fin.getline(tmpComment,190);
2923  fin.getline(tmpComment,190);
2924  fin.getline(tmpComment,190);
2925 
2926  fin >> mNbPoint;
2927  fin.getline(tmpComment,190);
2928  REAL min,step;
2929  fin >> min >> step;
2930  min *= DEG2RAD;
2931  step *= DEG2RAD;
2932  this->SetPowderPatternPar(min,step,mNbPoint);
2933  VFN_DEBUG_MESSAGE("PowderPattern::ImportPowderPatternILL_D1AD2B() :"\
2934  << " 2Theta min=" << min*RAD2DEG << " 2Theta max=" << min*RAD2DEG+mNbPoint*step*RAD2DEG \
2935  << " NbPoints=" << mNbPoint,5)
2936  mPowderPatternObs.resize (mNbPoint);
2937  mPowderPatternObsSigma.resize (mNbPoint);
2938  mPowderPatternWeight.resize(mNbPoint);
2939 
2940  //if(""==mName) mName.append(tmpComment);
2941 
2942  for(unsigned long i=0;i<mNbPoint;i++) fin >> mPowderPatternObs(i);
2943  for(unsigned long i=0;i<mNbPoint;i++) fin >> mPowderPatternObsSigma(i);
2944  fin.close();
2945  this->SetWeightToInvSigmaSq();
2947  this->UpdateDisplay();
2948  {
2949  char buf [200];
2950  sprintf(buf,"Imported powder pattern: %d points, 2theta=%7.3f -> %7.3f, step=%6.3f",
2951  (int)mNbPoint,min*RAD2DEG,(min+step*(mNbPoint-1))*RAD2DEG,
2952  step*RAD2DEG);
2953  (*fpObjCrystInformUser)((string)buf);
2954  }
2955  VFN_DEBUG_MESSAGE("PowderPattern::ImportPowderPatternILL_D1AD2B():finished",5)
2956 }
2957 
2958 void PowderPattern::ImportPowderPatternXdd(const string &filename)
2959 {
2960  VFN_DEBUG_MESSAGE("PowderPattern::ImportPowderPatternXdd():from file :" \
2961  +filename,5)
2962  ifstream fin (filename.c_str());
2963  if(!fin)
2964  {
2965  throw ObjCrystException("PowderPattern::ImportPowderPatternXdd() : \
2966 Error opening file for input:"+filename);
2967  }
2968  char tmpComment[200];
2969  fin.getline(tmpComment,100);
2970  //if(""==mName) mName.append(tmpComment);
2971  REAL min,max,step,tmp;
2972  fin >> min >> step >> max;
2973  min *= DEG2RAD;
2974  max *= DEG2RAD;
2975  step *= DEG2RAD;
2976  this->SetPowderPatternPar(min,step,(long)((max-min)/step+1.001));
2977  mPowderPatternObs.resize (mNbPoint);
2979  mPowderPatternCalc.resize(mNbPoint);
2982 
2983  fin >> tmp; //Count time
2984  fin >> tmp; //unused
2985  fin >> tmp; //unused (wavelength?)
2986 
2987  for(unsigned long i=0;i<mNbPoint;i++) fin >> mPowderPatternObs(i);
2988  fin.close();
2989  this->SetSigmaToSqrtIobs();
2990  this->SetWeightToInvSigmaSq();
2991  this->UpdateDisplay();
2992  {
2993  char buf [200];
2994  sprintf(buf,"Imported powder pattern: %d points, 2theta=%7.3f -> %7.3f, step=%6.3f",
2995  (int)mNbPoint,min*RAD2DEG,max*RAD2DEG,step*RAD2DEG);
2996  (*fpObjCrystInformUser)((string)buf);
2997  }
2998  VFN_DEBUG_MESSAGE("DiffractionDataPowder::ImportXddPattern() :finished",5)
2999 }
3000 
3002 {
3003  VFN_DEBUG_ENTRY("PowderPattern::ImportPowderPatternSietronicsCPI():from file :" \
3004  +filename,5)
3005  ifstream fin (filename.c_str());
3006  if(!fin)
3007  {
3008  throw ObjCrystException("PowderPattern::ImportPowderPatternSietronicsCPI() : \
3009 Error opening file for input:"+filename);
3010  }
3011  char tmpComment[300];
3012  fin.getline(tmpComment,100);
3013  VFN_DEBUG_MESSAGE(" ->Discarded comment :"<<tmpComment,1)
3014  REAL min,max,step;
3015  fin >> min >> max >> step;
3016  min *= DEG2RAD;
3017  max *= DEG2RAD;
3018  step *= DEG2RAD;
3019  this->SetPowderPatternPar(min,step,(long)((max-min)/step+1.001));
3020  mPowderPatternObs.resize (mNbPoint);
3022  mPowderPatternCalc.resize(mNbPoint);
3025 
3026  //Following lines are ignored (no fixed format ?)
3027  string str;
3028  do
3029  {
3030  fin>>str;
3031  VFN_DEBUG_MESSAGE(" ->Read :"<<str,1)
3032  } while ("SCANDATA"!=str);
3033 
3034  for(unsigned long i=0;i<mNbPoint;i++) fin >> mPowderPatternObs(i);
3035  fin.close();
3036  this->SetSigmaToSqrtIobs();
3037  this->SetWeightToInvSigmaSq();
3039  this->UpdateDisplay();
3040  {
3041  char buf [200];
3042  sprintf(buf,"Imported powder pattern: %d points, 2theta=%7.3f -> %7.3f, step=%6.3f",
3043  (int)mNbPoint,min*RAD2DEG,max*RAD2DEG,step*RAD2DEG);
3044  (*fpObjCrystInformUser)((string)buf);
3045  }
3046  VFN_DEBUG_EXIT("DiffractionDataPowder::ImportPowderPatternSietronicsCPI()",5)
3047 }
3048 
3049 void PowderPattern::ImportPowderPattern2ThetaObsSigma(const string &filename,const int nbSkip)
3050 {
3051  VFN_DEBUG_MESSAGE("DiffractionDataPowder::ImportPowderPattern2ThetaObsSigma():from:" \
3052  +filename,5)
3053  ifstream fin (filename.c_str());
3054  if(!fin)
3055  {
3056  throw ObjCrystException("PowderPattern::ImportPowderPattern2ThetaObsSigma():\
3057 Error opening file for input:"+filename);
3058  }
3059  {//Get rid of first lines
3060  char tmpComment[200];
3061  for(int i=0;i<nbSkip;i++) fin.getline(tmpComment,150);
3062  }
3063  mPowderPatternObs.resize (500);
3064  mPowderPatternObsSigma.resize(500);
3065  mX.resize(500);
3066  mNbPoint=0;
3067 
3068  do
3069  {
3070  fin >> mX (mNbPoint);
3071  fin >> mPowderPatternObs (mNbPoint);
3073  //cout << mX (mNbPoint)<<" "
3074  // << mPowderPatternObs (mNbPoint)<<" "
3075  // << mPowderPatternObsSigma(mNbPoint)<<endl;
3076  if(!fin) break;
3077  mNbPoint++;
3078  if( (mNbPoint%500)==0)
3079  {
3080  mX.resizeAndPreserve(mNbPoint+500);
3081  mPowderPatternObs.resizeAndPreserve(mNbPoint+500);
3082  mPowderPatternObsSigma.resizeAndPreserve(mNbPoint+500);
3083  }
3084  } while(fin.eof() == false);
3085  fin.close();
3086 
3087  mX.resizeAndPreserve (mNbPoint);
3088  mPowderPatternObs.resizeAndPreserve (mNbPoint);
3089  mPowderPatternObsSigma.resizeAndPreserve(mNbPoint);
3091 
3092  mX *= DEG2RAD;
3093 
3094  this->SetWeightToInvSigmaSq();
3096  this->UpdateDisplay();
3097  {
3098  char buf [200];
3099  sprintf(buf,"Imported powder pattern: %d points, 2theta=%7.3f -> %7.3f, step=%6.3f",
3100  (int)mNbPoint,this->GetPowderPatternXMin()*RAD2DEG,
3101  this->GetPowderPatternXMax()*RAD2DEG,
3102  this->GetPowderPatternXStep()*RAD2DEG);
3103  (*fpObjCrystInformUser)((string)buf);
3104  }
3105  VFN_DEBUG_MESSAGE("PowderPattern::ImportPowderPattern2ThetaObsSigma()\
3106 :finished: "<<mNbPoint<<" points",5)
3107 }
3108 
3109 void PowderPattern::ImportPowderPattern2ThetaObs(const string &filename,const int nbSkip)
3110 {
3111  VFN_DEBUG_MESSAGE("PowderPattern::ImportPowderPattern2ThetaObs():from:" \
3112  +filename,5)
3113  ifstream fin (filename.c_str());
3114  if(!fin)
3115  {
3116  throw ObjCrystException("PowderPattern::ImportPowderPattern2ThetaObs():\
3117 Error opening file for input:"+filename);
3118  }
3119  {//Get rid of first lines
3120  char tmpComment[200];
3121  for(int i=0;i<nbSkip;i++) fin.getline(tmpComment,150);
3122  }
3123  mPowderPatternObs.resize (500);
3124  mPowderPatternObsSigma.resize(500);
3125  mX.resize(500);
3126  mNbPoint=0;
3127 
3128  do
3129  {
3130  fin >> mX(mNbPoint);
3131  fin >> mPowderPatternObs (mNbPoint);
3133  =sqrt(mPowderPatternObs(mNbPoint));
3134  if(!fin) break;
3135  mNbPoint++;
3136  if( (mNbPoint%500)==0)
3137  {
3138  mX.resizeAndPreserve(mNbPoint+500);
3139  mPowderPatternObs.resizeAndPreserve(mNbPoint+500);
3140  mPowderPatternObsSigma.resizeAndPreserve(mNbPoint+500);
3141  }
3142  } while(fin.eof() == false);
3143  fin.close();
3144 
3145  mX.resizeAndPreserve (mNbPoint);
3146  mPowderPatternObs.resizeAndPreserve (mNbPoint);
3147  mPowderPatternObsSigma.resizeAndPreserve(mNbPoint);
3149 
3150  mX *= DEG2RAD;
3151 
3152  this->SetSigmaToSqrtIobs();
3153  this->SetWeightToInvSigmaSq();
3155  this->UpdateDisplay();
3156  {
3157  char buf [200];
3158  sprintf(buf,"Imported powder pattern: %d points, 2theta=%7.3f -> %7.3f, step=%6.3f",
3159  (int)mNbPoint,this->GetPowderPatternXMin()*RAD2DEG,
3160  this->GetPowderPatternXMax()*RAD2DEG,
3161  this->GetPowderPatternXStep()*RAD2DEG);
3162  (*fpObjCrystInformUser)((string)buf);
3163  }
3164  VFN_DEBUG_MESSAGE("PowderPattern::ImportPowderPattern2ThetaObs():finished",5)
3165 }
3166 
3168 {
3169 
3170  //Sample 4: NaY + CF2=CCL2 T=20K, Lambda: 2.343 A.
3171  // 100 0 0.100 70 0 0
3172  // 3.000
3173  // 500000. 12000. 0.00 0.00
3174  //70 570369 569668 562868 532769 527469 495669 481969 452767 429468 4132
3175  //68 393269 372067 353769 337068 328268 310270 299469 296470 294668 2780
3176  //...
3177  //14 255814 282714 274014 281314 300314 302714 301114 298014 313214 3097
3178  //14 286914 295714 305214 300714 288311 295511 288511 3024 8 2937 7 2883
3179  // 7 2905 7 2895 7 2767 7 2777 7 2758 7 2495 7 2507 7 2496 7 2382 7 2329
3180  // 7 2542 7 2415 7 2049 7 2389 7 2270 7 2157 6 2227 6 2084 3 1875 2 2094
3181  // 1 1867
3182  // -1000
3183  // -10000
3184  VFN_DEBUG_MESSAGE("PowderPattern::ImportPowderPatternMultiDetectorLLBG42() : \
3185 from file : "+filename,5)
3186  ifstream fin(filename.c_str());
3187  if(!fin)
3188  {
3189  throw ObjCrystException("PowderPattern::ImportPowderPatternMultiDetectorLLBG42() : \
3190 Error opening file for input:"+filename);
3191  }
3192 
3193  string str;
3194  getline(fin,str);
3195  float junk;
3196  REAL min,step;
3197  fin >>junk>>junk>>step>>junk>>junk>>junk>>min>>junk>>junk>>junk>>junk;
3198  min *= DEG2RAD;
3199  step *= DEG2RAD;
3200  VFN_DEBUG_MESSAGE("PowderPattern::ImportPowderPatternMultiDetectorLLBG42() :"\
3201  << " 2Theta min=" << min*RAD2DEG << " 2Theta step=" << step*RAD2DEG,5)
3202  mPowderPatternObs.resize (500);
3203  mPowderPatternObsSigma.resize (500);
3204 
3205  getline(fin,str);//finish reading line
3206 
3207  float tmp;
3208  string sub;
3209  float ct,iobs;
3210  mNbPoint=0;
3211  for(;;)
3212  {
3213  getline(fin,str);
3214  sscanf(str.c_str(),"%f",&tmp);
3215  if(tmp<0) break;
3216  const unsigned int nb=str.length()/8;
3217  for(unsigned int i=0;i<nb;i++)
3218  {
3219  if(mNbPoint==(unsigned int)mPowderPatternObs.numElements())
3220  {
3221  mPowderPatternObs.resizeAndPreserve(mNbPoint+500);
3222  mPowderPatternObsSigma.resizeAndPreserve(mNbPoint+500);
3223  }
3224  sub=str.substr(i*8,8);
3225  sscanf(sub.c_str(),"%2f%6f",&ct,&iobs);
3227  mPowderPatternObsSigma(mNbPoint++)=sqrt(iobs/ct);
3228  }
3229  }
3230  this->SetPowderPatternPar(min,step,mNbPoint);
3231  //exit(1);
3232  mPowderPatternObs.resizeAndPreserve (mNbPoint);
3233  mPowderPatternObsSigma.resizeAndPreserve (mNbPoint);
3234  mPowderPatternWeight.resizeAndPreserve(mNbPoint);
3235 
3236  fin.close();
3237  this->SetWeightToInvSigmaSq();
3239  this->UpdateDisplay();
3240  {
3241  char buf [200];
3242  sprintf(buf,"Imported powder pattern: %d points, 2theta=%7.3f -> %7.3f, step=%6.3f",
3243  (int)mNbPoint,min*RAD2DEG,mX(mNbPoint-1)*RAD2DEG,step*RAD2DEG);
3244  (*fpObjCrystInformUser)((string)buf);
3245  }
3246  VFN_DEBUG_MESSAGE("PowderPattern::ImportPowderPatternMultiDetectorLLBG42():finished:"<<mNbPoint<<" points",5)
3247 }
3248 
3250 {
3251  //1.550 0.005 66.000
3252  // 213.135 193.243 208.811 185.873 231.607 200.995 196.792 187.516 215.977 199.634
3253  // 17.402 16.570 12.180 11.491 18.141 11.950 11.824 11.542 17.518 11.909
3254  // 211.890 185.740 204.610 200.645 199.489 169.549 203.189 178.298 186.241 198.522
3255  // 12.269 11.487 17.051 11.939 11.905 10.975 16.992 11.255 11.503 11.876
3256  VFN_DEBUG_MESSAGE("PowderPattern::ImportPowderPatternFullprof4() : \
3257 from file : "+filename,5)
3258  ifstream fin(filename.c_str());
3259  if(!fin)
3260  {
3261  throw ObjCrystException("PowderPattern::ImportPowderPatternFullprof4() : \
3262 Error opening file for input:"+filename);
3263  }
3264  REAL min,step,max;
3265  fin >> min >> step >> max;
3266  min *= DEG2RAD;
3267  max *= DEG2RAD;
3268  step *= DEG2RAD;
3269  this->SetPowderPatternPar(min,step,(long)((max-min)/step+1.001));
3270  VFN_DEBUG_MESSAGE("PowderPattern::ImportPowderPatternFullprof4() :"\
3271  << " 2Theta min=" << min*RAD2DEG << " 2Theta max=" << max*RAD2DEG \
3272  << " NbPoints=" << mNbPoint,5)
3273  mPowderPatternObs.resize (mNbPoint);
3276 
3277  string str;
3278  getline(fin,str);//read end of first line
3279 
3280  unsigned ct=0;
3281  unsigned ctSig=0;
3282  float line[10];
3283  for(;ct<mNbPoint;)
3284  {
3285  getline(fin,str);
3286  for(unsigned int i=0;i<str.size();i++) if(' '==str[i]) str[i]='0';
3287  sscanf(str.c_str(),"%8f%8f%8f%8f%8f%8f%8f%8f%8f%8f",
3288  line+0,line+1,line+2,line+3,line+4,line+5,line+6,line+7,line+8,line+9);
3289  for(unsigned int j=0;j<10;j++)
3290  if(ct<mNbPoint) mPowderPatternObs(ct++)=line[j];
3291  getline(fin,str);
3292  for(unsigned int i=0;i<str.size();i++) if(' '==str[i]) str[i]='0';
3293  sscanf(str.c_str(),"%8f%8f%8f%8f%8f%8f%8f%8f%8f%8f",
3294  line+0,line+1,line+2,line+3,line+4,line+5,line+6,line+7,line+8,line+9);
3295  for(unsigned int j=0;j<10;j++)
3296  if(ctSig<mNbPoint) mPowderPatternObsSigma(ctSig++)=line[j];
3297  }
3298  fin.close();
3299  this->SetWeightToInvSigmaSq();
3301  this->UpdateDisplay();
3302  {
3303  char buf [200];
3304  sprintf(buf,"Imported powder pattern: %d points, 2theta=%7.3f -> %7.3f, step=%6.3f",
3305  (int)mNbPoint,min*RAD2DEG,max*RAD2DEG,
3306  step*RAD2DEG);
3307  (*fpObjCrystInformUser)((string)buf);
3308  }
3309  VFN_DEBUG_MESSAGE("PowderPattern::ImportFullProfPattern4():finished:"<<mNbPoint<<" points",5)
3310 }
3311 
3313 {
3314  VFN_DEBUG_MESSAGE("DiffractionDataPowder::ImportPowderPatternTOF_ISIS_XYSigma():from:" \
3315  +filename,5)
3316  ifstream fin (filename.c_str());
3317  if(!fin)
3318  {
3319  throw ObjCrystException("PowderPattern::ImportPowderPatternTOF_ISIS_XYSigma():\
3320 Error opening file for input:"+filename);
3321  }
3322  {//Get rid of first line
3323  char junk[400];
3324  fin.getline(junk,150);
3325  }
3326  mPowderPatternObs.resize (500);
3327  mPowderPatternObsSigma.resize(500);
3328  mX.resize(500);
3329  mNbPoint=0;
3330 
3331  do
3332  {
3333  fin >> mX (mNbPoint);
3334  fin >> mPowderPatternObs (mNbPoint);
3336  //cout << mX (mNbPoint)<<" "
3337  // << mPowderPatternObs (mNbPoint)<<" "
3338  // << mPowderPatternObsSigma(mNbPoint)<<endl;
3339  if(!fin) break;
3340  mNbPoint++;
3341  if( (mNbPoint%500)==0)
3342  {
3343  mX.resizeAndPreserve(mNbPoint+500);
3344  mPowderPatternObs.resizeAndPreserve(mNbPoint+500);
3345  mPowderPatternObsSigma.resizeAndPreserve(mNbPoint+500);
3346  }
3347  } while(fin.eof() == false);
3348  fin.close();
3349 
3350  mX.resizeAndPreserve (mNbPoint);
3351  mPowderPatternObs.resizeAndPreserve (mNbPoint);
3352  mPowderPatternObsSigma.resizeAndPreserve(mNbPoint);
3354 
3355  // Reverse order of arrays, so that we are in ascending order of sin(theta)/lambda
3356  REAL tmp;
3357  for(unsigned long i=0;i<(mNbPoint/2);i++)
3358  {
3359  tmp=mX(i);
3360  mX(i)=mX(mNbPoint-1-i);
3361  mX(mNbPoint-1-i)=tmp;
3362 
3363  tmp=mPowderPatternObs(i);
3365  mPowderPatternObs(mNbPoint-1-i)=tmp;
3366 
3367  tmp=mPowderPatternObsSigma(i);
3370  }
3371  this->SetPowderPatternX(mX);
3372 
3373  this->SetWeightToInvSigmaSq();
3374  this->SetRadiationType(RAD_NEUTRON);
3375  this->GetRadiation().SetWavelengthType(WAVELENGTH_TOF);
3377  this->UpdateDisplay();
3378  {
3379  char buf [200];
3380  sprintf(buf,"Imported TOF powder pattern: %d points, TOF=%7.3f -> %7.3f",
3381  (int)mNbPoint,this->GetPowderPatternXMin(),
3382  this->GetPowderPatternXMax());
3383  (*fpObjCrystInformUser)((string)buf);
3384  }
3385  VFN_DEBUG_MESSAGE("PowderPattern::ImportPowderPatternTOF_ISIS_XYSigma()\
3386 :finished: "<<mNbPoint<<" points",5)
3387 }
3388 
3389 void PowderPattern::ImportPowderPatternGSAS(const string &filename)
3390 {
3391  VFN_DEBUG_ENTRY("PowderPattern::ImportPowderPatternGSAS():file:"<<filename,5)
3392  ifstream fin (filename.c_str());
3393  if(!fin)
3394  {
3395  throw ObjCrystException("PowderPattern::ImportPowderPatternGSAS():\
3396 Error opening file for input:"+filename);
3397  }
3398  {//Get rid of title
3399  char title[81];
3400  fin.read(title,80);
3401  while(isprint(fin.peek())==false)
3402  {
3403  if(fin.eof()) break;
3404  fin.get();
3405  }
3406  cout<<"Title:"<<title<<endl;
3407  if(this->GetName()=="Change Me!") this->SetName(title);
3408  }
3409 
3410  //BANK 1 38101 7620 CONST 200.00 0.10 0 0 ESD
3411  int numBank,nbRecords;
3412  string binType, type;
3413  float bcoeff[4];
3414  //string line;
3415  char line[81];
3416  char bank[5];
3417  do
3418  {
3419  fin.getline(line,80);
3420  while(isprint(fin.peek())==false)
3421  {
3422  if(fin.eof()) break;
3423  fin.get();
3424  }
3425  sscanf(line,"%4s",bank);
3426  if(fin.eof())
3427  throw ObjCrystException("PowderPattern::ImportPowderPatternGSAS():\
3428 Could not find BANK statement !! In file: "+filename);
3429  }
3430  while(string(bank)!=string("BANK"));
3431 
3432  {
3433  line[80]='\0';
3434  char binTypeC[20],typeC[20];
3435  sscanf(line,"%4s%d %ld %d %s %f %f %f %f %s",bank,&numBank,&mNbPoint,&nbRecords,
3436  binTypeC,&bcoeff[0],&bcoeff[1],&bcoeff[2],&bcoeff[3],typeC);
3437  binType=binTypeC;
3438  type=typeC;
3439  }
3440  if(binType=="CONST") binType="CONS";
3441  if((type!="ALT")&&(type!="ESD")) type="STD";
3442 
3443  cout<<"BANK #"<<numBank<<endl;
3444  cout<<"Number of data points:"<<mNbPoint<<endl;
3445  cout<<"Number of records:"<<nbRecords<<endl;
3446  cout<<"BinType:"<<binType<<endl;
3447  cout<<"BCoeff[1-4]:"<<bcoeff[0]<<","<<bcoeff[1]<<","<<bcoeff[2]<<","<<bcoeff[3]<<endl;
3448  cout<<"Type:"<<type<<endl;
3449 
3450  mPowderPatternObs.resize (mNbPoint);
3452  mX.resize(mNbPoint);
3453  bool importOK=false;
3454  if((binType=="CONS") && (type=="ESD"))
3455  {
3456  this->SetPowderPatternPar(bcoeff[0]*DEG2RAD/100,bcoeff[1]*DEG2RAD/100,mNbPoint);
3457  string sub;
3458  unsigned long point=0;
3459  REAL iobs,isig;
3460  string substr;
3461  for(long i=0;i<nbRecords;i++)
3462  {
3463  fin.read(line,80);
3464  line[80]='\0';
3465  while(isprint(fin.peek())==false)
3466  {
3467  if(fin.eof()) break;
3468  fin.get();
3469  }
3470  for(unsigned int j=0;j<5;j++)
3471  {
3472  /*
3473  substr=string(line).substr(j*16,16);
3474  sscanf(substr.c_str(),"%8f%8f",&iobs,&isig);
3475  */
3476  substr=string(line).substr(j*16+0 ,8);
3477  istringstream(substr) >> iobs;
3478  substr=string(line).substr(j*16+8 ,8);
3479  istringstream(substr) >> isig;
3480 
3481  mPowderPatternObs(point)=iobs;
3482  mPowderPatternObsSigma(point++)=isig;
3483  if(point==mNbPoint) break;
3484  }
3485  if(point==mNbPoint) break;
3486  }
3487  importOK=true;
3488  }
3489  if((binType=="CONS") && (type=="STD"))
3490  {
3491  this->SetPowderPatternPar(bcoeff[0]*DEG2RAD/100,bcoeff[1]*DEG2RAD/100,mNbPoint);
3492  unsigned long point=0;
3493  REAL iobs;
3494  int nc;
3495  string substr;
3496  for(long i=0;i<nbRecords;i++)
3497  {
3498  fin.read(line,80);
3499  line[80]='\0';
3500  while(isprint(fin.peek())==false)
3501  {
3502  if(fin.eof()) break;
3503  fin.get();
3504  }
3505  for(unsigned int j=0;j<10;j++)
3506  {
3507  /*
3508  substr=string(line).substr(j*8,8);
3509  if(substr.substr(0,2)==string(" "))
3510  {
3511  nc=1;
3512  sscanf(substr.c_str(),"%8f",&iobs);
3513  }
3514  else sscanf(substr.c_str(),"%2d%6f",&nc,&iobs);
3515  */
3516  substr=string(line).substr(j*8+0 ,2);
3517  if(substr==" ") nc=1;
3518  else sscanf(substr.c_str(),"%d",&nc);
3519  substr=string(line).substr(j*8+2 ,6);
3520  istringstream(substr) >> iobs;
3521  mPowderPatternObs(point)=iobs;
3522  mPowderPatternObsSigma(point++)=sqrt(iobs)/sqrt((REAL)nc);
3523  if(point==mNbPoint) break;
3524  }
3525  if(point==mNbPoint) break;
3526  }
3527  importOK=true;
3528  }
3529  if((binType=="RALF") && (type=="ALT"))
3530  {
3531  this->SetRadiationType(RAD_NEUTRON);
3532  this->GetRadiation().SetWavelengthType(WAVELENGTH_TOF);
3534 
3535  unsigned long point=0;
3536  REAL x,iobs,iobssigma;
3537  string substr;
3538  for(long i=0;i<nbRecords;i++)
3539  {
3540  fin.read(line,80);
3541  line[80]='\0';
3542  while(isprint(fin.peek())==false)
3543  {
3544  if(fin.eof()) break;
3545  fin.get();
3546  }
3547  for(unsigned int j=0;j<4;j++)
3548  {//4 records per line
3549  /* Does not work because sscanf ignores the leading spaces and shifts the reading !
3550  substr=string(line).substr(j*20,20);
3551  sscanf(substr.c_str(),"%8f%7f%5f",&x,&iobs,&iobssigma);
3552  */
3553  substr=string(line).substr(j*20+0 ,8);
3554  istringstream(substr) >> x;
3555  substr=string(line).substr(j*20+8 ,7);
3556  istringstream(substr) >> iobs;
3557  substr=string(line).substr(j*20+15,5);
3558  istringstream(substr) >> iobssigma;
3559  mPowderPatternObs(point)=iobs;
3560  mPowderPatternObsSigma(point)=iobssigma;
3561  mX(point)=x/32;
3562  if(++point==mNbPoint) break;
3563  }
3564  if(point==mNbPoint) break;
3565  }
3566  // Reverse order of arrays, so that we are in ascending order of sin(theta)/lambda
3567  REAL tmp;
3568  for(unsigned long i=0;i<(mNbPoint/2);i++)
3569  {
3570  tmp=mX(i);
3571  mX(i)=mX(mNbPoint-1-i);
3572  mX(mNbPoint-1-i)=tmp;
3573 
3574  tmp=mPowderPatternObs(i);
3576  mPowderPatternObs(mNbPoint-1-i)=tmp;
3577 
3578  tmp=mPowderPatternObsSigma(i);
3581  }
3582  importOK=true;
3583  }
3584  fin.close();
3585  if(!importOK)
3586  {
3587  mNbPoint=0;
3588  mPowderPatternObs.resize (mNbPoint);
3590  mX.resize(mNbPoint);
3591  throw ObjCrystException("PowderPattern::ImportPowderPatternGSAS(): Sorry, \
3592 this type of format is not handled yet (send an example file to the Fox author)!:"+filename);
3593  }
3595  this->SetPowderPatternX(mX);
3596  this->SetWeightToInvSigmaSq();
3597 
3598  this->UpdateDisplay();
3599  if(this->GetRadiation().GetWavelengthType()==WAVELENGTH_TOF)
3600  {
3601  char buf [200];
3602  sprintf(buf,"Imported powder pattern: %d points, tof=%7.3f us-> %7.3f us",
3603  (int)mNbPoint,this->GetPowderPatternXMin(),
3604  this->GetPowderPatternXMax());
3605  (*fpObjCrystInformUser)((string)buf);
3606  }
3607  else
3608  {
3609  char buf [200];
3610  sprintf(buf,"Imported powder pattern: %d points, 2theta=%7.3f -> %7.3f, step=%6.3f",
3611  (int)mNbPoint,this->GetPowderPatternXMin()*RAD2DEG,
3612  this->GetPowderPatternXMax()*RAD2DEG,
3613  this->GetPowderPatternXStep()*RAD2DEG);
3614  (*fpObjCrystInformUser)((string)buf);
3615  }
3616  VFN_DEBUG_EXIT("PowderPattern::ImportPowderPatternGSAS():file:"<<filename,5)
3617 }
3618 
3620 {
3621  VFN_DEBUG_ENTRY("PowderPattern::ImportPowderPatternCIF():file:",5)
3622  for(map<string,CIFData>::const_iterator pos=cif.mvData.begin();pos!=cif.mvData.end();++pos)
3623  if(pos->second.mPowderPatternObs.size()>10)
3624  {
3625  mNbPoint=pos->second.mPowderPatternObs.size();
3626  mX.resize(mNbPoint);
3627  mPowderPatternObs.resize(mNbPoint);
3630  if(pos->second.mDataType==WAVELENGTH_TOF)
3631  {
3632  this->SetRadiationType(RAD_NEUTRON);
3633  this->GetRadiation().SetWavelengthType(WAVELENGTH_TOF);
3635  }
3636  else this->GetRadiation().SetWavelength(pos->second.mWavelength);
3637  for(unsigned long i=0;i<mNbPoint;++i)
3638  {
3639  mPowderPatternObs(i)=pos->second.mPowderPatternObs[i];
3640  mX(i)=pos->second.mPowderPatternX[i];
3641  mPowderPatternObsSigma(i)=pos->second.mPowderPatternSigma[i];
3642  }
3643  this->SetWeightToInvSigmaSq();
3644  this->SetPowderPatternX(mX);
3645  }
3646  VFN_DEBUG_EXIT("PowderPattern::ImportPowderPatternCIF():file:",5)
3647 }
3648 
3649 
3650 void PowderPattern::SetPowderPatternObs(const CrystVector_REAL& obs)
3651 {
3652  VFN_DEBUG_MESSAGE("PowderPattern::ImportPowderPatternObs()",5)
3653  if((unsigned long)obs.numElements() != mNbPoint)
3654  {
3655  cout << obs.numElements()<<" "<<mNbPoint<<" "<<this<<endl;
3656  throw(ObjCrystException("PowderPattern::SetPowderPatternObs(vect): The \
3657 supplied vector of observed intensities does not have the expected number of points!"));
3658  }
3659 
3660  mPowderPatternObs=obs;
3661  mPowderPatternObsSigma.resize(mPowderPatternObs.numElements());
3663 
3664  this->SetSigmaToSqrtIobs();
3665  this->SetWeightToInvSigmaSq();
3666  mClockIntegratedFactorsPrep.Reset();
3667  {
3668  char buf [200];
3669  sprintf(buf,"Changed powder pattern: %d points, 2theta=%7.3f -> %7.3f, step=%6.3f",
3670  (int)mNbPoint,mX(0)*RAD2DEG,mX(mNbPoint-1)*RAD2DEG,
3671  (mX(mNbPoint-1)-mX(0))/(mNbPoint-1)*RAD2DEG);
3672  (*fpObjCrystInformUser)((string)buf);
3673  }
3674 }
3675 void PowderPattern::SavePowderPattern(const string &filename) const
3676 {
3677  VFN_DEBUG_MESSAGE("PowderPattern::SavePowderPattern",5)
3678  this->CalcPowderPattern();
3679  ofstream out(filename.c_str());
3680  CrystVector_REAL ttheta;
3681  ttheta=mX;
3682  if(this->GetRadiation().GetWavelengthType()!=WAVELENGTH_TOF) ttheta *= RAD2DEG;
3683 
3684  CrystVector_REAL diff;
3685  diff=mPowderPatternObs;
3686  diff-=mPowderPatternCalc;
3687  out << "# 2Theta/TOF Iobs ICalc Iobs-Icalc Weight Comp0" << endl;
3688  out << FormatVertVector<REAL>(ttheta,
3691  diff,mPowderPatternWeight,
3692  mPowderPatternComponentRegistry.GetObj(0).mPowderPatternCalc,16,8);
3693  out.close();
3694  VFN_DEBUG_MESSAGE("DiffractionDataPowder::SavePowderPattern:End",3)
3695 }
3696 
3697 void PowderPattern::PrintObsCalcData(ostream&os)const
3698 {
3699  VFN_DEBUG_MESSAGE("DiffractionDataPowder::PrintObsCalcData()",5);
3700  CrystVector_REAL ttheta;
3701  ttheta=mX;
3702  if(this->GetRadiation().GetWavelengthType()!=WAVELENGTH_TOF) ttheta *= RAD2DEG;
3703  os << "PowderPattern : " << mName <<endl;
3704  os << " 2Theta/TOF Obs Sigma Calc Weight" <<endl;
3705  os << FormatVertVector<REAL>(ttheta,mPowderPatternObs,mPowderPatternObsSigma,
3707  // mPowderPatternComponentRegistry.GetObj(0).mPowderPatternCalc,12,4);
3708 }
3709 
3711 {
3712  if( (0==this->GetPowderPatternObs().numElements())
3713  ||(0==GetNbPowderPatternComponent()))
3714  {
3715  return 0;
3716  }
3717  this->CalcPowderPattern();
3718  TAU_PROFILE("PowderPattern::GetR()","void ()",TAU_DEFAULT);
3719 
3720  REAL tmp1=0.;
3721  REAL tmp2=0.;
3722 
3723  unsigned long maxPoints=mNbPointUsed;
3724  if( (true==mStatisticsExcludeBackground)
3725  &&(mPowderPatternBackgroundCalc.numElements()>0))
3726  {
3727  const REAL *p1, *p2, *p3;
3728  p1=mPowderPatternCalc.data();
3729  p2=mPowderPatternObs.data();
3730  p3=mPowderPatternBackgroundCalc.data();
3731  const long nbExclude=mExcludedRegionMinX.numElements();
3732  if(0==nbExclude)
3733  {
3734  VFN_DEBUG_MESSAGE("PowderPattern::GetR():Exclude Backgd",4);
3735  for(unsigned long i=0;i<maxPoints;i++)
3736  {
3737  tmp1 += ((*p1)-(*p2)) * ((*p1)-(*p2));
3738  tmp2 += ((*p2)-(*p3)) * ((*p2)-(*p3));
3739  p1++;p2++;p3++;
3740  }
3741  }
3742  else
3743  {
3744  VFN_DEBUG_MESSAGE("PowderPattern::GetR():Exclude Backgd,Exclude regions",4);
3745  unsigned long min,max;
3746  unsigned long i=0;
3747  for(int j=0;j<nbExclude;j++)
3748  {
3749  min=(unsigned long)floor(this->X2Pixel(mExcludedRegionMinX(j)));
3750  max=(unsigned long)ceil (this->X2Pixel(mExcludedRegionMaxX(j)));
3751  if(min>maxPoints) break;
3752  if(max>maxPoints)max=maxPoints;
3753  for(;i<min;i++)
3754  {
3755  tmp1 += ((*p1)-(*p2)) * ((*p1)-(*p2));
3756  tmp2 += ((*p2)-(*p3)) * ((*p2)-(*p3));
3757  p1++;p2++;p3++;
3758  }
3759  p1 += max-i;
3760  p2 += max-i;
3761  p3 += max-i;
3762  i += max-i;
3763  }
3764  for(;i<maxPoints;i++)
3765  {
3766  tmp1 += ((*p1)-(*p2)) * ((*p1)-(*p2));
3767  tmp2 += ((*p2)-(*p3)) * ((*p2)-(*p3));
3768  p1++;p2++;p3++;
3769  }
3770 
3771  }
3772  } // Exclude Background ?
3773  else
3774  {
3775  const REAL *p1, *p2;
3776  p1=mPowderPatternCalc.data();
3777  p2=mPowderPatternObs.data();
3778  const long nbExclude=mExcludedRegionMinX.numElements();
3779  if(0==nbExclude)
3780  {
3781  VFN_DEBUG_MESSAGE("PowderPattern::GetR()",4);
3782  for(unsigned long i=0;i<maxPoints;i++)
3783  {
3784  tmp1 += ((*p1)-(*p2))*((*p1)-(*p2));
3785  tmp2 += (*p2) * (*p2);
3786  //cout <<i<<":"<< tmp1 << " "<<tmp2 << " " << *p1 <<" "<<*p2<<endl;
3787  p1++;p2++;
3788  }
3789  }
3790  else
3791  {
3792  VFN_DEBUG_MESSAGE("PowderPattern::GetR(),Exclude regions",4);
3793  unsigned long min,max;
3794  unsigned long i=0;
3795  for(int j=0;j<nbExclude;j++)
3796  {
3797  min=(unsigned long)floor(this->X2Pixel(mExcludedRegionMinX(j)));
3798  max=(unsigned long)ceil (this->X2Pixel(mExcludedRegionMaxX(j)));
3799  if(min>maxPoints) break;
3800  if(max>maxPoints)max=maxPoints;
3801  for(;i<min;i++)
3802  {
3803  tmp1 += ((*p1)-(*p2))*((*p1)-(*p2));
3804  tmp2 += (*p2) * (*p2);
3805  p1++;p2++;
3806  }
3807  p1 += max-i;
3808  p2 += max-i;
3809  i += max-i;
3810  }
3811  for(;i<maxPoints;i++)
3812  {
3813  tmp1 += ((*p1)-(*p2))*((*p1)-(*p2));
3814  tmp2 += (*p2) * (*p2);
3815  p1++;p2++;
3816  }
3817  }
3818  }
3819 
3820  VFN_DEBUG_MESSAGE("PowderPattern::GetR()="<<sqrt(tmp1/tmp2),4);
3821  //cout << FormatVertVector<REAL>(mPowderPatternCalc,mPowderPatternObs);
3822  //this->SavePowderPattern("refinedPattern.out");
3823  //abort();
3824  return sqrt(tmp1/tmp2);
3825 }
3826 
3827 REAL PowderPattern::GetIntegratedR()const
3828 {
3829  if( (0==this->GetPowderPatternObs().numElements())
3830  ||(0==GetNbPowderPatternComponent()))
3831  {
3832  return 0;
3833  }
3834  this->CalcPowderPattern();
3835  this->PrepareIntegratedRfactor();
3836  VFN_DEBUG_ENTRY("PowderPattern::GetIntegratedR()",4);
3837  TAU_PROFILE("PowderPattern::GetIntegratedR()","void ()",TAU_DEFAULT);
3838 
3839  REAL tmp1=0.;
3840  REAL tmp2=0.;
3841  const long numInterval=mIntegratedPatternMin.numElements();
3842  if( (true==mStatisticsExcludeBackground)
3843  &&(mPowderPatternBackgroundCalc.numElements()>0))
3844  {
3845  const REAL *p1, *p2, *p3;
3846  CrystVector_REAL integratedCalc(numInterval);
3847  integratedCalc=0;
3848  CrystVector_REAL backgdCalc(numInterval);
3849  backgdCalc=0;
3850  REAL *pp1=integratedCalc.data();
3851  REAL *pp2=backgdCalc.data();
3852  for(int i=0;i<numInterval;i++)
3853  {
3854  const long max=mIntegratedPatternMax(i);
3855  p1=mPowderPatternCalc.data()+mIntegratedPatternMin(i);
3856  for(int j=mIntegratedPatternMin(i);j<=max;j++) *pp1 += *p1++;
3857  pp1++;
3858  p1=mPowderPatternBackgroundCalc.data()+mIntegratedPatternMin(i);
3859  for(int j=mIntegratedPatternMin(i);j<=max;j++) *pp2 += *p1++;
3860  pp2++;
3861  }
3862 
3863  p1=integratedCalc.data();
3864  p2=mIntegratedObs.data();
3865  p3=backgdCalc.data();
3866  VFN_DEBUG_MESSAGE("PowderPattern::GetIntegratedR():Exclude Backgd",2);
3867  for(long i=0;i<numInterval;i++)
3868  {
3869  tmp1 += ((*p1)-(*p2)) * ((*p1)-(*p2));
3870  tmp2 += ((*p2)-(*p3)) * ((*p2)-(*p3));
3871  p1++;p2++;p3++;
3872  }
3873  } // Exclude Background ?
3874  else
3875  {
3876  const REAL *p1, *p2;
3877  CrystVector_REAL integratedCalc(numInterval);
3878  integratedCalc=0;
3879  REAL *pp1=integratedCalc.data();
3880  for(int i=0;i<numInterval;i++)
3881  {
3882  const long max=mIntegratedPatternMax(i);
3883  p1=mPowderPatternCalc.data()+mIntegratedPatternMin(i);
3884  for(int j=mIntegratedPatternMin(i);j<=max;j++) *pp1 += *p1++;
3885  pp1++;
3886  }
3887  p1=integratedCalc.data();
3888  p2=mIntegratedObs.data();
3889  VFN_DEBUG_MESSAGE("PowderPattern::GetIntegratedR()",2);
3890  for(long i=0;i<numInterval;i++)
3891  {
3892  tmp1 += ((*p1)-(*p2))*((*p1)-(*p2));
3893  tmp2 += (*p2) * (*p2);
3894  //cout <<i<<":"<< tmp1 << " "<<tmp2 << " " << *p1 <<" "<<*p2<<endl;
3895  p1++;p2++;
3896  }
3897  }
3898 
3899  VFN_DEBUG_EXIT("PowderPattern::GetIntegratedR()="<<sqrt(tmp1/tmp2),4);
3900  return sqrt(tmp1/tmp2);
3901 }
3902 
3904 {
3905  if( (0==this->GetPowderPatternObs().numElements())
3906  ||(0==GetNbPowderPatternComponent()))
3907  {
3908  return 0;
3909  }
3910  this->CalcPowderPattern();
3911  TAU_PROFILE("PowderPattern::GetRw()","void ()",TAU_DEFAULT);
3912  VFN_DEBUG_MESSAGE("PowderPattern::GetRw()",3);
3913 
3914 
3915  //cout <<FormatVertVector<REAL>(mPowderPatternObs,
3916  // mPowderPatternCalc,
3917  // mPowderPatternWeight);
3918  REAL tmp1=0.;
3919  REAL tmp2=0.;
3920 
3921  unsigned long maxPoints=mNbPointUsed;
3922 
3923  if( (true==mStatisticsExcludeBackground)
3924  &&(mPowderPatternBackgroundCalc.numElements()>0))
3925  {
3926  VFN_DEBUG_MESSAGE("PowderPattern::GetRw():Exclude Backgd",3);
3927  const REAL *p1, *p2, *p3, *p4;
3928  p1=mPowderPatternCalc.data();
3929  p2=mPowderPatternObs.data();
3930  p3=mPowderPatternBackgroundCalc.data();
3931  p4=mPowderPatternWeight.data();
3932  const long nbExclude=mExcludedRegionMinX.numElements();
3933  if(0==nbExclude)
3934  {
3935  for(unsigned long i=0;i<maxPoints;i++)
3936  {
3937  tmp1 += *p4 * ((*p1)-(*p2)) * ((*p1)-(*p2));
3938  tmp2 += *p4++ * ((*p2)-(*p3)) * ((*p2)-(*p3));
3939  p1++;p2++;p3++;
3940  }
3941  }
3942  else
3943  {
3944  unsigned long min,max;
3945  unsigned long i=0;
3946  for(int j=0;j<nbExclude;j++)
3947  {
3948  min=(unsigned long)floor(this->X2Pixel(mExcludedRegionMinX(j)));
3949  max=(unsigned long)ceil (this->X2Pixel(mExcludedRegionMaxX(j)));
3950  if(min>maxPoints) break;
3951  if(max>maxPoints)max=maxPoints;
3952  for(;i<min;i++)
3953  {
3954  tmp1 += *p4 * ((*p1)-(*p2)) * ((*p1)-(*p2));
3955  tmp2 += *p4++ * ((*p2)-(*p3)) * ((*p2)-(*p3));
3956  p1++;p2++;p3++;
3957  }
3958  p1 += max-i;
3959  p2 += max-i;
3960  p3 += max-i;
3961  p4 += max-i;
3962  i += max-i;
3963  }
3964  for(;i<maxPoints;i++)
3965  {
3966  tmp1 += *p4 * ((*p1)-(*p2)) * ((*p1)-(*p2));
3967  tmp2 += *p4++ * ((*p2)-(*p3)) * ((*p2)-(*p3));
3968  p1++;p2++;p3++;
3969  }
3970 
3971  }
3972  }
3973  else
3974  {
3975  VFN_DEBUG_MESSAGE("PowderPattern::GetRw()",3);
3976  const REAL *p1, *p2, *p4;
3977  p1=mPowderPatternCalc.data();
3978  p2=mPowderPatternObs.data();
3979  p4=mPowderPatternWeight.data();
3980  const long nbExclude=mExcludedRegionMinX.numElements();
3981  if(0==nbExclude)
3982  {
3983  for(unsigned long i=0;i<maxPoints;i++)
3984  {
3985  tmp1 += *p4 * ((*p1)-(*p2))*((*p1)-(*p2));
3986  tmp2 += *p4++ * (*p2) * (*p2);
3987  p1++;p2++;
3988  }
3989  }
3990  else
3991  {
3992  unsigned long min,max;
3993  unsigned long i=0;
3994  for(int j=0;j<nbExclude;j++)
3995  {
3996  min=(unsigned long)floor(this->X2Pixel(mExcludedRegionMinX(j)));
3997  max=(unsigned long)ceil (this->X2Pixel(mExcludedRegionMaxX(j)));
3998  if(min>maxPoints) break;
3999  if(max>maxPoints)max=maxPoints;
4000  for(;i<min;i++)
4001  {
4002  tmp1 += *p4 * ((*p1)-(*p2))*((*p1)-(*p2));
4003  tmp2 += *p4++ * (*p2) * (*p2);
4004  p1++;p2++;
4005  }
4006  p1 += max-i;
4007  p2 += max-i;
4008  p4 += max-i;
4009  i += max-i;
4010  }
4011  for(;i<maxPoints;i++)
4012  {
4013  tmp1 += *p4 * ((*p1)-(*p2))*((*p1)-(*p2));
4014  tmp2 += *p4++ * (*p2) * (*p2);
4015  p1++;p2++;
4016  }
4017  }
4018  }
4019  VFN_DEBUG_MESSAGE("PowderPattern::GetRw()="<<sqrt(tmp1/tmp2),3);
4020  return sqrt(tmp1/tmp2);
4021 }
4022 REAL PowderPattern::GetIntegratedRw()const
4023 {
4024  if( (0==this->GetPowderPatternObs().numElements())
4025  ||(0==GetNbPowderPatternComponent()))
4026  {
4027  return 0;
4028  }
4029  this->CalcPowderPattern();
4030  this->PrepareIntegratedRfactor();
4031  TAU_PROFILE("PowderPattern::GetIntegratedRw()","void ()",TAU_DEFAULT);
4032 
4033  REAL tmp1=0.;
4034  REAL tmp2=0.;
4035  const long numInterval=mIntegratedPatternMin.numElements();
4036  if( (true==mStatisticsExcludeBackground)
4037  &&(mPowderPatternBackgroundCalc.numElements()>0))
4038  {
4039  const REAL *p1, *p2, *p3, *p4;
4040  CrystVector_REAL integratedCalc(numInterval);
4041  integratedCalc=0;
4042  CrystVector_REAL backgdCalc(numInterval);
4043  backgdCalc=0;
4044  REAL *pp1=integratedCalc.data();
4045  REAL *pp2=backgdCalc.data();
4046  for(int i=0;i<numInterval;i++)
4047  {
4048  const long max=mIntegratedPatternMax(i);
4049  p1=mPowderPatternCalc.data()+mIntegratedPatternMin(i);
4050  for(int j=mIntegratedPatternMin(i);j<=max;j++) *pp1 += *p1++;
4051  pp1++;
4052  p1=mPowderPatternBackgroundCalc.data()+mIntegratedPatternMin(i);
4053  for(int j=mIntegratedPatternMin(i);j<=max;j++) *pp2 += *p1++;
4054  pp2++;
4055  }
4056 
4057  p1=integratedCalc.data();
4058  p2=mIntegratedObs.data();
4059  p3=backgdCalc.data();
4060  if(mIntegratedWeight.numElements()==0) p4=mIntegratedWeightObs.data();
4061  else p4=mIntegratedWeight.data();
4062  VFN_DEBUG_MESSAGE("PowderPattern::GetIntegratedRw():Exclude Backgd",4);
4063  for(long i=0;i<numInterval;i++)
4064  {
4065  tmp1 += *p4 * ((*p1)-(*p2)) * ((*p1)-(*p2));
4066  tmp2 += *p4++ * ((*p2)-(*p3)) * ((*p2)-(*p3));
4067  //cout <<i<<": " <<mIntegratedPatternMin(i)<<"->"<<mIntegratedPatternMax(i)
4068  // <<" "<< tmp1 << " "<<tmp2 << " " << *p1 <<" "<<*p2<<" "<<*p3<<" "<<*(p4-1) <<endl;
4069  p1++;p2++;p3++;
4070  }
4071  } // Exclude Background ?
4072  else
4073  {
4074  const REAL *p1, *p2, *p4;
4075  CrystVector_REAL integratedCalc(numInterval);
4076  integratedCalc=0;
4077  REAL *pp1=integratedCalc.data();
4078  for(int i=0;i<numInterval;i++)
4079  {
4080  const long max=mIntegratedPatternMax(i);
4081  p1=mPowderPatternCalc.data()+mIntegratedPatternMin(i);
4082  for(int j=mIntegratedPatternMin(i);j<=max;j++) *pp1 += *p1++;
4083  pp1++;
4084  }
4085  p1=integratedCalc.data();
4086  p2=mIntegratedObs.data();
4087  if(mIntegratedWeight.numElements()==0) p4=mIntegratedWeightObs.data();
4088  else p4=mIntegratedWeight.data();
4089  VFN_DEBUG_MESSAGE("PowderPattern::GetIntegratedRw()",4);
4090  for(long i=0;i<numInterval;i++)
4091  {
4092  tmp1 += *p4 * ((*p1)-(*p2))*((*p1)-(*p2));
4093  tmp2 += *p4++ * (*p2) * (*p2);
4094  //cout <<i<<":"<< tmp1 << " "<<tmp2 << " " << *p1 <<" "<<*p2<<endl;
4095  p1++;p2++;
4096  }
4097  }
4098 
4099  VFN_DEBUG_MESSAGE("PowderPattern::GetIntegratedRw()="<<sqrt(tmp1/tmp2),4);
4100  //cout << FormatVertVector<REAL>(mPowderPatternCalc,mPowderPatternObs);
4101  //this->SavePowderPattern("refinedPattern.out");
4102  //abort();
4103  return sqrt(tmp1/tmp2);
4104 }
4105 
4107 {
4108  if( (0==this->GetPowderPatternObs().numElements())
4109  ||(0==GetNbPowderPatternComponent()))
4110  {
4111  mChi2=0.;
4112  return mChi2;
4113  }
4114  this->CalcNbPointUsed();
4115  if(mClockChi2>mClockMaster) return mChi2;
4116 
4117  this->CalcPowderPattern();
4118  if( (mClockChi2>mClockPowderPatternPar)
4119  &&(mClockChi2>mClockScaleFactor)
4120  &&(mClockChi2>mClockPowderPatternCalc)) return mChi2;
4121  // We want the best scale factor
4122  this->FitScaleFactorForRw();
4123 
4124  TAU_PROFILE("PowderPattern::GetChi2()","void ()",TAU_DEFAULT);
4125 
4126  VFN_DEBUG_ENTRY("PowderPattern::GetChi2()",3);
4127 
4128  const unsigned long maxPoints=mNbPointUsed;
4129 
4130  mChi2=0.;
4131  mChi2LikeNorm=0.;
4132  VFN_DEBUG_MESSAGE("PowderPattern::GetChi2()Integrated profiles",3);
4133  const REAL * RESTRICT p1, * RESTRICT p2, * RESTRICT p3;
4134  p1=mPowderPatternCalc.data();
4135  p2=mPowderPatternObs.data();
4136  p3=mPowderPatternWeight.data();
4137  const long nbExclude=mExcludedRegionMinX.numElements();
4138  if(0==nbExclude)
4139  {
4140  for(unsigned long i=0;i<maxPoints;i++)
4141  {
4142  mChi2 += *p3 * ((*p1)-(*p2))*((*p1)-(*p2));
4143  if(*p3<=0) p3++;
4144  else mChi2LikeNorm -= log(*p3++);
4145  p1++;p2++;
4146  }
4147  }
4148  else
4149  {
4150  unsigned long min,max;
4151  unsigned long i=0;
4152  for(int j=0;j<nbExclude;j++)
4153  {
4154  min=(unsigned long)floor(this->X2Pixel(mExcludedRegionMinX(j)));
4155  max=(unsigned long)ceil (this->X2Pixel(mExcludedRegionMaxX(j)));
4156  if(min>maxPoints) break;
4157  if(max>maxPoints)max=maxPoints;
4158  for(;i<min;i++)
4159  {
4160  mChi2 += *p3 * ((*p1)-(*p2))*((*p1)-(*p2));
4161  if(*p3<=0) p3++;
4162  else mChi2LikeNorm -= log(*p3++);
4163  p1++;p2++;
4164  }
4165  p1 += max-i;
4166  p2 += max-i;
4167  p3 += max-i;
4168  i += max-i;
4169  }
4170  for(;i<maxPoints;i++)
4171  {
4172  mChi2 += *p3 * ((*p1)-(*p2))*((*p1)-(*p2));
4173  if(*p3<=0) p3++;
4174  else mChi2LikeNorm -= log(*p3++);
4175  p1++;p2++;
4176  }
4177  }
4178  mChi2LikeNorm/=2;
4179  VFN_DEBUG_MESSAGE("Chi^2="<<mChi2<<", log(norm)="<<mChi2LikeNorm,3)
4180  mClockChi2.Click();
4181  VFN_DEBUG_EXIT("PowderPattern::GetChi2()="<<mChi2,3);
4182  return mChi2;
4183 }
4184 
4186 {
4187  if( (0==this->GetPowderPatternObs().numElements())
4188  ||(0==GetNbPowderPatternComponent()))
4189  {
4190  mIntegratedChi2=0.;
4191  return mIntegratedChi2;
4192  }
4193  this->CalcNbPointUsed();
4194  if(mClockIntegratedChi2>mClockMaster) return mIntegratedChi2;
4195 
4197  if( (mClockChi2>mClockPowderPatternPar)
4198  &&(mClockChi2>mClockScaleFactor)
4199  &&(mClockChi2>mClockPowderPatternIntegratedCalc)) return mIntegratedChi2;
4200 
4201  // We want the best scale factor
4202  this->FitScaleFactorForIntegratedRw();
4203 
4204  TAU_PROFILE("PowderPattern::GetIntegratedChi2()","void ()",TAU_DEFAULT);
4205 
4206  VFN_DEBUG_ENTRY("PowderPattern::GetIntegratedChi2()",3);
4207 
4208  mIntegratedChi2=0.;
4209  mIntegratedChi2LikeNorm=0.;
4210  VFN_DEBUG_MESSAGE("PowderPattern::GetIntegratedChi2()",3);
4211  const REAL * RESTRICT p1, * RESTRICT p2, * RESTRICT p3;
4212  p1=mPowderPatternIntegratedCalc.data();
4213  p2=mIntegratedObs.data();
4214  if(mIntegratedWeight.numElements()==0) p3=mIntegratedWeightObs.data();
4215  else p3=mIntegratedWeight.data();
4216  double weightProd=1.;
4217  VFN_DEBUG_MESSAGE("PowderPattern::GetIntegratedIntegratedRw()",4);
4218  for(unsigned long i=0;i<mNbIntegrationUsed;)
4219  {
4220  // group weights to avoid computing too many log()
4221  // group only a limited number to avoid underflow...
4222  for(unsigned long j=0;j<32;++j)
4223  {
4224  mIntegratedChi2 += *p3 * ((*p1)-(*p2))*((*p1)-(*p2));
4225  if(*p3>0) weightProd *= *p3;
4226  p1++;p2++;p3++;
4227  if(++i == mNbIntegrationUsed) break;
4228  }
4229  mIntegratedChi2LikeNorm -= log(weightProd);
4230  weightProd=1.;
4231  }
4232  mIntegratedChi2LikeNorm/=2;
4233  VFN_DEBUG_MESSAGE("Chi^2="<<mIntegratedChi2<<", log(norm)="<<mIntegratedChi2LikeNorm,3)
4234  mClockIntegratedChi2.Click();
4235  VFN_DEBUG_EXIT("PowderPattern::GetChi2()="<<mIntegratedChi2,3);
4236  return mIntegratedChi2;
4237 }
4238 
4240 {
4241  if(0 == mOptProfileIntegration.GetChoice()) return this->GetIntegratedChi2();
4242  else return this->GetChi2();
4243 }
4244 
4246 {
4247  if( (0==this->GetPowderPatternObs().numElements())
4248  ||(0==GetNbPowderPatternComponent()))
4249  {
4250  return ;
4251  }
4252  this->CalcPowderPattern();
4253  TAU_PROFILE("PowderPattern::FitScaleFactorForR()","void ()",TAU_DEFAULT);
4254  VFN_DEBUG_ENTRY("PowderPattern::FitScaleFactorForR()",3);
4255  // Which components are scalable ?
4256  mScalableComponentIndex.resize(mPowderPatternComponentRegistry.GetNb());
4257  int nbScale=0;
4258  for(int i=0;i<mPowderPatternComponentRegistry.GetNb();i++)
4259  {
4260  if(mPowderPatternComponentRegistry.GetObj(i).IsScalable())
4261  mScalableComponentIndex(nbScale++)=i;
4262  }
4263  VFN_DEBUG_MESSAGE("-> Number of Scale Factors:"<<nbScale<<":Index:"<<endl<<mScalableComponentIndex,3);
4264  if(0==nbScale)
4265  {
4266  VFN_DEBUG_EXIT("PowderPattern::FitScaleFactorForR(): No scalable component!",3);
4267  return;
4268  }
4269  mScalableComponentIndex.resizeAndPreserve(nbScale);
4270  // prepare matrices
4271  mFitScaleFactorM.resize(nbScale,nbScale);
4272  mFitScaleFactorB.resize(nbScale,1);
4273  mFitScaleFactorX.resize(nbScale,1);
4274  // Build Matrix & Vector for LSQ
4275  const long nbExclude=mExcludedRegionMinX.numElements();
4276  if(0==nbExclude)
4277  {
4278  for(int i=0;i<nbScale;i++)
4279  {
4280  for(int j=i;j<nbScale;j++)
4281  {
4282  // Here use a direct access to the powder spectrum, since
4283  // we know it has just been recomputed
4284  const REAL *p1=mPowderPatternComponentRegistry.GetObj(mScalableComponentIndex(i))
4285  .mPowderPatternCalc.data();
4286  const REAL *p2=mPowderPatternComponentRegistry.GetObj(mScalableComponentIndex(j))
4287  .mPowderPatternCalc.data();
4288  REAL m=0.;
4289  for(unsigned long k=0;k<mNbPointUsed;k++) m += *p1++ * *p2++;
4290  mFitScaleFactorM(i,j)=m;
4291  mFitScaleFactorM(j,i)=m;
4292  }
4293  }
4294  for(int i=0;i<nbScale;i++)
4295  {
4296  const REAL *p1=mPowderPatternObs.data();
4297  const REAL *p2=mPowderPatternComponentRegistry.GetObj(mScalableComponentIndex(i))
4298  .mPowderPatternCalc.data();
4299  REAL b=0.;
4300  if(mPowderPatternBackgroundCalc.numElements()<=1)
4301  for(unsigned long k=0;k<mNbPointUsed;k++) b += *p1++ * *p2++;
4302  else
4303  {
4304  const REAL *p3=mPowderPatternBackgroundCalc.data();
4305  for(unsigned long k=0;k<mNbPointUsed;k++) b += (*p1++ - *p3++) * *p2++;
4306  }
4307  mFitScaleFactorB(i,0) =b;
4308  }
4309  }
4310  else
4311  {
4312  unsigned long min,max;
4313  for(int i=0;i<nbScale;i++)
4314  {
4315  for(int j=i;j<nbScale;j++)
4316  {
4317  // Here use a direct access to the powder spectrum, since
4318  // we know it has just been recomputed
4319  const REAL *p1=mPowderPatternComponentRegistry.GetObj(mScalableComponentIndex(i))
4320  .mPowderPatternCalc.data();
4321  const REAL *p2=mPowderPatternComponentRegistry.GetObj(mScalableComponentIndex(j))
4322  .mPowderPatternCalc.data();
4323  REAL m=0.;
4324  unsigned long l=0;
4325  for(int k=0;k<nbExclude;k++)
4326  {
4327  min=(unsigned long)floor(this->X2Pixel(mExcludedRegionMinX(j)));
4328  max=(unsigned long)ceil (this->X2Pixel(mExcludedRegionMaxX(j)));
4329  if(min>mNbPointUsed) break;
4330  if(max>mNbPointUsed)max=mNbPointUsed;
4332  for(;l<min;l++) m += *p1++ * *p2++;
4333  p1 += max-l;
4334  p2 += max-l;
4335  l = max;
4336  }
4337  for(;l<mNbPointUsed;l++) m += *p1++ * *p2++;
4338  mFitScaleFactorM(i,j)=m;
4339  mFitScaleFactorM(j,i)=m;
4340  }
4341  }
4342  for(int i=0;i<nbScale;i++)
4343  {
4344  const REAL *p1=mPowderPatternObs.data();
4345  const REAL *p2=mPowderPatternComponentRegistry
4346  .GetObj(mScalableComponentIndex(i))
4347  .mPowderPatternCalc.data();
4348  REAL b=0.;
4349  unsigned long l=0;
4350  if(mPowderPatternBackgroundCalc.numElements()<=1)
4351  {
4352  for(int k=0;k<nbExclude;k++)
4353  {
4354  min=(unsigned long)floor(this->X2Pixel(mExcludedRegionMinX(k)));
4355  max=(unsigned long)ceil (this->X2Pixel(mExcludedRegionMaxX(k)));
4356  if(min>mNbPointUsed) break;
4357  if(max>mNbPointUsed)max=mNbPointUsed;
4359  for(;l<min;l++) b += *p1++ * *p2++;
4360  p1 += max-l;
4361  p2 += max-l;
4362  l = max;
4363  }
4364  for(;l<mNbPointUsed;l++) b += *p1++ * *p2++;
4365  }
4366  else
4367  {
4368  const REAL *p3=mPowderPatternBackgroundCalc.data();
4369  for(int k=0;k<nbExclude;k++)
4370  {
4371  min=(unsigned long)floor(this->X2Pixel(mExcludedRegionMinX(k)));
4372  max=(unsigned long)ceil (this->X2Pixel(mExcludedRegionMaxX(k)));
4373  if(min>mNbPointUsed) break;
4374  if(max>mNbPointUsed)max=mNbPointUsed;
4376  for(;l<min;l++) b += (*p1++ - *p3++) * *p2++;
4377  p1 += max-l;
4378  p2 += max-l;
4379  l = max;
4380  }
4381  for(;l<mNbPointUsed;l++) b += (*p1++ - *p3++) * *p2++;
4382  }
4383  mFitScaleFactorB(i,0) =b;
4384  }
4385  }
4386  if(1==nbScale) mFitScaleFactorX=mFitScaleFactorB(0)/mFitScaleFactorM(0);
4387  else
4388  mFitScaleFactorX=product(InvertMatrix(mFitScaleFactorM),mFitScaleFactorB);
4389  VFN_DEBUG_MESSAGE("B, M, X"<<endl<<mFitScaleFactorB<<endl<<mFitScaleFactorM<<endl<<mFitScaleFactorX,2)
4390  for(int i=0;i<nbScale;i++)
4391  {
4392  const REAL * p1=mPowderPatternComponentRegistry.GetObj(mScalableComponentIndex(i))
4393  .mPowderPatternCalc.data();
4394  REAL * p0 = mPowderPatternCalc.data();
4395  const REAL s = mFitScaleFactorX(i)
4396  -mScaleFactor(mScalableComponentIndex(i));
4397  if(ISNAN_OR_INF(s))
4398  {
4399  (*fpObjCrystInformUser)("Warning: working around NaN scale factor...");
4400  continue;
4401  }
4402  for(unsigned long j=0;j<mNbPointUsed;j++) *p0++ += s * *p1++;
4403  VFN_DEBUG_MESSAGE("-> Old:"<<mScaleFactor(mScalableComponentIndex(i)) <<" Change:"<<mFitScaleFactorX(i),2);
4404  mScaleFactor(mScalableComponentIndex(i)) = mFitScaleFactorX(i);
4406  mClockPowderPatternCalc.Click();//we *did* correct the spectrum
4407  }
4408  VFN_DEBUG_EXIT("PowderPattern::FitScaleFactorForR():End",3);
4409 }
4410 
4411 void PowderPattern::FitScaleFactorForIntegratedR()const
4412 {
4413  if( (0==this->GetPowderPatternObs().numElements())
4414  ||(0==GetNbPowderPatternComponent()))
4415  {
4416  return ;
4417  }
4418  this->CalcPowderPattern();
4419  this->PrepareIntegratedRfactor();
4420  VFN_DEBUG_ENTRY("PowderPattern::FitScaleFactorForIntegratedR()",3);
4421  TAU_PROFILE("PowderPattern::FitScaleFactorForIntegratedR()","void ()",TAU_DEFAULT);
4422  // Which components are scalable ?
4423  mScalableComponentIndex.resize(mPowderPatternComponentRegistry.GetNb());
4424  int nbScale=0;
4425  for(int i=0;i<mPowderPatternComponentRegistry.GetNb();i++)
4426  {
4427  if(mPowderPatternComponentRegistry.GetObj(i).IsScalable())
4428  mScalableComponentIndex(nbScale++)=i;
4429  }
4430  VFN_DEBUG_MESSAGE("-> Number of Scale Factors:"<<nbScale<<":Index:"<<endl<<mScalableComponentIndex,2);
4431  if(0==nbScale)
4432  {
4433  VFN_DEBUG_EXIT("PowderPattern::FitScaleFactorForIntegratedR(): No scalable component!",3);
4434  return;
4435  }
4436  mScalableComponentIndex.resizeAndPreserve(nbScale);
4437  // prepare matrices
4438  mFitScaleFactorM.resize(nbScale,nbScale);
4439  mFitScaleFactorB.resize(nbScale,1);
4440  mFitScaleFactorX.resize(nbScale,1);
4441  // Build Matrix & Vector for LSQ
4442  VFN_DEBUG_MESSAGE("PowderPattern::FitScaleFactorForIntegratedR():1",2);
4443  const long numInterval=mIntegratedPatternMin.numElements();
4444  CrystVector_REAL *integratedCalc= new CrystVector_REAL[nbScale];
4445  for(int i=0;i<nbScale;i++)
4446  {
4447  integratedCalc[i].resize(numInterval);
4448 
4449  // Here use a direct access to the powder spectrum, since
4450  // we know it has just been recomputed
4451  const REAL *p1=mPowderPatternComponentRegistry.GetObj(mScalableComponentIndex(i))
4452  .mPowderPatternCalc.data();
4453 
4454  REAL *p2=integratedCalc[i].data();
4455  for(int j=0;j<numInterval;j++)
4456  {
4457  const long max=mIntegratedPatternMax(j);
4458  p1=mPowderPatternComponentRegistry.GetObj(mScalableComponentIndex(i))
4459  .mPowderPatternCalc.data()+mIntegratedPatternMin(j);
4460  *p2=0;
4461  for(int k=mIntegratedPatternMin(j);k<=max;k++) *p2 += *p1++;
4462  //cout <<"Calc#"<<i<<":"<< mIntegratedPatternMin(j) << " "
4463  // <<mIntegratedPatternMax(j)<<" "
4464  // << *p2<<endl;
4465  p2++;
4466  }
4467  }
4468  VFN_DEBUG_MESSAGE("PowderPattern::FitScaleFactorForIntegratedR():2",2);
4469  CrystVector_REAL backdIntegrated(numInterval);
4470  if(mPowderPatternBackgroundCalc.numElements()>1)
4471  {
4472  const REAL *p1;
4473  REAL *p2=backdIntegrated.data();
4474  for(int j=0;j<numInterval;j++)
4475  {
4476  const long max=mIntegratedPatternMax(j);
4477  p1=mPowderPatternBackgroundCalc.data()+mIntegratedPatternMin(j);
4478  *p2=0;
4479  for(int k=mIntegratedPatternMin(j);k<=max;k++) *p2 += *p1++;
4480  //cout <<"Backgd:"<< mIntegratedPatternMin(j) << " "
4481  // <<mIntegratedPatternMax(j)<<" "
4482  // << *p2<<endl;
4483  p2++;
4484  }
4485  }
4486 
4487  //if(mPowderPatternBackgroundCalc.numElements()<=1)
4488  // cout<< FormatVertVector<REAL>(integratedCalc[0],mIntegratedObs,mIntegratedWeight,backdIntegrated)<<endl;
4489  //else
4490  // cout<< FormatVertVector<REAL>(integratedCalc[0],mIntegratedObs,mIntegratedWeight)<<endl;
4491  VFN_DEBUG_MESSAGE("PowderPattern::FitScaleFactorForIntegratedR():3",2);
4492  for(int i=0;i<nbScale;i++)
4493  {
4494  for(int j=i;j<nbScale;j++)
4495  {
4496  const REAL *p1=integratedCalc[i].data();
4497  const REAL *p2=integratedCalc[j].data();
4498  REAL m=0.;
4499  for(long k=0;k<numInterval;k++)
4500  {
4501  m += *p1++ * *p2++;
4502  //cout <<"M:"<< mIntegratedPatternMin(k) << " "<<mIntegratedPatternMax(k)<<" "<<m<<endl;
4503  }
4504  mFitScaleFactorM(i,j)=m;
4505  mFitScaleFactorM(j,i)=m;
4506  }
4507  }
4508  VFN_DEBUG_MESSAGE("PowderPattern::FitScaleFactorForIntegratedR():4",2);
4509  for(int i=0;i<nbScale;i++)
4510  {
4511  const REAL *p1=mIntegratedObs.data();
4512  const REAL *p2=integratedCalc[i].data();
4513  REAL b=0.;
4514  if(mPowderPatternBackgroundCalc.numElements()<=1)
4515  for(long k=0;k<numInterval;k++)
4516  {
4517  b += *p1++ * *p2++;
4518  //cout<<"B:"<<mIntegratedPatternMin(k)<<" "<<mIntegratedPatternMax(k)<<" "<<b<<endl;
4519  }
4520  else
4521  {
4522  const REAL *p3=backdIntegrated.data();
4523  for(long k=0;k<numInterval;k++)
4524  {
4525  //cout<<"B(minus backgd):"<<mIntegratedPatternMin(k)<<" "
4526  // <<mIntegratedPatternMax(k)<<" "
4527  // <<*p1<<" "<<*p2<<" "<<*p3<<" "<<b<<endl;
4528  b += (*p1++ - *p3++) * *p2++;
4529  }
4530  }
4531  mFitScaleFactorB(i,0) =b;
4532  }
4533  VFN_DEBUG_MESSAGE("PowderPattern::FitScaleFactorForIntegratedR():5",2);
4534 
4535  if(1==nbScale) mFitScaleFactorX=mFitScaleFactorB(0)/mFitScaleFactorM(0);
4536  else
4537  mFitScaleFactorX=product(InvertMatrix(mFitScaleFactorM),mFitScaleFactorB);
4538  VFN_DEBUG_MESSAGE("B, M, X"<<endl<<mFitScaleFactorB<<endl<<mFitScaleFactorM<<endl<<mFitScaleFactorX,3)
4539  for(int i=0;i<nbScale;i++)
4540  {
4541  const REAL * p1=mPowderPatternComponentRegistry.GetObj(mScalableComponentIndex(i))
4542  .mPowderPatternCalc.data();
4543  REAL * p0 = mPowderPatternCalc.data();
4544  const REAL s = mFitScaleFactorX(i)
4545  -mScaleFactor(mScalableComponentIndex(i));
4546  if(ISNAN_OR_INF(s))
4547  {
4548  (*fpObjCrystInformUser)("Warning: working around NaN scale factor...");
4549  continue;
4550  }
4551  for(unsigned long j=0;j<mNbPointUsed;j++) *p0++ += s * *p1++;
4552  VFN_DEBUG_MESSAGE("-> Old:"<<mScaleFactor(mScalableComponentIndex(i)) <<" New:"<<mFitScaleFactorX(i),3);
4553  mScaleFactor(mScalableComponentIndex(i)) = mFitScaleFactorX(i);
4555  mClockPowderPatternCalc.Click();//we *did* correct the spectrum
4556  }
4557  delete[] integratedCalc;
4558  VFN_DEBUG_EXIT("PowderPattern::FitScaleFactorForIntegratedR():End",3);
4559 }
4560 
4562 {
4563  if( (0==this->GetPowderPatternObs().numElements())
4564  ||(0==GetNbPowderPatternComponent()))
4565  {
4566  return ;
4567  }
4568  TAU_PROFILE("PowderPattern::FitScaleFactorForRw()","void ()",TAU_DEFAULT);
4569  VFN_DEBUG_ENTRY("PowderPattern::FitScaleFactorForRw()",3);
4570  this->CalcPowderPattern();
4571  //:TODO: take into account excluded regions...
4572  // Which components are scalable ?
4573  mScalableComponentIndex.resize(mPowderPatternComponentRegistry.GetNb());
4574  int nbScale=0;
4575  for(int i=0;i<mPowderPatternComponentRegistry.GetNb();i++)
4576  {
4577  if(mPowderPatternComponentRegistry.GetObj(i).IsScalable())
4578  mScalableComponentIndex(nbScale++)=i;
4579  }
4580  VFN_DEBUG_MESSAGE("-> Number of Scale Factors:"<<nbScale<<":Index:"<<endl<<mScalableComponentIndex,2);
4581  if(0==nbScale)
4582  {
4583  VFN_DEBUG_EXIT("PowderPattern::FitScaleFactorForRw(): No scalable component!",3);
4584  return;
4585  }
4586  mScalableComponentIndex.resizeAndPreserve(nbScale);
4587  // prepare matrices
4588  mFitScaleFactorM.resize(nbScale,nbScale);
4589  mFitScaleFactorB.resize(nbScale,1);
4590  mFitScaleFactorX.resize(nbScale,1);
4591  // Build Matrix & Vector for LSQ
4592  const long nbExclude=mExcludedRegionMinX.numElements();
4593  if(0==nbExclude)
4594  {
4595  for(int i=0;i<nbScale;i++)
4596  {
4597  for(int j=i;j<nbScale;j++)
4598  {
4599  // Here use a direct access to the powder spectrum, since
4600  // we know it has just been recomputed
4601  const REAL *p1=mPowderPatternComponentRegistry.GetObj(mScalableComponentIndex(i))
4602  .mPowderPatternCalc.data();
4603  const REAL *p2=mPowderPatternComponentRegistry.GetObj(mScalableComponentIndex(j))
4604  .mPowderPatternCalc.data();
4605  const REAL *p3=mPowderPatternWeight.data();
4606  REAL m=0.;
4607  for(unsigned long k=0;k<mNbPointUsed;k++) m += *p1++ * *p2++ * *p3++;
4608  mFitScaleFactorM(i,j)=m;
4609  mFitScaleFactorM(j,i)=m;
4610  }
4611  }
4612  for(int i=0;i<nbScale;i++)
4613  {
4614  const REAL *p1=mPowderPatternObs.data();
4615  const REAL *p2=mPowderPatternComponentRegistry.GetObj(mScalableComponentIndex(i))
4616  .mPowderPatternCalc.data();
4617  const REAL *p3=mPowderPatternWeight.data();
4618  REAL b=0.;
4619  if(mPowderPatternBackgroundCalc.numElements()<=1)
4620  for(unsigned long k=0;k<mNbPointUsed;k++) b += *p1++ * *p2++ * *p3++;
4621  else
4622  {
4623  const REAL *p4=mPowderPatternBackgroundCalc.data();
4624  for(unsigned long k=0;k<mNbPointUsed;k++)
4625  b += (*p1++ - *p4++) * *p2++ * *p3++;
4626  }
4627  mFitScaleFactorB(i,0) =b;
4628  }
4629  }
4630  else
4631  {
4632  unsigned long min,max;
4633  for(int i=0;i<nbScale;i++)
4634  {
4635  for(int j=i;j<nbScale;j++)
4636  {
4637  // Here use a direct access to the powder spectrum, since
4638  // we know it has just been recomputed
4639  const REAL *p1=mPowderPatternComponentRegistry.GetObj(mScalableComponentIndex(i))
4640  .mPowderPatternCalc.data();
4641  const REAL *p2=mPowderPatternComponentRegistry.GetObj(mScalableComponentIndex(j))
4642  .mPowderPatternCalc.data();
4643  const REAL *p3=mPowderPatternWeight.data();
4644  REAL m=0.;
4645  unsigned long l=0;
4646  for(int k=0;k<nbExclude;k++)
4647  {
4648  min=(unsigned long)floor(this->X2Pixel(mExcludedRegionMinX(j)));
4649  max=(unsigned long)ceil (this->X2Pixel(mExcludedRegionMaxX(j)));
4650  if(min>mNbPointUsed) break;
4651  if(max>mNbPointUsed)max=mNbPointUsed;
4653  for(;l<min;l++) m += *p1++ * *p2++ * *p3++;
4654  p1 += max-l;
4655  p2 += max-l;
4656  p3 += max-l;
4657  l = max;
4658  }
4659  for(;l<mNbPointUsed;l++) m += *p1++ * *p2++ * *p3++;
4660  mFitScaleFactorM(i,j)=m;
4661  mFitScaleFactorM(j,i)=m;
4662  }
4663  }
4664  for(int i=0;i<nbScale;i++)
4665  {
4666  const REAL *p1=mPowderPatternObs.data();
4667  const REAL *p2=mPowderPatternComponentRegistry
4668  .GetObj(mScalableComponentIndex(i))
4669  .mPowderPatternCalc.data();
4670  const REAL *p3=mPowderPatternWeight.data();
4671  REAL b=0.;
4672  unsigned long l=0;
4673  if(mPowderPatternBackgroundCalc.numElements()<=1)
4674  {
4675  for(int k=0;k<nbExclude;k++)
4676  {
4677  min=(unsigned long)floor(this->X2Pixel(mExcludedRegionMinX(k)));
4678  max=(unsigned long)ceil (this->X2Pixel(mExcludedRegionMaxX(k)));
4679  if(min>mNbPointUsed) break;
4680  if(max>mNbPointUsed)max=mNbPointUsed;
4682  for(;l<min;l++) b += *p1++ * *p2++ * *p3++;
4683  p1 += max-l;
4684  p2 += max-l;
4685  p3 += max-l;
4686  l = max;
4687  }
4688  for(;l<mNbPointUsed;l++) b += *p1++ * *p2++ * *p3++;
4689  }
4690  else
4691  {
4692  const REAL *p4=mPowderPatternBackgroundCalc.data();
4693  for(int k=0;k<nbExclude;k++)
4694  {
4695  min=(unsigned long)floor(this->X2Pixel(mExcludedRegionMinX(k)));
4696  max=(unsigned long)ceil (this->X2Pixel(mExcludedRegionMaxX(k)));
4697  if(min>mNbPointUsed) break;
4698  if(max>mNbPointUsed)max=mNbPointUsed;
4700  for(;l<min;l++) b += (*p1++ - *p4++) * *p2++ * *p3++;
4701  p1 += max-l;
4702  p2 += max-l;
4703  p3 += max-l;
4704  l = max;
4705  }
4706  for(;l<mNbPointUsed;l++) b += (*p1++ - *p4++) * *p2++ * *p3++;
4707  }
4708  mFitScaleFactorB(i,0) =b;
4709  }
4710  }
4711  if(1==nbScale) mFitScaleFactorX=mFitScaleFactorB(0)/mFitScaleFactorM(0);
4712  else
4713  mFitScaleFactorX=product(InvertMatrix(mFitScaleFactorM),mFitScaleFactorB);
4714  VFN_DEBUG_MESSAGE("B, M, X"<<endl<<mFitScaleFactorB<<endl<<mFitScaleFactorM<<endl<<mFitScaleFactorX,2)
4715  for(int i=0;i<nbScale;i++)
4716  {
4717  const REAL * p1=mPowderPatternComponentRegistry.GetObj(mScalableComponentIndex(i))
4718  .mPowderPatternCalc.data();
4719  REAL * p0 = mPowderPatternCalc.data();
4720  const REAL s = mFitScaleFactorX(i)
4721  -mScaleFactor(mScalableComponentIndex(i));
4722  if(ISNAN_OR_INF(s))
4723  {
4724  (*fpObjCrystInformUser)("Warning: working around NaN scale factor...");
4725  continue;
4726  }
4727  for(unsigned long j=0;j<mNbPointUsed;j++) *p0++ += s * *p1++;
4728  VFN_DEBUG_MESSAGE("-> Old:"<<mScaleFactor(mScalableComponentIndex(i)) <<" Change:"<<mFitScaleFactorX(i),3);
4729  mScaleFactor(mScalableComponentIndex(i)) = mFitScaleFactorX(i);
4731  mClockPowderPatternCalc.Click();//we *did* correct the spectrum
4732  }
4733  VFN_DEBUG_EXIT("PowderPattern::FitScaleFactorForRw():End",3);
4734 }
4735 
4736 void PowderPattern::FitScaleFactorForIntegratedRw()const
4737 {
4738  if( (0==this->GetPowderPatternObs().numElements())
4739  ||(0==GetNbPowderPatternComponent()))
4740  {
4741  return ;
4742  }
4745  VFN_DEBUG_ENTRY("PowderPattern::FitScaleFactorForIntegratedRw()",3);
4746  TAU_PROFILE("PowderPattern::FitScaleFactorForIntegratedRw()","void ()",TAU_DEFAULT);
4747  // Which components are scalable ? Which scalable calculated components have a non-null variance ?
4748  mScalableComponentIndex.resize(mPowderPatternComponentRegistry.GetNb());
4749  CrystVector_REAL vScalableVarianceIndex(mPowderPatternComponentRegistry.GetNb());
4750  int nbScale=0;
4751  int nbVarCalc=0;
4752  for(int i=0;i<mPowderPatternComponentRegistry.GetNb();i++)
4753  {
4754  if(mPowderPatternComponentRegistry.GetObj(i).IsScalable())
4755  mScalableComponentIndex(nbScale++)=i;
4756  if(mPowderPatternComponentRegistry.GetObj(i).HasPowderPatternCalcVariance())
4757  {
4758  if(0!=mPowderPatternComponentRegistry.GetObj(i)
4759  .GetPowderPatternIntegratedCalcVariance().first->numElements())
4760  vScalableVarianceIndex(nbVarCalc++)=i;
4761  }
4762  }
4763  VFN_DEBUG_MESSAGE("-> Number of Scale Factors:"<<nbScale<<"("<<nbVarCalc
4764  <<"with variance). Index:"<<endl<<mScalableComponentIndex,2);
4765  if(0==nbScale)
4766  {
4767  VFN_DEBUG_EXIT("PowderPattern::FitScaleFactorForIntegratedRw(): No scalable component!",3);
4768  return;
4769  }
4770  bool again;
4771  unsigned int ctagain=0;
4772  RECALC_SCALE_FACTOR_VARIANCE_FitScaleFactorForIntegratedRw:
4773  if(false)//if((nbScale==1)&&(nbVarCalc==1))
4774  {// Special case when using ML error, the scale appears also in the weight so we have to
4775  // use a 2nd-degree equation (ax^2+bx+c=0) to get the scale factor.
4776  double a=0.,b=0.,c=0.;// possible overflows...
4777  REAL newscale;
4778  {
4779  const REAL * RESTRICT p1=mIntegratedObs.data();
4780  const REAL * RESTRICT p2=mPowderPatternComponentRegistry.GetObj(mScalableComponentIndex(0))
4781  .GetPowderPatternIntegratedCalc().first->data();
4782  const REAL * RESTRICT p1v=mIntegratedVarianceObs.data();
4783  const REAL * RESTRICT p2v=mPowderPatternComponentRegistry.GetObj(mScalableComponentIndex(0))
4784  .GetPowderPatternIntegratedCalcVariance().first->data();
4785  if(mPowderPatternBackgroundIntegratedCalc.numElements()<=1)
4786  {
4787  for(unsigned long k=mNbIntegrationUsed;k>0;k--)
4788  {
4789  a += *p2v * *p1 * *p2;
4790  b += *p2 * *p2 * *p1v - *p1 * *p1 * *p2v;
4791  c -= *p1 * *p2 * *p1v;
4792  p1++;p2++;p1v++;p2v++;
4793  }
4794  }
4795  else
4796  {
4797  const REAL * RESTRICT p3=mPowderPatternBackgroundIntegratedCalc.data();
4798  for(unsigned long k=mNbIntegrationUsed;k>0;k--)
4799  {
4800  a += *p2v * (*p1 - *p3) * *p2;
4801  b += *p2 * *p2 * *p1v - (*p1 - *p3) * (*p1 - *p3) * *p2v;
4802  c -= (*p1 - *p3) * *p2 * *p1v;
4803  p1++;p2++;p1v++;p2v++;p3++;
4804  }
4805  }
4806  // Only one >0 solution to the 2nd degree equation
4807  newscale=(REAL)((-b+sqrt(b*b-4*a*c))/(2*a));
4808  }
4809  {// Store new scale, and correct old calculated pattern
4810  const REAL s = newscale-mScaleFactor(mScalableComponentIndex(0));
4811  const REAL s2 = newscale*newscale
4812  -mScaleFactor(mScalableComponentIndex(0))
4813  *mScaleFactor(mScalableComponentIndex(0));
4814  REAL * RESTRICT p0 = mPowderPatternIntegratedCalc.data();
4815  const REAL * RESTRICT p1=mPowderPatternComponentRegistry.GetObj(mScalableComponentIndex(0))
4816  .GetPowderPatternIntegratedCalc().first->data();
4817  REAL * RESTRICT p0v = mPowderPatternVarianceIntegrated.data();
4818  const REAL * RESTRICT p1v=mPowderPatternComponentRegistry.GetObj(mScalableComponentIndex(0))
4819  .GetPowderPatternIntegratedCalcVariance().first->data();
4820  REAL * RESTRICT p0w = mIntegratedWeight.data();
4821  for(unsigned long j=mNbIntegrationUsed;j>0;j--)
4822  {
4823  *p0++ += s * *p1++;
4824  *p0v += s2 * *p1v++;
4825  if(*p0v <=0) {*p0w++ =0;p0v++;}else *p0w++ = 1. / *p0v++;
4826  }
4827  }
4828  mScaleFactor(mScalableComponentIndex(0)) = newscale;
4830  mClockPowderPatternIntegratedCalc.Click();//we *did* correct the spectrum
4831  }
4832  else
4833  {
4834  again=false;
4835  mScalableComponentIndex.resizeAndPreserve(nbScale);
4836  // prepare matrices
4837  mFitScaleFactorM.resize(nbScale,nbScale);
4838  mFitScaleFactorB.resize(nbScale,1);
4839  mFitScaleFactorX.resize(nbScale,1);
4840  // Build Matrix & Vector for LSQ
4841  VFN_DEBUG_MESSAGE("PowderPattern::FitScaleFactorForIntegratedRw():1",2);
4842  vector<const CrystVector_REAL*> integratedCalc;
4843  for(int i=0;i<nbScale;i++)
4844  {
4845  integratedCalc.push_back(mPowderPatternComponentRegistry.
4846  GetObj(mScalableComponentIndex(i))
4847  .GetPowderPatternIntegratedCalc().first);
4848  }
4849  VFN_DEBUG_MESSAGE("PowderPattern::FitScaleFactorForIntegratedRw():3",2);
4850  for(int i=0;i<nbScale;i++)
4851  {
4852  for(int j=i;j<nbScale;j++)
4853  {
4854  const REAL * RESTRICT p1=integratedCalc[i]->data();
4855  const REAL * RESTRICT p2=integratedCalc[j]->data();
4856  const REAL * RESTRICT p3;
4857  if(mIntegratedWeight.numElements()==0)
4858  {
4859  p3=mIntegratedWeightObs.data();
4860  if(ctagain>5) VFN_DEBUG_MESSAGE("ctagain="<<ctagain<<", using mIntegratedWeightObs",5);
4861  }
4862  else
4863  {
4864  p3=mIntegratedWeight.data();
4865  if(ctagain>5) VFN_DEBUG_MESSAGE("ctagain="<<ctagain<<", using mIntegratedWeight",5);
4866  }
4867  REAL m=0.;
4868  if(j==i)
4869  for(unsigned long k=mNbIntegrationUsed;k>0;k--)
4870  {m += *p1 * *p1 * *p3++; p1++;}
4871  else
4872  for(unsigned long k=mNbIntegrationUsed;k>0;k--)
4873  m += *p1++ * *p2++ * *p3++;
4874  mFitScaleFactorM(i,j)=m;
4875  mFitScaleFactorM(j,i)=m;
4876  }
4877  }
4878  VFN_DEBUG_MESSAGE("PowderPattern::FitScaleFactorForIntegratedRw():4",2);
4879  for(int i=0;i<nbScale;i++)
4880  {
4881  const REAL * RESTRICT p1=mIntegratedObs.data();
4882  const REAL * RESTRICT p2=integratedCalc[i]->data();
4883  const REAL * RESTRICT p4;
4884  if(mIntegratedWeight.numElements()==0) p4=mIntegratedWeightObs.data();
4885  else p4=mIntegratedWeight.data();
4886  REAL b=0.;
4887  if(mPowderPatternBackgroundIntegratedCalc.numElements()<=1)
4888  {
4889  for(unsigned long k=mNbIntegrationUsed;k>0;k--)
4890  {
4891  b += *p1++ * *p2++ * *p4++;
4892  //cout<<"B:"<<mIntegratedPatternMin(k)<<" "<<mIntegratedPatternMax(k)<<" "<<b<<endl;
4893  }
4894  }
4895  else
4896  {
4897  const REAL * RESTRICT p3=mPowderPatternBackgroundIntegratedCalc.data();
4898  for(unsigned long k=mNbIntegrationUsed;k>0;k--)
4899  {
4900  //cout<<"B(minus backgd):"<<mIntegratedPatternMin(k)<<" "
4901  // <<mIntegratedPatternMax(k)<<" "
4902  // <<*p1<<" "<<*p2<<" "<<*p3<<" "<<b<<endl;
4903  b += (*p1++ - *p3++) * *p2++ * *p4++;
4904  }
4905  }
4906  mFitScaleFactorB(i,0) =b;
4907  }
4908  VFN_DEBUG_MESSAGE("PowderPattern::FitScaleFactorForIntegratedRw():5",2);
4909 
4910  if(1==nbScale) mFitScaleFactorX=mFitScaleFactorB(0)/mFitScaleFactorM(0);
4911  else if(1<nbScale)
4912  mFitScaleFactorX=product(InvertMatrix(mFitScaleFactorM),mFitScaleFactorB);
4913  VFN_DEBUG_MESSAGE("B, M, X"<<endl<<mFitScaleFactorB<<endl<<mFitScaleFactorM<<endl<<mFitScaleFactorX,3)
4914  //Correct the variance, if necessary
4915  if(nbVarCalc>0)
4916  {
4917  //if(ctagain>0) cout <<"ctgain, sumvar="<<log(mPowderPatternVarianceIntegrated.sum());
4918  for(int i=0;i<nbScale;i++)
4919  {
4920  if(mPowderPatternComponentRegistry.GetObj(mScalableComponentIndex(i)).HasPowderPatternCalcVariance())
4921  {
4922  if(0!=mPowderPatternComponentRegistry.GetObj(mScalableComponentIndex(i))
4923  .GetPowderPatternIntegratedCalcVariance().first->numElements())
4924  {
4925  const REAL * RESTRICT p1=mPowderPatternComponentRegistry.GetObj(mScalableComponentIndex(i))
4926  .GetPowderPatternIntegratedCalcVariance().first->data();
4927  //cout <<",sumvar(i)="<<log(mPowderPatternComponentRegistry.GetObj(mScalableComponentIndex(i))
4928  // .GetPowderPatternIntegratedCalcVariance().first->sum());
4929  REAL * RESTRICT p0 = mPowderPatternVarianceIntegrated.data();
4930  const REAL s2 = mFitScaleFactorX(i)*mFitScaleFactorX(i)
4931  -mScaleFactor(mScalableComponentIndex(i))
4932  *mScaleFactor(mScalableComponentIndex(i));
4933  for(unsigned long j=mNbIntegrationUsed;j>0;j--) *p0++ += s2 * *p1++;
4934  }
4935  }
4936  }
4937  //if(ctagain>0) cout <<" ->"<<log(mPowderPatternVarianceIntegrated.sum())
4938  // <<" , sumobsvar="<<log(mIntegratedVarianceObs.sum())<<endl;
4939  REAL * RESTRICT p0 = mIntegratedWeight.data();
4940  const REAL * RESTRICT p1=mPowderPatternVarianceIntegrated.data();
4941  for(unsigned long j=mNbIntegrationUsed;j>0;j--)
4942  if(*p1 <=0) {*p0++ =0;p1++;}
4943  else *p0++ = 1. / *p1++;
4944 
4945  }
4946  // Correct the calculated integrated pattern
4947  for(int i=0;i<nbScale;i++)
4948  {
4949  const REAL * RESTRICT p1=mPowderPatternComponentRegistry.GetObj(mScalableComponentIndex(i))
4950  .GetPowderPatternIntegratedCalc().first->data();
4951  REAL * RESTRICT p0 = mPowderPatternIntegratedCalc.data();
4952  const REAL s = mFitScaleFactorX(i)
4953  -mScaleFactor(mScalableComponentIndex(i));
4954  if(ISNAN_OR_INF(s))
4955  {
4956  (*fpObjCrystInformUser)("Warning: working around NaN scale factor...");
4957  continue;
4958  }
4959  if(nbVarCalc>0)
4960  {
4961  if(abs(s/mFitScaleFactorX(i))>0.001)
4962  {
4963  again=true;
4964  //cout<<"log(scale) :"<<log(mScaleFactor(mScalableComponentIndex(i)))
4965  // <<"->"<<log(mFitScaleFactorX(i))<<" Again="<<ctagain++<<endl;
4966  }
4967  if((!again)&&(ctagain>0))
4968  {
4969  VFN_DEBUG_MESSAGE("log(scale) :"<<log(mScaleFactor(mScalableComponentIndex(i)))
4970  <<"->"<<log(mFitScaleFactorX(i))<<" Again="<<ctagain,5);
4971  }
4972  }
4973  for(unsigned long j=mNbIntegrationUsed;j>0;j--) *p0++ += s * *p1++;
4974  if(ctagain>5)
4975  {
4976  VFN_DEBUG_MESSAGE("->log(scale) Old :"<<log(mScaleFactor(mScalableComponentIndex(i))) <<" New:"<<log(mFitScaleFactorX(i)),10);
4977  }
4978  mScaleFactor(mScalableComponentIndex(i)) = mFitScaleFactorX(i);
4979  }
4980 
4982  mClockPowderPatternIntegratedCalc.Click();//we *did* correct the spectrum
4983  }
4984  if(again && (ctagain<20))
4985  {
4986  ctagain++;
4987  if(ctagain>5)
4988  {
4989  VFN_DEBUG_MESSAGE("PowderPattern::FitScaleFactorForIntegratedRw(), scaling again #"<<ctagain,10)
4990  }
4991  goto RECALC_SCALE_FACTOR_VARIANCE_FitScaleFactorForIntegratedRw;
4992  }
4993  VFN_DEBUG_EXIT("PowderPattern::FitScaleFactorForIntegratedRw():End",3);
4994 }
4995 
4997 {
4998  VFN_DEBUG_MESSAGE("PowderPattern::SetSigmaToSqrtIobs()",5);
4999  for(long i=0;i<mPowderPatternObs.numElements();i++)
5000  mPowderPatternObsSigma(i)=sqrt(fabs(mPowderPatternObs(i)));
5001 }
5002 void PowderPattern::SetWeightToInvSigmaSq(const REAL minRelatSigma)
5003 {
5004  VFN_DEBUG_MESSAGE("PowderPattern::SetWeightToInvSigmaSq()",5);
5005  //:KLUDGE: If less than 1e-4*max, set to 0.... Do not give weight to unobserved points
5006  const REAL min=MaxAbs(mPowderPatternObsSigma)*minRelatSigma;
5007  //mPowderPatternWeight.resize(mPowderPatternObsSigma.numElements());
5008  REAL tmp;
5009  for(long i=0;i<mPowderPatternObsSigma.numElements();i++)
5010  {
5011  tmp = mPowderPatternObsSigma(i);
5012  if(tmp<min) mPowderPatternWeight(i)= 1./min/min;
5013  else mPowderPatternWeight(i) =1./tmp/tmp;
5014  }
5015 }
5017 {
5018  VFN_DEBUG_MESSAGE("PowderPattern::SetWeightToSinTheta()",5);
5019  //mPowderPatternWeight.resize(mPowderPatternObs.numElements());
5021 }
5022 void PowderPattern::SetWeightPolynomial(const REAL a, const REAL b,
5023  const REAL c,
5024  const REAL minRelatIobs)
5025 {
5026  VFN_DEBUG_MESSAGE("PowderPattern::SetWeightPolynomial()",5);
5027  const REAL min=MaxAbs(mPowderPatternObs)*minRelatIobs;
5028  REAL tmp;
5029  for(long i=0;i<mPowderPatternWeight.numElements();i++)
5030  {
5031  tmp=mPowderPatternObs(i);
5032  if(tmp<min)
5033  mPowderPatternWeight(i) = 1./(a+min+b*min*min+c*min*min*min);
5034  else mPowderPatternWeight(i) = 1./(a+tmp+b*tmp*tmp+c*tmp*tmp*tmp);
5035  }
5036 }
5037 
5038 void PowderPattern::BeginOptimization(const bool allowApproximations,
5039  const bool enableRestraints)
5040 {
5041  this->Prepare();
5042  if(0 == mOptProfileIntegration.GetChoice()) this->FitScaleFactorForIntegratedRw();
5043  else this->FitScaleFactorForRw();
5044  this->RefinableObj::BeginOptimization(allowApproximations,enableRestraints);
5045 }
5046 
5047 //void PowderPattern::SetApproximationFlag(const bool allow)
5048 //{// Do we need this ?
5049 // this->Prepare();
5050 // if(0 == mOptProfileIntegration.GetChoice()) this->FitScaleFactorForIntegratedRw();
5051 // else this->FitScaleFactorForRw();
5052 // this->RefinableObj::SetApproximationFlag(allow);
5053 //}
5054 
5055 void PowderPattern::GlobalOptRandomMove(const REAL mutationAmplitude,
5056  const RefParType *type)
5057 {
5058  if(mRandomMoveIsDone) return;
5059  this->RefinableObj::GlobalOptRandomMove(mutationAmplitude,type);
5060 }
5061 
5062 void PowderPattern::AddExcludedRegion(const REAL min,const REAL max)
5063 {
5064  VFN_DEBUG_MESSAGE("PowderPattern::AddExcludedRegion()",5)
5065  const int num=mExcludedRegionMinX.numElements();
5066  if(num>0)
5067  {
5068  mExcludedRegionMinX.resizeAndPreserve(num+1);
5069  mExcludedRegionMaxX.resizeAndPreserve(num+1);
5070  }
5071  else
5072  {
5073  mExcludedRegionMinX.resize(1);
5074  mExcludedRegionMaxX.resize(1);
5075  }
5076  mExcludedRegionMinX(num)=min;
5077  mExcludedRegionMaxX(num)=max;
5078 
5079  //ensure regions are sorted by ascending 2thetamin
5080  CrystVector_long subs;
5081  subs=SortSubs(mExcludedRegionMinX);
5082  CrystVector_REAL tmp1,tmp2;
5083  tmp1=mExcludedRegionMinX;
5084  tmp2=mExcludedRegionMaxX;
5085  for(int i=0;i<mExcludedRegionMinX.numElements();i++)
5086  {
5087  mExcludedRegionMinX(i)=tmp1(subs(i));
5088  mExcludedRegionMaxX(i)=tmp2(subs(i));
5089  }
5091  VFN_DEBUG_MESSAGE("PowderPattern::Add2ThetaExcludedRegion():End",5)
5092 }
5093 
5095 {
5096  REAL tmp=this->GetChi2_Option();
5097  if(mOptProfileIntegration.GetChoice()==0) tmp+=mIntegratedChi2LikeNorm;
5098  else tmp+=mChi2LikeNorm;
5099  return tmp;
5100 }
5101 
5102 unsigned int PowderPattern::GetNbLSQFunction()const{return 2;}
5103 
5104 const CrystVector_REAL&
5105  PowderPattern::GetLSQCalc(const unsigned int idx) const
5106 {
5107  TAU_PROFILE("PowderPattern::GetLSQCalc()","void ()",TAU_DEFAULT);
5108  switch(idx)
5109  {
5110  case 1:
5111  {
5114  break;
5115  }
5116  default:
5117  {
5119  mPowderPatternUsedCalc.resizeAndPreserve(mNbPointUsed);
5120  break;
5121  }
5122  }
5123  return mPowderPatternUsedCalc;
5124 }
5125 
5126 const CrystVector_REAL&
5127  PowderPattern::GetLSQObs(const unsigned int idx) const
5128 {
5129  TAU_PROFILE("PowderPattern::GetLSQObs()","void ()",TAU_DEFAULT);
5130  switch(idx)
5131  {
5132  case 1:
5133  {
5134  this->PrepareIntegratedRfactor();
5135  mPowderPatternUsedObs=mIntegratedObs;
5136  break;
5137  }
5138  default:
5139  {
5141  mPowderPatternUsedObs.resizeAndPreserve(mNbPointUsed);
5142  break;
5143  }
5144  }
5145  return mPowderPatternUsedObs;
5146 }
5147 
5148 const CrystVector_REAL&
5149  PowderPattern::GetLSQWeight(const unsigned int idx) const
5150 {
5151  TAU_PROFILE("PowderPattern::GetLSQWeight()","void ()",TAU_DEFAULT);
5152  switch(idx)
5153  {
5154  case 1:
5155  {
5156  this->PrepareIntegratedRfactor();
5157  // :KLUDGE: When variance is used, mIntegratedWeight will change at each powder pattern calculation,
5158  // so this might be quite wrong...
5159  if(mIntegratedWeight.numElements()==0)
5160  mPowderPatternUsedWeight=mIntegratedWeightObs;
5161  else mPowderPatternUsedWeight=mIntegratedWeight;
5162  break;
5163  }
5164  default:
5165  {
5167  mPowderPatternUsedWeight.resizeAndPreserve(mNbPointUsed);
5168  break;
5169  }
5170  }
5171  return mPowderPatternUsedWeight;
5172 }
5173 
5174 std::map<RefinablePar*, CrystVector_REAL>& PowderPattern::GetLSQ_FullDeriv(const unsigned int idx,std::set<RefinablePar *> &vPar)
5175 {
5176  TAU_PROFILE("PowderPattern::GetLSQ_FullDeriv()","void ()",TAU_DEFAULT);
5177  //return this->RefinableObj::GetLSQ_FullDeriv(idx,vPar);
5178  if(idx==1)
5179  {
5180  this->CalcPowderPatternIntegrated_FullDeriv(vPar);
5181  #if 0
5182  std::map<RefinablePar*, CrystVector_REAL> fullderiv_old;
5183  std::vector<const CrystVector_REAL*> v;
5184  int n=0;
5185  //cout<<"PowderPattern::GetLSQ_FullDeriv(integrated):scales:"<<mScaleFactor<<endl;
5186  cout<<"PowderPattern::GetLSQ_FullDeriv(integrated):parameters:"<<endl;
5187  for(std::set<RefinablePar*>::iterator par=vPar.begin();par!=vPar.end();++par)
5188  {
5189  v.push_back(&(mPowderPatternIntegrated_FullDeriv[*par]));
5190  fullderiv_old[*par]=this->GetLSQDeriv(idx,*(*par));
5191  v.push_back(&(fullderiv_old[*par]));
5192  cout<<(*par)->GetName()<<":"<<mPowderPatternIntegrated_FullDeriv[*par].size()<<","<<fullderiv_old[*par].size()<<endl;
5193  if(++n>8) break;
5194  }
5195  cout<<"PowderPattern::GetLSQ_FullDeriv(integrated):"<<endl<<FormatVertVector<REAL>(v,12,1,20)<<endl;
5196  //exit(0);
5197  #endif
5198  return mPowderPatternIntegrated_FullDeriv;
5199  }
5200  mPowderPattern_FullDeriv=this->GetPowderPattern_FullDeriv(vPar);
5201  return mPowderPattern_FullDeriv;
5202 }
5203 
5205 {
5206  VFN_DEBUG_MESSAGE("PowderPattern::Prepare()",5);
5207  for(int i=0;i<mPowderPatternComponentRegistry.GetNb();i++)
5208  {
5209  mPowderPatternComponentRegistry.GetObj(i).SetMaxSinThetaOvLambda(mMaxSinThetaOvLambda);
5210  mPowderPatternComponentRegistry.GetObj(i).Prepare();
5211  }
5212 }
5214  CrystVector_uint & groupIndex,
5215  unsigned int &first) const
5216 {
5217  // One group for scales, one for theta error parameters
5218  unsigned int scaleIndex=0;
5219  unsigned int thetaIndex=0;
5220  VFN_DEBUG_MESSAGE("PowderPattern::GetGeneGroup()",4)
5221  for(long i=0;i<obj.GetNbPar();i++)
5222  for(long j=0;j<this->GetNbPar();j++)
5223  if(&(obj.GetPar(i)) == &(this->GetPar(j)))
5224  {
5226  {
5227  if(scaleIndex==0) scaleIndex=first++;
5228  groupIndex(i)=scaleIndex;
5229  }
5230  else //gpRefParTypeScattDataCorrPos
5231  {
5232  if(thetaIndex==0) thetaIndex=first++;
5233  groupIndex(i)=thetaIndex;
5234  }
5235  }
5236 }
5237 
5239 
5241 
5242 const CrystVector_long& PowderPattern::GetIntegratedProfileMin()const
5243 {
5244  this->PrepareIntegratedRfactor();
5245  return mIntegratedPatternMin;
5246 }
5247 
5248 const CrystVector_long& PowderPattern::GetIntegratedProfileMax()const
5249 {
5250  this->PrepareIntegratedRfactor();
5251  return mIntegratedPatternMax;
5252 }
5253 
5255 {
5256  return mClockIntegratedFactorsPrep;
5257 }
5258 
5259 REAL PowderPattern::STOL2X(const REAL stol)const
5260 {
5261  REAL x;
5262  if(this->GetRadiation().GetWavelengthType()==WAVELENGTH_TOF)
5263  {
5264  if(stol>0) x = 1.0/(2*stol); else return 0;
5265  x = mDIFC*x+mDIFA*x*x;
5266  VFN_DEBUG_MESSAGE("PowderPattern::STOL2X("<<stol<<","<<1.0/(2*stol)<<")="<<x,2)
5267  }
5268  else
5269  {
5270  x=stol*this->GetWavelength();
5271  if(abs(x)<1.0) x=2*asin(x); else x=2*M_PI;
5272  }
5273  return x;
5274 }
5275 
5276 REAL PowderPattern::X2STOL(const REAL x)const
5277 {
5278  REAL stol;
5279  if(this->GetRadiation().GetWavelengthType()==WAVELENGTH_TOF)
5280  {
5281  if(abs(mDIFA)>abs(mDIFC*1e-6))
5282  {
5283  const REAL delta=mDIFC*mDIFC+4.0*mDIFA*x;
5284  stol = (-mDIFC+sqrt(delta))/(2.0*mDIFA);
5285  stol = 1/(2.0*stol);
5286  }
5287  else stol=mDIFC/(2.0*x);
5288  }
5289  else
5290  stol=sin(x/2.0)/this->GetWavelength();
5291  VFN_DEBUG_MESSAGE("PowderPattern::X2STOL("<<x<<")="<<stol,2)
5292  return stol;
5293 }
5294 
5295 REAL PowderPattern::STOL2Pixel(const REAL stol)const
5296 {
5297  return this->X2Pixel(this->STOL2X(stol));
5298 }
5299 
5300 PeakList PowderPattern::FindPeaks(const float dmin,const float maxratio,const unsigned int maxpeak)
5301 {
5302  const long nb=this->GetNbPoint() ;
5303  // Limit peak detection to 1.5A resolution
5304  long start,finish;
5305  if(this->GetRadiation().GetWavelengthType()!=WAVELENGTH_TOF)
5306  {
5307  start=1;// do not start at 0, if this is a simulation that really start at theta=0...
5308  for(finish=0;finish<nb;++finish)
5309  {
5310  const REAL d=1/(this->X2STOL(this->GetPowderPatternX()(finish))*2);
5311  if(d<dmin) break;
5312  }
5313  }
5314  else
5315  {
5316  finish=nb-1;
5317  for(start=nb-1;start>=0;--start)
5318  {
5319  const REAL d=1/(this->X2STOL(this->GetPowderPatternX()(start))*2);
5320  if(d<dmin) break;
5321  }
5322  }
5323  // First evaluate approximate width (in number of pixels) of reflections
5324  unsigned int width_golay=10;
5325  {
5326  CrystVector_REAL obs;
5327  const int nbwidth=9;
5328  CrystVector_long width(nbwidth);
5329  width=0;
5330  obs=this->GetPowderPatternObs();
5331  const long nb=obs.numElements();
5332  for(int j=0;j<nbwidth;j++)
5333  {
5334  const long imax=obs.imax(nb/10,nb-1);
5335  const REAL iobs_max=obs(imax);
5336  REAL thres=iobs_max;
5337  long i;
5338  for(i=imax-100;i<(imax+100);++i)
5339  {
5340  if(i<0){i=0;continue;}
5341  if(i>=nb) break;
5342  if(obs(i)<thres) thres=obs(i);
5343  }
5344  thres=(iobs_max+thres)/2;
5345  i=imax;
5346  while(obs(i)>=thres)
5347  {
5348  cout<<obs(i)<<" ";
5349  obs(i--)=0;
5350  width(j)+=1;
5351  if(i<0) break;
5352  }
5353  i=imax+1;
5354  while(obs(i)>=thres)
5355  {
5356  cout<<obs(i)<<" ";
5357  obs(i++)=0;
5358  width(j)+=1;
5359  if(i>=nb) break;
5360  }
5361  cout<<endl<<" => "<<width(j)<<endl;
5362  for(i=imax-width(j)*2;i<(imax+width(j)*2);++i)
5363  {
5364  if(i<0) continue;
5365  if(i>=nb) break;
5366  obs(i)=0;
5367  }
5368  }
5369  cout<<"Width of "<<nbwidth<<" strongest peaks:"<<endl<<width;
5370  width_golay=width(SortSubs(width)(nbwidth/2));
5371  cout<<"median width:"<<width_golay<<endl;
5372  if(width_golay<=4)width_golay=4;
5373  if(width_golay>=16)width_golay=16;
5374  }
5375 
5376  // get 2nd derivative
5377  CrystVector_REAL obsd2;
5378  obsd2=SavitzkyGolay(this->GetPowderPatternObs(),width_golay,2);
5379  const float norm=-obsd2.min();
5380  // Normalize, so that the derivative has the same extent as the observed pattern
5381  obsd2 *= mPowderPatternObs.max()/(-norm);
5382 
5383  REAL min_iobs;
5384  if(maxratio<0)
5385  {//Automatic discrimination - get an idea of noise from the distribution of the scond derivative
5386  CrystVector_REAL tmp;
5387  tmp=obsd2;
5388  tmp.resizeAndPreserve(tmp.numElements()/4);// First quarter, avoid too many peaks
5389  CrystVector<long> sub(tmp.numElements());
5390  QuickSortSubs(tmp,sub,tmp.numElements()-1,0);
5391  min_iobs=5*(tmp(tmp.numElements()/2)-tmp(tmp.numElements()/4));
5392  //cout<<__FILE__<<":"<<__LINE__<<" MIN_IOBS (automatic)="<<min_iobs<<endl;
5393  }else min_iobs=-1;// This will be set after highest peak is found
5394  PeakList pl;
5395  int nbav_min=0;//minimum numerb of points over which the peak is integrated
5396  while(true)
5397  {// Start from max
5398  const long imax=obsd2.imax(start,finish);
5399  REAL iobs=obsd2(imax);
5400  if(iobs<=0) break;
5401  REAL xmax=mX(imax)*iobs;
5402  long nbav=1;
5403  long i=imax;
5404  REAL lastiobs=obsd2(i);
5405  const REAL iobs_max=lastiobs;
5406  //cout<<i<<":"<<lastiobs<<":"<<mX(i)*RAD2DEG<<endl;
5407  while(true)
5408  {
5409  if(i<=1) break;
5410  if(obsd2(--i)>=lastiobs) break;
5411  lastiobs=obsd2(i);
5412  obsd2(i)=0;
5413  iobs+=lastiobs;
5414  xmax+=mX(i)*lastiobs;nbav++;
5415  if(lastiobs<=0) break;
5416  //cout<<i<<":"<<lastiobs<<":"<<mX(i)*RAD2DEG<<endl;
5417  }
5418  float dleft=mX(i+1);
5419  i=imax;
5420  lastiobs=obsd2(i);
5421  while(true)
5422  {
5423  if(i>=(nb-2)) break;
5424  if(obsd2(++i)>=lastiobs) break;
5425  lastiobs=obsd2(i);
5426  obsd2(i)=0;
5427  iobs+=lastiobs;
5428  xmax+=mX(i)*lastiobs;nbav++;
5429  if(lastiobs<=0) break;
5430  //cout<<i<<":"<<lastiobs<<":"<<mX(i)*RAD2DEG<<endl;
5431  }
5432  float dright=mX(i-1);
5433  xmax/=iobs;
5434  obsd2(imax)=0;
5435  REAL dmax=this->X2STOL(xmax)*2;
5436  dright=this->X2STOL(dright)*2;
5437  dleft =this->X2STOL(dleft)*2;
5438  //TODO : evaluate min intensity ratio from noise ?
5439  if(min_iobs<0)
5440  {
5441  cout<<this->GetPowderPatternObsSigma()(imax)<<","<<this->GetPowderPatternObs()(imax)<<endl;
5442  min_iobs=iobs/nbav*maxratio;
5443  //cout<<__FILE__<<":"<<__LINE__<<" MIN_IOBS="<<min_iobs<<endl;
5444  }
5445  if(nbav_min==0)
5446  {
5447  nbav_min=nbav/2;
5448  if(nbav_min<3)nbav_min=3;
5449  }
5450  if(pl.GetPeakList().size()>maxpeak) break;
5451  float sigma=this->GetPowderPatternObsSigma()(imax);
5452  if(this->GetRadiation().GetWavelengthType()!=WAVELENGTH_TOF)
5453  xmax *= RAD2DEG;
5454  cout<<"Peak #"<<pl.GetPeakList().size()<<"imax="<<imax<<", x="<<xmax<<",d="<<1/dmax<<", d2iobs_max="<<iobs_max
5455  <<", d2Iobs="<<iobs<<", nbav="<<nbav<<", min_iobs="<<min_iobs<<",sigma="<<sigma<<endl;
5456  if( ((nbav>=nbav_min) &&(iobs_max>min_iobs)&&((iobs/nbav)>min_iobs))
5457  ||((nbav>=nbav_min) &&(iobs_max>min_iobs)&&((iobs/nbav)>min_iobs*.2)&&((iobs/nbav)>3*sigma))
5458  ||((nbav>=nbav_min/2)&&(iobs_max>min_iobs)&&((iobs/nbav)>min_iobs*2 )&&((iobs/nbav)>6*sigma)))
5459  {
5460  pl.AddPeak(dmax,iobs,abs(dright-dleft)*.25);
5461  if((pl.GetPeakList().size()==1)&&(maxratio<0)&&(min_iobs<0.005*iobs/nbav)) min_iobs=0.005*iobs/nbav;
5462  }
5463  }
5464  pl.Print(cout);
5465  return pl;
5466 }
5467 
5468 const CrystVector_REAL& PowderPattern::GetScaleFactor() const{return mScaleFactor;}
5469 
5470 CrystVector_REAL& PowderPattern::GetScaleFactor(){return mScaleFactor;}
5471 
5472 //Local structures to export atoms, bond and angle restraints
5473 struct exportAtom
5474 {
5475  exportAtom(string n,REAL X, REAL Y, REAL Z,REAL b,REAL o,const ScatteringPower *pow):
5476  name(n),x(X),y(Y),z(Z),biso(b),occ(o),occMult(1),mpScattPow(pow){}
5477  string name;
5478  REAL x,y,z,biso,occ;
5480  int occMult;
5481  const ScatteringPower *mpScattPow;
5482 };
5483 
5484 struct exportBond
5485 {
5486  exportBond(const string &a1,const string &a2, REAL d, REAL s):
5487  at1(a1),at2(a2),dist(d),sigma(s){}
5488  string at1,at2;
5489  REAL dist,sigma;
5490 };
5491 
5492 struct exportAngle
5493 {
5494  exportAngle(const string &a1,const string &a2,const string &a3, REAL a, REAL s):
5495  at1(a1),at2(a2),at3(a3),ang(a),sigma(s){}
5496  string at1,at2,at3;
5497  REAL ang,sigma;
5498 };
5499 
5500 void PowderPattern::ExportFullprof(const std::string &prefix)const
5501 {
5502  // Analyze our data - background ? number of crystalline phases ?
5503  const PowderPatternBackground *pBackground=0;
5504  vector<const PowderPatternDiffraction*> vDiff;
5505  for(unsigned int i=0;i<this->GetNbPowderPatternComponent();i++)
5506  {
5507  if(this->GetPowderPatternComponent(i).GetClassName()=="PowderPatternBackground")
5508  pBackground=dynamic_cast<const PowderPatternBackground*>(&(this->GetPowderPatternComponent(i)));
5509  if(this->GetPowderPatternComponent(i).GetClassName()=="PowderPatternDiffraction")
5510  vDiff.push_back(dynamic_cast<const PowderPatternDiffraction*>(&(this->GetPowderPatternComponent(i))));
5511  }
5512  if((pBackground==0)||vDiff.size()==0) return;
5513 
5514  // Powder data file
5515  ofstream dat((prefix+".dat").c_str());
5516  dat<<"XYDATA"<<endl
5517  <<"INTER 1.0 1.0 0"<<endl<<endl<<endl<<endl<<endl;
5518 
5519  CrystVector_REAL ttheta;
5520  ttheta=mX;
5521  if(this->GetRadiation().GetWavelengthType()!=WAVELENGTH_TOF) ttheta *= RAD2DEG;
5522  dat << FormatVertVector<REAL>(ttheta,mPowderPatternObs,mPowderPatternObsSigma,12,4);
5523  dat.close();
5524 
5525  // PCR file
5526  ofstream pcr((prefix+".pcr").c_str());
5527  // if(!pcr) :TODO:
5528 
5529  // Title
5530  pcr<<"Fox/ObjCryst exported file:"<<this->GetName()<<endl;
5531  // Number of patterns
5532  pcr<<"NPATT 1"<<endl;
5533  // Weight of each pattern
5534  pcr<<"W_PAT 1.0"<<endl;
5535  // Multi-pattern format
5536  pcr<<"! Nph Dum Ias Nre Cry Opt Aut"<<endl;
5537  pcr<<" "<<vDiff.size()
5538  <<" 0 0 0 0 1 1 "<<endl;
5539  // For each phase
5540  {
5541  int job=0;
5542  if(this->GetRadiation().GetRadiationType()==RAD_XRAY) job=0;
5543  if(this->GetRadiation().GetRadiationType()==RAD_NEUTRON) job=1;
5544  //:TODO: TOF:
5545  pcr<<"! Job Npr Nba Nex Nsc Nor Iwg Ilo Res Ste Uni Cor Anm"<<endl
5546  <<" "<<job
5547  <<" 5 "<<pBackground->GetInterpPoints().first->numElements()
5548  <<" 0 0 1 0 0 0 1 0 0 0"<<endl;
5549  }
5550  // Names of data files
5551  string shortName;
5552  {// Strip path for data file
5553  std::string::size_type idx =prefix.rfind("/");
5554  std::string::size_type idx2=prefix.rfind("\\");
5555  std::string::size_type idx3=prefix.rfind(":");
5556  if(((long)idx2!=(long)string::npos)&&((long)idx2>(long)idx))idx=idx2;
5557  if(((long)idx3!=(long)string::npos)&&((long)idx3>(long)idx))idx=idx3;
5558  if(idx==string::npos)
5559  shortName=prefix;
5560  else
5561  shortName=prefix.substr(idx+1);
5562  }
5563  pcr<<"! File names of data files"<<endl;
5564  pcr<<shortName<<".dat"<<endl;
5565  // Output options...
5566  pcr<<"! Mat Pcr Syo Rpa Sym Sho"<<endl
5567  <<" 1 1 0 -1 0 0 "<<endl;
5568  // Output options... For each pattern
5569  pcr<<"! Ipr Ppl Ioc Ls1 Ls2 Ls3 Prf Ins Hkl Fou Ana"<<endl
5570  <<" 0 0 0 0 0 0 1 10 0 0 1 "<<endl;
5571  // Fixed experimental parameters For each 2-theta pattern :TODO: Check !
5572  int wdt=16;
5573  if(this->GetRadiation().GetRadiationType()==RAD_NEUTRON) wdt=10;
5574  pcr<<"!lambda1 lambda2 Ratio Bkpos Wdt Cthm muR AsyLim Rpolarz -> Patt #1"<<endl
5575  <<this->GetRadiation().GetWavelength()(0)<<" "<<this->GetRadiation().GetWavelength()(0)
5576  << " 0 0 "<<wdt
5577  <<" 0 0 0 0.95"<<endl;
5578  // Refinement parameters - changes are damped !!
5579  pcr<<"!NCY Eps R_at R_an R_pr R_gl"<<endl
5580  <<" 5 0.2 1.0 1.0 1.0 1.0"<<endl;
5581  // Refinement parameters & powder data range, for each 2theta pattern
5582  pcr<<"! Thmin Step Thmax PSD Sent0 -> Patt #1"<<endl
5583  <<" 0 0 0 0 0"<<endl;
5584  // Background points
5585  pcr<<"!2Theta Background for Pattern #1"<<endl;
5586  for(long i=0;i<pBackground->GetInterpPoints().first->numElements();i++)
5587  pcr<<(*(pBackground->GetInterpPoints().first))(i)*RAD2DEG<<" "
5588  <<(*(pBackground->GetInterpPoints().second))(i)<<" 0.0"<<endl;
5589  // Number of refined parameters - just use one for the scale factor !
5590  pcr<<"!"<<endl<<"!"<<endl<<"1 !Number of refined parameters"<<endl;
5591  // Powder data experimental set-up II (refinable parameters)
5592  pcr<<"! Zero Code Sycos Code Sysin Code Lambda Code More -> Patt #1"<<endl;
5593  pcr<<" "<<mXZero*RAD2DEG <<" 0.0 "
5594  <<m2ThetaDisplacement*RAD2DEG <<" 0.0 "
5595  <<m2ThetaTransparency*RAD2DEG <<" 0.0 "
5596  <<"0.000 0.0 0"<<endl;
5597  // PHASE DESCRIPTIONS
5598  for(unsigned int i=0;i<vDiff.size();++i)
5599  {
5600  pcr<<"!-------------------------------------------------------------------------------"<<endl
5601  <<"! Data for PHASE number: "<<i<<" ==> Current R_Bragg for Pattern# 1: 0.00 "<<endl
5602  <<"!-------------------------------------------------------------------------------"<<endl;
5603  //Phase name
5604  pcr<<vDiff[i]->GetCrystal().GetName()<<endl;
5605 
5606  // List all atoms, remove overlapping ones
5607  map<int,exportAtom> vExportAtom;
5608  list<exportBond> vExportBond;
5609  list<exportAngle> vExportAngle;
5610  {
5611  CrystMatrix_REAL minDistTable;
5612  minDistTable=vDiff[i]->GetCrystal().GetMinDistanceTable(-1.);
5613  unsigned long k=0;
5614  // list0 is the full scattering component list with all atoms except dummies,
5615  // and a correct mDynPopCorr
5616  const ScatteringComponentList list0=vDiff[i]->GetCrystal().GetScatteringComponentList();
5617  for(int s=0;s<vDiff[i]->GetCrystal().GetScattererRegistry().GetNb();s++)
5618  {
5619  const ScatteringComponentList list=vDiff[i]->GetCrystal().GetScatt(s).GetScatteringComponentList();
5620 
5621  // If we have a Molecule, remember the names used for the atoms to describe restraints
5622  // We can't use the original atom names as they might not be unique in the crystal
5623  const Molecule *pMol=0;
5624  if(vDiff[i]->GetCrystal().GetScatt(s).GetClassName()=="Molecule")
5625  pMol=dynamic_cast<const Molecule*>(&(vDiff[i]->GetCrystal().GetScatt(s)));
5626  map<const MolAtom*,string> vMolAtomName;
5627 
5628  for(int j=0;j<list.GetNbComponent();j++)
5629  {
5630  if(0==list(j).mpScattPow) continue;//Can this happen ?
5631  bool redundant=false;
5632  for(unsigned long l=0;l<k;++l)
5633  if(abs(minDistTable(l,k))<0.5)
5634  {
5635  map<int,exportAtom>::iterator pos=vExportAtom.find(l);
5636  if(pos!=vExportAtom.end()) pos->second.occMult+=1;
5637  redundant=true;//-1 means dist > 10A
5638  }
5639  if(!redundant)
5640  {
5641  //:TODO: avoid non-alphanumeric characters in name
5642  stringstream name;
5643  name<<list(j).mpScattPow->GetName()<<k+1;
5644  vExportAtom.insert(make_pair(k,exportAtom(name.str(),
5645  list(j).mX,list(j).mY,list(j).mZ,
5646  list(j).mpScattPow->GetBiso(),
5647  list(j).mOccupancy*list0(k).mDynPopCorr,
5648  list(j).mpScattPow)));
5649  if(pMol!=0) vMolAtomName.insert(make_pair(pMol->GetAtomList()[j],name.str()));
5650  }
5651  k++;
5652  }
5653  if(pMol!=0)
5654  {
5655  for(vector<MolBond*>::const_iterator pos=pMol->GetBondList().begin();
5656  pos!=pMol->GetBondList().end();++pos)
5657  {
5658  map<const MolAtom*,string>::const_iterator p1,p2;
5659  p1=vMolAtomName.find(&((*pos)->GetAtom1()));
5660  p2=vMolAtomName.find(&((*pos)->GetAtom2()));
5661  if( (p1!=vMolAtomName.end()) && (p2!=vMolAtomName.end()))
5662  vExportBond.push_back(exportBond(p1->second, p2->second,
5663  (*pos)->GetLength0(),(*pos)->GetLengthSigma()));
5664  }
5665 
5666  for(vector<MolBondAngle*>::const_iterator pos=pMol->GetBondAngleList().begin();
5667  pos!=pMol->GetBondAngleList().end();++pos)
5668  {
5669  map<const MolAtom*,string>::const_iterator p1,p2,p3;
5670  p1=vMolAtomName.find(&((*pos)->GetAtom1()));
5671  p2=vMolAtomName.find(&((*pos)->GetAtom2()));
5672  p3=vMolAtomName.find(&((*pos)->GetAtom3()));
5673  if( (p1!=vMolAtomName.end()) && (p2!=vMolAtomName.end()) && (p3!=vMolAtomName.end()))
5674  vExportAngle.push_back(exportAngle(p1->second, p2->second,p3->second,
5675  (*pos)->GetAngle0(),(*pos)->GetAngleSigma()));
5676  }
5677  }
5678  }
5679  // :TODO: recognize special positions, and move the atoms on them.
5680  // :TODO: list atoms excluded, commented out
5681  }
5682 
5683  // Main control codes line for the phase
5684  //:TODO: extract distance (Dis) and bond angle (Ang) restraints whenever possible
5685  //const ScatteringComponentList *pSC=&(vDiff[i]->GetCrystal().GetScatteringComponentList());
5686  pcr<<"!Nat Dis Ang Jbt Isy Str Furth ATZ Nvk More"<<endl
5687  << vExportAtom.size() <<" "<<vExportBond.size()<<" "<<vExportAngle.size()
5688  <<" 0 0 0 0 1.0 0 1"<<endl;
5689  pcr<<"!Jvi Jdi Hel Sol Mom Ter N_Domains"<<endl
5690  <<" 0 3 0 0 0 0 0"<<endl;
5691  // Contribution to the patterns
5692  pcr<<"!Contributions (0/1) of this phase to the patterns"<<endl
5693  <<" 1"<<endl;
5694  //
5695  pcr<<"!Irf Npr Jtyp Nsp_Ref Ph_Shift for Pattern#"<<i<<endl
5696  <<" 0 0 0 0 0"<<endl;
5697  pcr<<"! Pr1 Pr2 Pr3 Brind. Rmua Rmub Rmuc for Pattern#"<<i<<endl
5698  <<" 1.0 1.0 1.0 1.0 0.0 0.0 0.0"<<endl;
5699 
5700  // Limits for distance calculations
5701  pcr<<"! Max_dst(dist) (angles) Bond-Valence Calc."<<endl
5702  <<" 2.7000 1.5000 0"<<endl;
5703 
5704  // Space group symbol
5705  pcr<<vDiff[i]->GetCrystal().GetSpaceGroup().GetCCTbxSpg().match_tabulated_settings().hermann_mauguin()
5706  <<" <- Space Group Symbol"<<endl;
5707  // Atomic parameters
5708  pcr<<"!Atom Typ X Y Z Biso Occ In Fin N_t Spc / Codes"<<endl;
5709  for(map<int,exportAtom>::const_iterator pos=vExportAtom.begin();pos!=vExportAtom.end();++pos)
5710  {
5711  pcr<<pos->second.name
5712  <<" "<<pos->second.mpScattPow->GetSymbol()<<" "
5713  <<pos->second.x<<" "<<pos->second.y<<" "<<pos->second.z<<" "
5714  <<pos->second.biso<<" "
5715  <<pos->second.occ*pos->second.occMult<<" 0 0 0 0"<<endl
5716  <<" 0 0 0 0 0"<<endl;
5717  }
5718  // POWDER DATA-I: PROFILE PARAMETERS FOR EACH PATTERN
5719  REAL eta0=vDiff[0]->GetProfile().GetPar("Eta0").GetHumanValue();
5720  if(eta0<.01) eta0=.01;
5721  else if(eta0>.99) eta0=.99;
5722  pcr<<"!Scale Shape1 Bov Str1 Str2 Str3 Strain-Model"<<endl
5723  <<" 1.0 "<<eta0
5724  <<" 0.0 0.0 0.0 0.0 0"<<endl
5725  <<" 1.0 0.0 0.0 0.0 0.0 0.0 0"<<endl;
5726 
5727  // :TODO: make sure the profile used corrseponds to pseudo-Voigt first !
5728  pcr<<"! U V W X Y GauSiz LorSiz Size-Model"<<endl;
5729  // :NOTE: we need to separate the 3 next lines, or the same number is generated
5730  // three times when compiled with Visual c++ express 2008 (compiler bug ?)
5731  pcr<<vDiff[i]->GetProfile().GetPar("U").GetHumanValue()<<" ";
5732  pcr<<vDiff[i]->GetProfile().GetPar("V").GetHumanValue()<<" ";
5733  pcr<<vDiff[i]->GetProfile().GetPar("W").GetHumanValue()<<" ";
5734  pcr<< " 0.0 0.0 0.0 0.0 "<<endl
5735  << " 0.0 0.0 0.0 0.0 0.0 0.0 0.0 "<<endl;
5736  // Cell parameters
5737  pcr<<"! a b c alpha beta gamma #Cell Info"<<endl
5738  <<vDiff[i]->GetCrystal().GetLatticePar(0)<<" "
5739  <<vDiff[i]->GetCrystal().GetLatticePar(1)<<" "
5740  <<vDiff[i]->GetCrystal().GetLatticePar(2)<<" "
5741  <<vDiff[i]->GetCrystal().GetLatticePar(3)*RAD2DEG<<" "
5742  <<vDiff[i]->GetCrystal().GetLatticePar(4)*RAD2DEG<<" "
5743  <<vDiff[i]->GetCrystal().GetLatticePar(5)*RAD2DEG<<endl
5744  <<" 0.0 0.0 0.0 0.0 0.0 0.0"<<endl;
5745  pcr<<"! Pref1 Pref2 alpha0 beta0 alpha1 beta1 ?"<<endl
5746  <<" 0.0 0.0 0.0 0.0 0.0 0.0"<<endl
5747  <<" 0.0 0.0 0.0 0.0 0.0 0.0"<<endl;
5748  // ??
5749  //pcr<<"! ??"<<endl
5750  // <<"0.00 0.00 0.00000 0.00"<<endl;
5751  if(vExportBond.size()>0)
5752  {
5753  pcr<<"!Soft distance constraints"<<endl;
5754  for(list<exportBond>::const_iterator pos=vExportBond.begin();pos!=vExportBond.end();++pos)
5755  {
5756  pcr<<pos->at1<<" "<<pos->at2<<" 1 0 0 0 "<<pos->dist<<" "<<pos->sigma<<endl;
5757  }
5758  }
5759  if(vExportBond.size()>0)
5760  {
5761  pcr<<"!Soft angle constraints"<<endl;
5762  for(list<exportAngle>::const_iterator pos=vExportAngle.begin();pos!=vExportAngle.end();++pos)
5763  {
5764  pcr<<pos->at1<<" "<<pos->at2<<" "<<pos->at3<<" 1 1 0 0 0 0 0 0 "
5765  <<pos->ang*RAD2DEG<<" "<<pos->sigma*RAD2DEG<<endl;
5766  }
5767  }
5768  }
5769  pcr.close();
5770 }
5771 
5773 {
5774  this->CalcNbPointUsed();
5776 
5777  TAU_PROFILE("PowderPattern::CalcPowderPattern()","void ()",TAU_DEFAULT);
5778  VFN_DEBUG_ENTRY("PowderPattern::CalcPowderPattern()",3);
5779  if(mPowderPatternComponentRegistry.GetNb()==0)
5780  {
5781  mPowderPatternCalc.resize(mNbPoint);
5783 
5787 
5788  const REAL *p0 = mPowderPatternVariance.data();
5789  REAL *p1=mPowderPatternWeight.data();
5790  for(unsigned long j=0;j<mNbPoint;j++)
5791  {
5792  if(*p0 <=0) {*p1 =0.;}
5793  else *p1 = 1. / *p0;
5794  p0++;p1++;
5795  }
5796  VFN_DEBUG_EXIT("PowderPattern::CalcPowderPattern():no components!",3);
5797  return;
5798  }
5799  TAU_PROFILE_TIMER(timer1,"PowderPattern::CalcPowderPattern1()Calc components","", TAU_FIELD);
5800  TAU_PROFILE_TIMER(timer2,"PowderPattern::CalcPowderPattern2(Add spectrums-scaled)"\
5801  ,"", TAU_FIELD);
5802  TAU_PROFILE_TIMER(timer3,"PowderPattern::CalcPowderPattern2(Add spectrums-backgd1)"\
5803  ,"", TAU_FIELD);
5804  TAU_PROFILE_TIMER(timer4,"PowderPattern::CalcPowderPattern2(Add spectrums-backgd2)"\
5805  ,"", TAU_FIELD);
5806  TAU_PROFILE_TIMER(timer5,"PowderPattern::CalcPowderPattern3(Variance)","", TAU_FIELD);
5807  TAU_PROFILE_START(timer1);
5808  for(int i=0;i<mPowderPatternComponentRegistry.GetNb();i++)
5809  mPowderPatternComponentRegistry.GetObj(i).CalcPowderPattern();
5810  TAU_PROFILE_STOP(timer1);
5811  VFN_DEBUG_MESSAGE("PowderPattern::CalcPowderPattern():Calculated components..",3);
5812  bool b=false;
5814  {
5815  b=true;
5816  }
5817  else
5818  for(int i=0;i<mPowderPatternComponentRegistry.GetNb();i++)
5820  mPowderPatternComponentRegistry.GetObj(i).GetClockPowderPatternCalc())
5821  {
5822  b=true;
5823  break;
5824  }
5825 
5826  if(false==b)
5827  {
5828  VFN_DEBUG_EXIT("PowderPattern::CalcPowderPattern():no need to recalc",3);
5829  return;
5830  }
5831  mPowderPatternCalc.resize(mNbPoint);
5832  int nbBackgd=0;//count number of background phases
5833  for(int i=0;i<mPowderPatternComponentRegistry.GetNb();i++)
5834  {//THIS SHOULD GO FASTER (PRE-FETCHING ARRAY DATA?)
5835  VFN_DEBUG_MESSAGE("PowderPattern::CalcPowderPattern():Adding "<< mPowderPatternComponentRegistry.GetObj(i).GetName(),3);
5836  if(true==mPowderPatternComponentRegistry.GetObj(i).IsScalable())
5837  {
5838  TAU_PROFILE_START(timer2);
5839  if(0==i)
5840  {
5841  const REAL * p1=mPowderPatternComponentRegistry.GetObj(i)
5842  .mPowderPatternCalc.data();
5843  REAL * p0 = mPowderPatternCalc.data();
5844  const REAL s = mScaleFactor(i);
5845  for(unsigned long j=0;j<mNbPointUsed;j++) *p0++ = s * *p1++;
5846  if(!(this->IsBeingRefined())) for(unsigned long j=mNbPointUsed;j<mNbPoint;j++) *p0++ = 0;
5847  }
5848  else
5849  {
5850  const REAL * p1=mPowderPatternComponentRegistry.GetObj(i)
5851  .mPowderPatternCalc.data();
5852  REAL * p0 = mPowderPatternCalc.data();
5853  const REAL s = mScaleFactor(i);
5854  for(unsigned long j=0;j<mNbPointUsed;j++) *p0++ += s * *p1++;
5855  }
5856  TAU_PROFILE_STOP (timer2);
5857  }
5858  else
5859  {// This is a background phase
5860  TAU_PROFILE_START(timer3);
5861  if(0==i)
5862  {
5863  const REAL * p1=mPowderPatternComponentRegistry.GetObj(i)
5864  .mPowderPatternCalc.data();
5865  REAL * p0 = mPowderPatternCalc.data();
5866  for(unsigned long j=0;j<mNbPointUsed;j++) *p0++ = *p1++;
5867  if(!(this->IsBeingRefined())) for(unsigned long j=mNbPointUsed;j<mNbPoint;j++) *p0++ = 0;
5868  }
5869  else
5870  {
5871  const REAL * p1=mPowderPatternComponentRegistry.GetObj(i)
5872  .mPowderPatternCalc.data();
5873  REAL * p0 = mPowderPatternCalc.data();
5874  for(unsigned long j=0;j<mNbPointUsed;j++) *p0++ += *p1++;
5875 
5876  }
5877  TAU_PROFILE_STOP(timer3);
5878  TAU_PROFILE_START(timer4);
5879  // The following is useless if there is only one background phase...
5880  if(0==nbBackgd)
5881  {
5883  REAL *p0 = mPowderPatternBackgroundCalc.data();
5884  const REAL *p1=mPowderPatternComponentRegistry.GetObj(i)
5885  .mPowderPatternCalc.data();
5886  for(unsigned long j=0;j<mNbPointUsed;j++) *p0++ = *p1++;
5887  if(!(this->IsBeingRefined())) for(unsigned long j=mNbPointUsed;j<mNbPoint;j++) *p0++ = 0;
5888  }
5889  else
5890  {
5891  REAL *p0 = mPowderPatternBackgroundCalc.data();
5892  const REAL *p1=mPowderPatternComponentRegistry.GetObj(i)
5893  .mPowderPatternCalc.data();
5894  for(unsigned long j=0;j<mNbPointUsed;j++) *p0++ += *p1++;
5895  }
5896  nbBackgd++;
5897  TAU_PROFILE_STOP(timer4);
5898  }
5899  }
5900  if(0==nbBackgd) mPowderPatternBackgroundCalc.resize(0);//:KLUDGE:
5901 
5902  TAU_PROFILE_START(timer5);
5903  // Calc variance
5904  {
5905  VFN_DEBUG_MESSAGE("PowderPattern::CalcPowderPattern():variance",2);
5908  for(int i=0;i<mPowderPatternComponentRegistry.GetNb();i++)
5909  {
5910  if(mPowderPatternComponentRegistry.GetObj(i).HasPowderPatternCalcVariance())
5911  {
5912  if(0==mPowderPatternComponentRegistry.GetObj(i).GetPowderPatternCalcVariance()
5913  .numElements()) break;
5914 
5915  REAL *p0 = mPowderPatternVariance.data();
5916  const REAL *p1=mPowderPatternComponentRegistry.GetObj(i)
5917  .GetPowderPatternCalcVariance().data();
5918 
5919  if(true==mPowderPatternComponentRegistry.GetObj(i).IsScalable())
5920  {
5921  const REAL s2 = mScaleFactor(i) * mScaleFactor(i);
5922  for(unsigned long j=0;j<mNbPointUsed;j++) *p0++ += s2 * *p1++;
5923  }
5924  else for(unsigned long j=0;j<mNbPointUsed;j++) *p0++ += *p1++;
5925  }
5926  }
5927  REAL *p0 = mPowderPatternWeight.data();
5928  const REAL *p1=mPowderPatternVariance.data();
5929  for(unsigned long j=0;j<mNbPointUsed;j++)
5930  if(*p1 <=0) {*p0++ =0;p1++;}
5931  else *p0++ = 1. / *p1++;
5932  }
5934  TAU_PROFILE_STOP(timer5);
5935  VFN_DEBUG_EXIT("PowderPattern::CalcPowderPattern():End",3);
5936 }
5937 
5938 void PowderPattern::CalcPowderPattern_FullDeriv(std::set<RefinablePar*> &vPar)
5939 {
5940  TAU_PROFILE("PowderPattern::CalcPowderPattern_FullDeriv()","void ()",TAU_DEFAULT);
5941  this->CalcPowderPattern();
5942  mPowderPattern_FullDeriv.clear();
5943  if(mPowderPatternComponentRegistry.GetNb()==0) return;
5944  std::vector<std::map<RefinablePar*,CrystVector_REAL>*> comps;
5945  for(int i=0;i<mPowderPatternComponentRegistry.GetNb();i++)
5946  comps.push_back(&(mPowderPatternComponentRegistry.GetObj(i).GetPowderPattern_FullDeriv(vPar)));
5947 
5948  for(std::set<RefinablePar *>::iterator par=vPar.begin();par!=vPar.end();++par)
5949  {
5950  if(*par==0) continue;
5951  for(int i=0;i<mPowderPatternComponentRegistry.GetNb();i++)
5952  {
5953  if((*par)->GetPointer()==mScaleFactor.data()+i)
5954  {
5955  mPowderPattern_FullDeriv[*par]=mPowderPatternComponentRegistry.GetObj(i).GetPowderPatternCalc();
5956  continue;
5957  }
5958  else
5959  {
5960  if((*(comps[i]))[*par].size()==0) continue;
5961  }
5962  if(true==mPowderPatternComponentRegistry.GetObj(i).IsScalable())
5963  {
5964  if(mPowderPattern_FullDeriv[*par].size()==0)
5965  {
5966  mPowderPattern_FullDeriv[*par].resize(mNbPoint);// :TODO: only use mNbPointUsed
5967  const REAL * p1=(*comps[i])[*par].data();
5968  REAL * p0 = mPowderPattern_FullDeriv[*par].data();
5969  const REAL s = mScaleFactor(i);
5970  for(unsigned long j=0;j<mNbPointUsed;j++) *p0++ = s * *p1++;
5971  for(unsigned long j=mNbPointUsed;j<mNbPoint;j++) *p0++ = 0;// :TODO: only use mNbPointUsed
5972  }
5973  else
5974  {
5975  const REAL * p1=(*comps[i])[*par].data();
5976  REAL * p0 = mPowderPattern_FullDeriv[*par].data();
5977  const REAL s = mScaleFactor(i);
5978  for(unsigned long j=0;j<mNbPointUsed;j++) *p0++ += s * *p1++;
5979  }
5980  }
5981  else
5982  {// This is a background phase
5983  if(mPowderPattern_FullDeriv[*par].size()==0)
5984  {
5985  mPowderPattern_FullDeriv[*par].resize(mNbPoint);// :TODO: only use mNbPointUsed
5986  const REAL * p1=(*comps[i])[*par].data();
5987  REAL * p0 = mPowderPattern_FullDeriv[*par].data();
5988  for(unsigned long j=0;j<mNbPointUsed;j++) *p0++ = *p1++;
5989  for(unsigned long j=mNbPointUsed;j<mNbPoint;j++) *p0++ = 0;
5990  }
5991  else
5992  {
5993  const REAL * p1=(*comps[i])[*par].data();
5994  REAL * p0 = mPowderPattern_FullDeriv[*par].data();
5995  for(unsigned long j=0;j<mNbPointUsed;j++) *p0++ += *p1++;
5996 
5997  }
5998  }
5999  }
6000  }
6001  #if 0
6002  std::map<RefinablePar*, CrystVector_REAL> oldDeriv;
6003  std::vector<const CrystVector_REAL*> v;
6004  int n=0;
6005  for(std::map<RefinablePar*, CrystVector_REAL>::reverse_iterator pos=mPowderPattern_FullDeriv.rbegin();pos!=mPowderPattern_FullDeriv.rend();++pos)
6006  {
6007  if(pos->first==0) continue;
6008  if(pos->second.size()==0) continue;
6009 
6010  const REAL step=pos->first->GetDerivStep();
6011  pos->first->Mutate(step);
6012  this->CalcPowderPattern();
6013  oldDeriv[pos->first]=mPowderPatternCalc;
6014  pos->first->Mutate(-2*step);
6015  this->CalcPowderPattern();
6016  oldDeriv[pos->first]-=mPowderPatternCalc;
6017  oldDeriv[pos->first]/=2*step;
6018  pos->first->Mutate(step);
6019 
6020  v.push_back(&(pos->second));
6021  v.push_back(&(oldDeriv[pos->first]));
6022  cout<<pos->first->GetName()<<":"<<pos->second.size()<<","<<oldDeriv[pos->first].size()<<endl;
6023  if(++n>8) break;
6024  }
6025  cout<<FormatVertVector<REAL>(v,16)<<endl;
6026  exit(0);
6027  #endif
6028 }
6029 
6031 {
6032  this->CalcNbPointUsed();
6034 
6035  this->PrepareIntegratedRfactor();
6036  TAU_PROFILE("PowderPattern::CalcPowderPatternIntegrated()","void ()",TAU_DEFAULT);
6037  VFN_DEBUG_ENTRY("PowderPattern::CalcPowderPatternIntegrated()",4);
6038  if(mPowderPatternComponentRegistry.GetNb()==0)
6039  {
6040  mPowderPatternIntegratedCalc.resize(0);
6042  VFN_DEBUG_EXIT("PowderPattern::CalcPowderPatternIntegrated():no components!",4);
6043  return;
6044  }
6045  TAU_PROFILE_TIMER(timer1,"PowderPattern::CalcPowderPatternIntegrated()1:Calc components",\
6046  "", TAU_FIELD);
6047  TAU_PROFILE_TIMER(timer2,"PowderPattern::CalcPowderPatternIntegrated()2:Add comps-scaled"\
6048  ,"", TAU_FIELD);
6049  TAU_PROFILE_TIMER(timer3,"PowderPattern::CalcPowderPatternIntegrated()2:Add backgd1"\
6050  ,"", TAU_FIELD);
6051  TAU_PROFILE_TIMER(timer4,"PowderPattern::CalcPowderPatternIntegrated()2:Add backgd2"\
6052  ,"", TAU_FIELD);
6053  TAU_PROFILE_TIMER(timer5,"PowderPattern::CalcPowderPatternIntegrated()3:Variance"
6054  ,"", TAU_FIELD);
6055  TAU_PROFILE_START(timer1);
6056  vector< pair<const CrystVector_REAL*,const RefinableObjClock*> > comps;
6057  for(int i=0;i<mPowderPatternComponentRegistry.GetNb();i++)
6058  {
6059  comps.push_back(mPowderPatternComponentRegistry.GetObj(i).
6060  GetPowderPatternIntegratedCalc());
6061  }
6062  TAU_PROFILE_STOP(timer1);
6063  bool b=false;
6065  {
6066  b=true;
6067  }
6068  else
6069  for(vector< pair<const CrystVector_REAL*,const RefinableObjClock*> >::iterator
6070  pos=comps.begin();pos!=comps.end();++pos)
6071  if(mClockPowderPatternCalc < *(pos->second) )
6072  {
6073  b=true;
6074  break;
6075  }
6076 
6077  if(false==b)
6078  {
6079  VFN_DEBUG_EXIT("PowderPattern::CalcPowderPatternIntegrated():no need to recalc",4);
6080  return;
6081  }
6082  VFN_DEBUG_MESSAGE("PowderPattern::CalcPowderPatternIntegrated():Recalc",3);
6083  mPowderPatternIntegratedCalc.resize(mNbIntegrationUsed);
6084  int nbBackgd=0;//count number of background phases
6085  for(int i=0;i<mPowderPatternComponentRegistry.GetNb();i++)
6086  {
6087  VFN_DEBUG_MESSAGE("PowderPattern::CalcPowderPatternIntegrated():Adding "\
6088  << mPowderPatternComponentRegistry.GetObj(i).GetName(),3);
6089  if(true==mPowderPatternComponentRegistry.GetObj(i).IsScalable())
6090  {
6091  TAU_PROFILE_START(timer2);
6092  if(0==i)
6093  {
6094  const REAL * RESTRICT p1= comps[i].first->data();
6095  REAL * RESTRICT p0 = mPowderPatternIntegratedCalc.data();
6096  const REAL s = mScaleFactor(i);
6097  for(unsigned long j=mNbIntegrationUsed;j>0;j--) *p0++ = s * *p1++;
6098  }
6099  else
6100  {
6101  const REAL * RESTRICT p1= comps[i].first->data();
6102  REAL * RESTRICT p0 = mPowderPatternIntegratedCalc.data();
6103  const REAL s = mScaleFactor(i);
6104  for(unsigned long j=mNbIntegrationUsed;j>0;j--) *p0++ += s * *p1++;
6105  }
6106  TAU_PROFILE_STOP (timer2);
6107  }
6108  else
6109  {// This is a background phase
6110  TAU_PROFILE_START(timer3);
6111  if(0==i)
6112  {
6113  const REAL * RESTRICT p1= comps[i].first->data();
6114  REAL * RESTRICT p0 = mPowderPatternIntegratedCalc.data();
6115  for(unsigned long j=mNbIntegrationUsed;j>0;j--) *p0++ = *p1++;
6116  }
6117  else
6118  {
6119  const REAL * RESTRICT p1= comps[i].first->data();
6120  REAL * RESTRICT p0 = mPowderPatternIntegratedCalc.data();
6121  for(unsigned long j=mNbIntegrationUsed;j>0;j--) *p0++ += *p1++;
6122 
6123  }
6124  TAU_PROFILE_STOP(timer3);
6125  TAU_PROFILE_START(timer4);
6126  // The following is useless if there is only one background phase...
6127  if(0==nbBackgd)
6128  {
6129  mPowderPatternBackgroundIntegratedCalc.resize(mNbIntegrationUsed);
6130  const REAL * RESTRICT p1= comps[i].first->data();
6131  REAL * RESTRICT p0 = mPowderPatternBackgroundIntegratedCalc.data();
6132  for(unsigned long j=mNbIntegrationUsed;j>0;j--) *p0++ = *p1++;
6133  }
6134  else
6135  {
6136  const REAL * RESTRICT p1= comps[i].first->data();
6137  REAL * RESTRICT p0 = mPowderPatternBackgroundIntegratedCalc.data();
6138  for(unsigned long j=mNbIntegrationUsed;j>0;j--) *p0++ += *p1++;
6139  }
6140  nbBackgd++;
6141  TAU_PROFILE_STOP(timer4);
6142  }
6143  }
6144  TAU_PROFILE_START(timer5);
6145  if(0==nbBackgd) mPowderPatternBackgroundIntegratedCalc.resize(0);
6146  // Calc variance
6147  bool useCalcVariance=false;
6148  for(int i=0;i<mPowderPatternComponentRegistry.GetNb();i++)
6149  if(mPowderPatternComponentRegistry.GetObj(i).HasPowderPatternCalcVariance())
6150  if(mPowderPatternComponentRegistry.GetObj(i)
6151  .GetPowderPatternIntegratedCalcVariance().first->numElements() !=0)
6152  useCalcVariance=true;
6153  if(useCalcVariance)
6154  {
6155  VFN_DEBUG_MESSAGE("PowderPattern::CalcPowderPatternIntegrated():variance",3);
6156  {
6157  mPowderPatternVarianceIntegrated.resize(mNbIntegrationUsed);
6158  mIntegratedWeight.resize(mNbIntegrationUsed);
6159  const REAL * RESTRICT p1= mIntegratedVarianceObs.data();
6160  REAL * RESTRICT p0 = mPowderPatternVarianceIntegrated.data();
6161  for(unsigned long j=mNbIntegrationUsed;j>0;j--) *p0++ = *p1++;
6162  }
6163  //cout <<"PowderPattern::CalcPowderPatternIntegrated():variance"
6164  // <<"obsvarsum="<<log(mIntegratedVarianceObs.sum());
6165  for(int i=0;i<mPowderPatternComponentRegistry.GetNb();i++)
6166  {
6167  if(mPowderPatternComponentRegistry.GetObj(i).HasPowderPatternCalcVariance())
6168  {
6169  if(0==mPowderPatternComponentRegistry.GetObj(i)
6170  .GetPowderPatternIntegratedCalcVariance().first->numElements()) break;
6171 
6172  const REAL * RESTRICT p1= mPowderPatternComponentRegistry.GetObj(i)
6173  .GetPowderPatternIntegratedCalcVariance().first->data();
6174  //cout <<",sumvar(i)="<<log(mPowderPatternComponentRegistry.GetObj(i)
6175  // .GetPowderPatternIntegratedCalcVariance().first->sum());
6176  REAL * RESTRICT p0 = mPowderPatternVarianceIntegrated.data();
6177 
6178  if(true==mPowderPatternComponentRegistry.GetObj(i).IsScalable())
6179  {
6180  const REAL s2 = mScaleFactor(i) * mScaleFactor(i);
6181  for(unsigned long j=mNbIntegrationUsed;j>0;j--) *p0++ += s2 * *p1++;
6182  }
6183  else for(unsigned long j=mNbIntegrationUsed;j>0;j--) *p0++ += *p1++;
6184 
6185  }
6186  }
6187  //cout <<endl;
6188  REAL *p0 = mIntegratedWeight.data();
6189  const REAL *p1=mPowderPatternVarianceIntegrated.data();
6190  for(unsigned long j=mNbIntegrationUsed;j>0;j--)
6191  if(*p1 <=0) {*p0++ =0;p1++;}
6192  else *p0++ = 1. / *p1++;
6193  }
6194  else mIntegratedWeight.resize(0);
6196  TAU_PROFILE_STOP(timer5);
6197  /*
6198  // Compare-DEBUG ONLY
6199  {
6200  VFN_DEBUG_MESSAGE("PowderPattern::CalcPowderPatternIntegrated():Check",10);
6201  this->CalcPowderPattern();
6202  CrystVector_REAL integr(mNbIntegrationUsed);
6203  for(int k=0;k<mPowderPatternComponentRegistry.GetNb();k++)
6204  {
6205  VFN_DEBUG_MESSAGE("PowderPattern::CalcPowderPatternIntegrated():Check #"<<k,10);
6206  integr=0;
6207  const CrystVector_REAL *v
6208  =&(mPowderPatternComponentRegistry.GetObj(k).GetPowderPatternCalc());
6209  for(unsigned int i=0;i<mNbIntegrationUsed;i++)
6210  {
6211  integr(i)=0;
6212  for(int j=mIntegratedPatternMin(i);j<=mIntegratedPatternMax(i);j++)
6213  {
6214  integr(i) += (*v)(j);
6215  }
6216  }
6217  cout << "Integrated intensities, Component #"<<k<<endl
6218  << FormatVertVector<REAL> (integr,*(comps[k].first))<<endl;
6219  }
6220  }
6221  */
6222  VFN_DEBUG_EXIT("PowderPattern::CalcPowderPatternIntegrated():End",4);
6223 }
6224 
6225 void PowderPattern::CalcPowderPatternIntegrated_FullDeriv(std::set<RefinablePar *> &vPar)
6226 {
6227  TAU_PROFILE("PowderPattern::CalcPowderPatternIntegrated_FullDeriv()","void ()",TAU_DEFAULT);
6229  this->CalcNbPointUsed();
6230 
6231  this->PrepareIntegratedRfactor();
6232  mPowderPatternUsed_FullDeriv.clear();
6233  if(mPowderPatternComponentRegistry.GetNb()==0) return;
6234  std::vector<map<RefinablePar*,CrystVector_REAL>*> comps;
6235  for(int i=0;i<mPowderPatternComponentRegistry.GetNb();i++)
6236  {
6237  comps.push_back(&(mPowderPatternComponentRegistry.GetObj(i).
6238  GetPowderPatternIntegrated_FullDeriv(vPar)));
6239  }
6240  //RefinablePar *scalePar;
6241  //cout<<"PowderPattern::CalcPowderPatternIntegrated_FullDeriv():"<<endl;
6242  mPowderPatternIntegrated_FullDeriv.clear();
6243  for(std::set<RefinablePar *>::iterator par=vPar.begin();par!=vPar.end();++par)
6244  {
6245  if(*par==0) continue; //:TODO: store the calculated (non-derived) pattern here ?
6246  for(int i=0;i<mPowderPatternComponentRegistry.GetNb();i++)
6247  {
6248  if((*par)->GetPointer()==mScaleFactor.data()+i)
6249  {
6250  //scalePar=*par;
6251  mPowderPatternIntegrated_FullDeriv[*par]=*(mPowderPatternComponentRegistry.GetObj(i).GetPowderPatternIntegratedCalc().first);
6252  //cout<<"PowderPattern::CalcPowderPatternIntegrated_FullDeriv():scale #"<<i<<":"<<(*par)->GetName()<<":"<<(*par)->GetPointer()<<":"<<mPowderPatternUsed_FullDeriv[*par]<<endl;
6253  continue;
6254  }
6255  else
6256  {
6257  if((*(comps[i]))[*par].size()==0) continue;
6258  }
6259  if(true==mPowderPatternComponentRegistry.GetObj(i).IsScalable())
6260  {
6261  if(mPowderPatternIntegrated_FullDeriv[*par].size()==0)
6262  {
6263  mPowderPatternIntegrated_FullDeriv[*par].resize(mNbIntegrationUsed);
6264  const REAL * RESTRICT p1= (*(comps[i]))[*par].data();
6265  REAL * RESTRICT p0 = mPowderPatternIntegrated_FullDeriv[*par].data();
6266  const REAL s = mScaleFactor(i);
6267  for(unsigned long j=mNbIntegrationUsed;j>0;j--)
6268  {
6269  #if 1
6270  *p0++ = s * *p1++;
6271  #else
6272  *p0 = s * *p1;
6273  if((j==mNbIntegrationUsed)&&((*par)->GetName()=="Cimetidine_C11_x")) cout<<__FILE__<<":"<<__LINE__<<":"<<*p0<<","<<s<<"*"<<*p1<<endl;
6274  p0++;p1++;
6275  #endif
6276  }
6277  }
6278  else
6279  {
6280  const REAL * RESTRICT p1= (*(comps[i]))[*par].data();
6281  REAL * RESTRICT p0 = mPowderPatternIntegrated_FullDeriv[*par].data();
6282  const REAL s = mScaleFactor(i);
6283  for(unsigned long j=mNbIntegrationUsed;j>0;j--)
6284  {
6285  #if 1
6286  *p0++ += s * *p1++;
6287  #else
6288  *p0 += s * *p1;
6289  if((j==mNbIntegrationUsed)&&((*par)->GetName()=="Cimetidine_C11_x")) cout<<__FILE__<<":"<<__LINE__<<":"<<*p0<<","<<s<<"*"&l