Skip to content

Commit

Permalink
Refs #5772 Fixed some errors
Browse files Browse the repository at this point in the history
  • Loading branch information
@RuthFromDuluth
RuthFromDuluth committed Oct 9, 2012
1 parent 805a719 commit 60f1a1a
Show file tree
Hide file tree
Showing 3 changed files with 67 additions and 43 deletions.
12 changes: 8 additions & 4 deletions Code/Mantid/Framework/Crystal/inc/MantidCrystal/IntegratePeakTimeSlices.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -139,7 +139,7 @@ namespace Crystal
double col, double &Varx, double &Vary, double &Varxy,
std::vector<double>&ParameterValues);

bool IsEnoughData(const double *ParameterValues, Kernel::Logger& g_log);
bool IsEnoughData(const double *ParameterValues, Kernel::Logger& );

double getNewRCRadius();

Expand Down Expand Up @@ -186,14 +186,16 @@ namespace Crystal

bool CalcVariances( );



double StatBaseVals( int index)
{
return StatBase[index];
}

double CalcISAWIntensity( const double* params) const;
double CalcISAWIntensity( const double* params) ;

double CalcISAWIntensityVariance( const double* params, const double* errs, double chiSqOvDOF) const;
double CalcISAWIntensityVariance( const double* params, const double* errs, double chiSqOvDOF) ;

double CalcSampleIntensityMultiplier( const double* params) const;

Expand All @@ -214,6 +216,7 @@ namespace Crystal
std::vector<double> StatBase;

double EdgeX,EdgeY;
double lastISAWIntensity, lastISAWVariance;
bool CalcVariance;
bool case4;//if true result of successful merge of dir =1 chan=0 and chan=1
double back_calc,
Expand All @@ -238,7 +241,8 @@ namespace Crystal
CalcVariance = true;
CellWidth = CellHeight = 0;
currentRadius = -1;

lastISAWIntensity = -1;
lastISAWVariance = -1;
currentPosition = Kernel::V3D();
HalfWidthAtHalfHeightRadius = -1;
case4 = false;
Expand Down
90 changes: 55 additions & 35 deletions Code/Mantid/Framework/Crystal/src/IntegratePeakTimeSlices.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -236,6 +236,8 @@ namespace Mantid
this->EdgeY = handler.EdgeY;
this->CalcVariance = handler.CalcVariance;

this->lastISAWIntensity=handler.lastISAWIntensity;
this->lastISAWVariance = handler.lastISAWIntensity;
this->back_calc = handler.back_calc;
this->Intensity_calc = handler.Intensity_calc;
this->row_calc = handler.row_calc;
Expand Down Expand Up @@ -280,6 +282,7 @@ namespace Mantid
double Variance = StatBase[IVariance];
double TotBoundaryIntensities= StatBase[ITotBoundary ];
double TotBoundaryVariances = StatBase[IVarBoundary];

double nBoundaryCells= StatBase[INBoundary];
double back = TotBoundaryIntensities/nBoundaryCells;
double backVar= TotBoundaryVariances/nBoundaryCells/nBoundaryCells;
Expand All @@ -299,6 +302,7 @@ namespace Mantid
str<< max<double>(0.0, back_calc-(1+2*relError)*sqrt(backVar))<<"<Background<"<<(back+(1+2*relError)*sqrt(backVar))
<<","<< max<double>(0.0, Intensity_calc-(1+2*relError)*sqrt(IntensVar))<<
"<Intensity<"<<TotIntensity-ncells*back+(1+2*relError)*sqrt(IntensVar);

double min =max<double>(0.0, back_calc-(1+2*relError)*sqrt(backVar));
double maxx =back+(1+2*relError)*sqrt(backVar);
Bounds.push_back( pair<double,double>(min,maxx ));
Expand Down Expand Up @@ -803,6 +807,9 @@ namespace Mantid
lastRow = Row0;
lastCol = Col0;
lastAttributeList = origAttributeList;
if( TabWS->rowCount() >0 )
LastTableRow = 0;

}
else if( Chan+dir*chan <0 || Chan+dir*chan >= (int)X.size())
done = true;
Expand All @@ -827,9 +834,7 @@ namespace Mantid
MatrixWorkspace_sptr Data = WorkspaceFactory::Instance().create(
std::string("Workspace2D"), 3, NN,NN);//neighbors.size(), neighbors.size());

// sprintf(logInfo, string(
// " A:chan= %d time=%7.2f Radius=%7.3f row= %5.2f col=%5.2f \n").c_str(),
// xchan, time, Radius, lastRow, lastCol);//std::string(logInfo)

g_log.debug()<<" A:chan="<< xchan<<" time="<<time<<" Radius="<<Radius
<<"row= "<<lastRow<<" col="<<lastCol<<std::endl;

Expand Down Expand Up @@ -951,10 +956,10 @@ namespace Mantid
LastTableRow =UpdateOutputWS(TabWS, dir, (chanMin+chanMax)/2.0, params, errs,
names, chisqOverDOF, AttributeValues->time, spec_idList);

if( lastAttributeList->CellHeight > 0 && lastAttributeList->StatBase.size()>=NAttributes)
if( lastAttributeList->lastISAWVariance > 0 )
{
TotIntensity -=lastAttributeList->StatBaseVals(IIntensities);
TotVariance -= lastAttributeList->StatBaseVals(IVariance);
TotIntensity -=lastAttributeList->lastISAWIntensity;
TotVariance -= lastAttributeList->lastISAWVariance;
}

double TotSliceIntensity = AttributeValues->StatBaseVals(IIntensities);
Expand Down Expand Up @@ -1716,7 +1721,7 @@ namespace Mantid

double x = params[IIntensity] / (AttributeValues->StatBaseVals(IIntensities) - params[Ibk] * ncells);

if ((x < .8 || x > 1.25)&& !EdgePeak)// The fitted intensity should be close to tot intensity - background
if ((x < .7 || x > 1.45)&& !EdgePeak)// The fitted intensity should be close to tot intensity - background
{
g_log.debug()<< " Bad Slice. Fitted Intensity & Observed Intensity(-back) too different. ratio="
<<x<<std::endl;
Expand Down Expand Up @@ -1761,14 +1766,13 @@ namespace Mantid

GoodNums = true;

// std::cout<<"xxxx"<<params[Ibk]<<","<<(params[Ibk] < -.001)<<","<<(params[IIntensity]<0)<<std::endl;

if (params[Ibk] < -.001)
GoodNums = false;

if (params[IIntensity] < 0)
GoodNums = false;

//double sqrChi = SQRT(chisqOverDOF);

double IsawIntensity= AttributeValues->CalcISAWIntensity( params.data());
double IsawVariance = AttributeValues->CalcISAWIntensityVariance(params.data(), errs.data(),chisqOverDOF);
Expand Down Expand Up @@ -1895,11 +1899,15 @@ namespace Mantid

TabWS->getRef<double> (std::string("TotIntensity"), TableRow) = AttributeValues->StatBaseVals(IIntensities);
TabWS->getRef<double> (std::string("BackgroundError"), TableRow) = errs[Ibk] * SQRT(chisq);
TabWS->getRef<double> (std::string("ISAWIntensity"), TableRow) = AttributeValues->StatBaseVals(IIntensities)
- params[Ibk] * ncells;
TabWS->getRef<double> (std::string("ISAWIntensity"), TableRow) = AttributeValues->CalcISAWIntensity(params.data());
//AttributeValues->StatBaseVals(IIntensities)
// - params[Ibk] * ncells;

TabWS->getRef<double> (std::string("ISAWIntensityError"), TableRow) = CalculateIsawIntegrateError(
params[Ibk], errs[Ibk], chisq, AttributeValues->StatBaseVals(IVariance), ncells);
TabWS->getRef<double> (std::string("ISAWIntensityError"), TableRow) = sqrt(AttributeValues->
CalcISAWIntensityVariance(params.data(),errs.data(), Chisq));

//CalculateIsawIntegrateError(
//params[Ibk], errs[Ibk], chisq, AttributeValues->StatBaseVals(IVariance), ncells);

TabWS->getRef<double> (std::string("Time"), TableRow) = time;

Expand Down Expand Up @@ -1932,28 +1940,31 @@ namespace Mantid
double const chisqdivDOF,
const int ncells)
{
int Ibk = find("Background", names);
UNUSED_ARG( TotSliceIntensity);
UNUSED_ARG( TotSliceVariance );
UNUSED_ARG( ncells);
// int Ibk = find("Background", names);
double err =0;
// double intensty=0; //for debugging only

if( !EdgePeak )
{
err = CalculateIsawIntegrateError(params[Ibk], errs[Ibk], chisqdivDOF, TotSliceVariance, ncells);
// intensty =TotSliceIntensity - params[IBACK] * ncells;
TotIntensity += TotSliceIntensity - params[IBACK] * ncells;
err = AttributeValues->CalcISAWIntensityVariance(params.data(),errs.data(), chisqdivDOF);
//CalculateIsawIntegrateError(params[Ibk], errs[Ibk], chisqdivDOF, TotSliceVariance, ncells);

TotIntensity += AttributeValues->CalcISAWIntensity( params.data());
//TotSliceIntensity - params[IBACK] * ncells;

TotVariance += err * err;
TotVariance += err ;

}else
{
int IInt = find("Intensity", names);
// intensty = params[IInt];

TotIntensity += params[IInt];
err = errs[IInt];
TotVariance += err * err;

}
// std::cout<<EdgePeak<<"("<<intensty<<","<<AttributeValues->CalcISAWIntensity(params.data())<<")("
// <<err*err<<","<< AttributeValues->CalcISAWIntensityVariance(params.data(),errs.data(), chisqdivDOF)<<")"<<std::endl;

}

Expand All @@ -1972,24 +1983,31 @@ namespace Mantid
return false;
}

bool DataModeHandler::IsEnoughData( const double *ParameterValues, Kernel::Logger& g_log )
bool DataModeHandler::IsEnoughData( const double *ParameterValues, Kernel::Logger& )
{


// Check if flat
/* double Varx,Vary, Cov; No Good. Eliminated too many peaks
double Varx,Vary, Cov;

std::vector<double>PP(ParameterValues, ParameterValues+IVXY);
CalcVariancesFromData( ParameterValues[IBACK],ParameterValues[IXMEAN],ParameterValues[IYMEAN],
Varx,Cov,Vary, PP);
double ncells = (int)StatBase[IIntensities];
if( ncells <=0)
return false;
double meanx = StatBase[ISSIx]/ncells;
double meany= StatBase[ISSIy]/ncells;

if( Varx <4 || Vary < 4)//<--- This one especially
CalcVariancesFromData( 0, meanx,meany, Varx,Cov,Vary, PP);

if( Varx <.6 || Vary <.6) //All data essentially the same.
return false;

if( Cov*Cov < .75*Varx*Vary)//All data on a line
if( Cov*Cov > .90*Varx*Vary)//All data on a obtuse line
return false;
*/
double VIx0_num = StatBase[ISSIxx] - 2 * ParameterValues[IXMEAN] * StatBase[ISSIx]



/* double VIx0_num = StatBase[ISSIxx] - 2 * ParameterValues[IXMEAN] * StatBase[ISSIx]
+ ParameterValues[IXMEAN] * ParameterValues[IXMEAN] * StatBase[IIntensities];
double VIy0_num = StatBase[ISSIyy] - 2 * ParameterValues[IYMEAN] * StatBase[ISSIy]
Expand Down Expand Up @@ -2032,12 +2050,12 @@ namespace Mantid
g_log.debug()<<" Not Enough Data, Peak Heigh "<<Z <<" is too low"<<std::endl;
return false;
}

*/
return true;

}

double DataModeHandler::CalcISAWIntensity(const double* params)const
double DataModeHandler::CalcISAWIntensity(const double* params)
{

double ExperimentalIntensity = StatBase[IIntensities]-
Expand All @@ -2049,11 +2067,12 @@ namespace Mantid
double alpha =(float)( 0+.5*( r-1 ) );
alpha = std::min<double>( 1.0f , alpha );

return ExperimentalIntensity*r*( 1-alpha )+ alpha * FitIntensity;
lastISAWIntensity = ExperimentalIntensity*r*( 1-alpha )+ alpha * FitIntensity;
return lastISAWIntensity;

}

double DataModeHandler::CalcISAWIntensityVariance( const double* params,const double* errs, double chiSqOvDOF)const
double DataModeHandler::CalcISAWIntensityVariance( const double* params,const double* errs, double chiSqOvDOF)
{

int ncells = (int)StatBase[ISS1];
Expand All @@ -2073,7 +2092,8 @@ namespace Mantid
double alpha =(float)( 0+.5*( r - 1 ));
alpha = std::min<double>( 1.0f , alpha );

return ExperimVar*r*r*( 1 - alpha ) + alpha * FitVar;
lastISAWVariance = ExperimVar*r*r*( 1 - alpha ) + alpha * FitVar;
return lastISAWVariance;
}

double DataModeHandler::CalcSampleIntensityMultiplier( const double* params)const
Expand Down
8 changes: 4 additions & 4 deletions Code/Mantid/Framework/Crystal/test/PeakIntegrationTest.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -202,12 +202,12 @@ class PeakIntegrationTest : public CxxTest::TestSuite

double intensity =peak.getIntensity();
double sigIntensity =peak.getSigmaIntensity();
// std::cout<<"Peak Intens,sig,slice="<<intensity;
// std::cout<<","<<sigIntensity<<","<<slices<<std::endl;
double intensity0=3141;
//std::cout<<"Peak Intens,sig,slice="<<intensity;
//std::cout<<","<<sigIntensity<<","<<IC<<std::endl;
double intensity0=2509;
TS_ASSERT_DELTA(intensity,intensity0, 400.0);

double sig0=115;
double sig0=100;


TS_ASSERT_DELTA( sigIntensity,sig0, 25.0);
Expand Down

0 comments on commit 60f1a1a

Please sign in to comment.