diff --git a/Rich/RichRecUtils/RichRecUtils/RichCKResolutionFitter.h b/Rich/RichRecUtils/RichRecUtils/RichCKResolutionFitter.h
index 522cc1d83c30912eb132cb16152c107b10a4144b..544f2bf29ab8f47e72dddf897b73fe34ad505e02 100644
--- a/Rich/RichRecUtils/RichRecUtils/RichCKResolutionFitter.h
+++ b/Rich/RichRecUtils/RichRecUtils/RichCKResolutionFitter.h
@@ -5,6 +5,7 @@
#include <string>
#include <vector>
#include <array>
+#include <cstdint>
// ROOT includes
#include <TH1.h>
@@ -33,26 +34,35 @@ namespace Rich
class FitParams final
{
+ private:
+
+ /// Number of Radiators
+ static const std::uint8_t NRads = 2;
+
public: // General fitting parameters
- using RichFitTypes_t = std::array< std::vector<std::string> , 2 >;
+ /// Type for list for fit types for each radiator
+ using RichFitTypes_t = std::array< std::vector<std::string> , NRads >;
- // RICH1 RICH2
+ // RICH1 RICH2
/// Starting resolution guesses
- std::array<double,2> RichStartRes { { 0.0017, 0.0007 } };
+ std::array<double,NRads> RichStartRes { { 0.0017, 0.0007 } };
/// Starting prefit deltas
- std::array<double,2> RichStartDelta { { 0.0025, 0.00105 } };
+ std::array<double,NRads> RichStartDelta { { 0.0025, 0.00105 } };
/// Fit Range Minimums
- std::array<double,2> RichFitMin { { -0.0079, -0.00395 } };
+ std::array<double,NRads> RichFitMin { { - 0.0079, -0.00395 } };
/// Fit range Maxmimums
- std::array<double,2> RichFitMax { { 0.0077, 0.00365 } };
+ std::array<double,NRads> RichFitMax { { 0.0077, 0.00365 } };
+ /// Default Gauss Asymmetry
+ std::array<double,NRads> RichStartAsym { { 0.0000, 0.0000 } };
/// Fit Type(s)
- RichFitTypes_t RichFitTypes { { {"FreeNPol"}, {"FreeNPol"} } };
+ RichFitTypes_t RichFitTypes { { {"AsymNormal:Hyperbolic","AsymNormal:FreeNPol"},
+ {"AsymNormal:Hyperbolic","AsymNormal:FreeNPol"} } };
public: // Free polynominal background fit options
/// Polynominal degree
- std::array<unsigned int,2> RichNPol { { 3 , 3 } };
+ std::array<unsigned int,NRads> RichNPol { { 3 , 3 } };
public: // Fixed polynominal background function
@@ -68,73 +78,203 @@ namespace Rich
private:
+ // polynominal function
+ inline std::string Pol( const int degree,
+ unsigned int iFirstParam ) const noexcept
+ {
+ return (
+ "auto pol = [](double *x, double *p) -> double "
+ "{ "
+ " const unsigned int iFP = "+std::to_string(iFirstParam)+"; "
+ " const unsigned int deg = "+std::to_string(degree)+"; "
+ " auto ret = p[iFP+deg]; "
+ " for ( int i = deg-1; i >= 0; --i ) { ret = (ret*x[0]) + p[iFP+i]; }; "
+ " return ret;"
+ "}; "
+ );
+ }
+
// hyperbolic fnction
- inline std::string Hyperbolic( const unsigned int iFirstParam ) const
+ inline std::string Hyperbolic( const unsigned int iFirstParam ) const noexcept
{
return (
+ "auto hyperbol = [](double *x, double *p) -> double "
+ "{ "
// The fit parameters
- "const double A = p["+std::to_string(iFirstParam )+"]; "
- "const double B = p["+std::to_string(iFirstParam+1)+"]; "
- "const double Cm = p["+std::to_string(iFirstParam+2)+"]; "
- "const double Cp = p["+std::to_string(iFirstParam+3)+"]; "
- "const double D = p["+std::to_string(iFirstParam+4)+"]; "
+ " const double A = p["+std::to_string(iFirstParam )+"]; "
+ " const double B = p["+std::to_string(iFirstParam+1)+"]; "
+ " const double Cm = p["+std::to_string(iFirstParam+2)+"]; "
+ " const double Cp = p["+std::to_string(iFirstParam+3)+"]; "
+ " const double D = p["+std::to_string(iFirstParam+4)+"]; "
// Which C
- "const double C = ( x[0] < D ? Cm : Cp ); "
+ " const double C = ( x[0] < D ? Cm : Cp ); "
// the function
- "const double hyperbol = ( A - ( std::sqrt( B + C*std::pow(x[0]-D,2) ) ) ); "
+ " return ( A - ( std::sqrt( B + C*std::pow(x[0]-D,2) ) ) ); "
+ "}; "
);
}
- // Gaussian function
- inline std::string Gauss( const unsigned int iFirstParam ) const
+ // Asymmetric Normal function
+ inline std::string AsymNormal( const unsigned int iFirstParam ) const noexcept
{
return (
+ "auto normal = [](double *x, double *p) -> double "
+ "{ "
// The fit parameters
- "const double N = p["+std::to_string(iFirstParam )+"]; "
- "const double mean = p["+std::to_string(iFirstParam+1)+"]; "
- "const double sigma = fabs(p["+std::to_string(iFirstParam+2)+"]); "
+ " const double N = p["+std::to_string(iFirstParam )+"]; "
+ " const double mean = p["+std::to_string(iFirstParam+1)+"]; "
+ " const double sigma = fabs(p["+std::to_string(iFirstParam+2)+"]); "
+ " const double alpha = p["+std::to_string(iFirstParam+3)+"]; "
+ " const double sigmaa = sigma * ( 1.0 + ( 0.5 * alpha ) ); "
+ " const double sigmab = sigma * ( 1.0 - ( 0.5 * alpha ) ); "
+ " const double s = ( x[0] > mean ? sigmaa : sigmab ); "
// the function
- "const double gaus = ( sigma > 0 ? N * std::exp( -0.5 * std::pow((x[0]-mean)/sigma,2) ) : 0 ); "
+ " return ( s > 0 ? N * std::exp( -0.5 * std::pow((x[0]-mean)/s,2) ) : 0 ); "
+ "}; "
+ );
+ }
+
+ // Normal function
+ inline std::string Normal( const unsigned int iFirstParam ) const noexcept
+ {
+ return (
+ "auto normal = [](double *x, double *p) -> double "
+ "{ "
+ // The fit parameters
+ " const double N = p["+std::to_string(iFirstParam )+"]; "
+ " const double mean = p["+std::to_string(iFirstParam+1)+"]; "
+ " const double sigma = fabs(p["+std::to_string(iFirstParam+2)+"]); "
+ // the function
+ " return ( sigma > 0 ? N * std::exp( -0.5 * std::pow((x[0]-mean)/sigma,2) ) : 0 ); "
+ "}; "
+ );
+ }
+
+ // Skewed Normal function
+ inline std::string SkewedNormal( const unsigned int iFirstParam ) const noexcept
+ {
+ return (
+ "auto normal = [](double *x, double *p) -> double "
+ "{ "
+ // The fit parameters
+ " const double N = p["+std::to_string(iFirstParam )+"]; "
+ " const double mean = p["+std::to_string(iFirstParam+1)+"]; "
+ " const double sigma = fabs(p["+std::to_string(iFirstParam+2)+"]); "
+ " const double alpha = p["+std::to_string(iFirstParam+3)+"]; "
+ // convert fit parameters to function parameters
+ " const double pi = 3.14159265358979; "
+ " const double delta = alpha / std::sqrt( 1.0 + std::pow(alpha,2) ); "
+ " const double w = sigma / std::sqrt( 1.0 - (2.0*std::pow(delta,2)/pi) ); "
+ " const double eta = mean - (w*delta*std::sqrt(2.0/pi)); "
+ " const double xx = ( x[0] - eta ) / w; "
+ // the function
+ " return N * std::exp(-0.5*xx*xx) * ( 1.0 + std::erf(alpha*xx/std::sqrt(2.0)) ); "
+ "}; "
);
}
public: // Hyperbolic functions
- // RICH1 RICH2
+ // RICH1 RICH2
/// Side band sizes
- std::array<double,2> SideBandSize { { 0.002, 0.001 } };
+ std::array<double,NRads> SideBandSize { { 0.002, 0.001 } };
/// Init values for parameter B
- std::array<double,2> HyperInitB { { 1e8, 1e8 } };
+ std::array<double,NRads> HyperInitB { { 1e8, 1e8 } };
/// Init values for parameter D
- std::array<double,2> HyperInitD { { 7e-4, 5e-4 } };
+ std::array<double,NRads> HyperInitD { { 7e-4, 5e-4 } };
/// lambda function for signal + background
- inline std::string GaussHyperbolFitFunc( const unsigned int iFirstParam ) const
+ inline std::string NormalHyperbolFitFunc( const unsigned int iFirstParam = 0 ) const noexcept
{
return ( "[&](double *x, double *p) { "
- + Gauss(iFirstParam)
+ + Normal(iFirstParam)
+ Hyperbolic(iFirstParam+3)
- + "return hyperbol + gaus; }" );
+ + "return hyperbol(x,p) + normal(x,p); }" );
}
- /// lambda function for background only
- inline std::string HyperbolFitFunc( const unsigned int iFirstParam ) const
+ /// lambda function for signal + background
+ inline std::string AsymNormalHyperbolFitFunc( const unsigned int iFirstParam = 0 ) const noexcept
+ {
+ return ( "[&](double *x, double *p) { "
+ + AsymNormal(iFirstParam)
+ + Hyperbolic(iFirstParam+4)
+ + "return hyperbol(x,p) + normal(x,p); }" );
+ }
+
+ /// lambda function for signal + background
+ inline std::string NormalPolFitFunc( const int degree,
+ const unsigned int iFirstParam = 0 ) const noexcept
+ {
+ return ( "[&](double *x, double *p) { "
+ + Normal(iFirstParam)
+ + Pol(degree,iFirstParam+3)
+ + "return pol(x,p) + normal(x,p); }" );
+ }
+
+ /// lambda function for signal + background
+ inline std::string AsymNormalPolFitFunc( const int degree,
+ const unsigned int iFirstParam = 0 ) const noexcept
+ {
+ return ( "[&](double *x, double *p) { "
+ + AsymNormal(iFirstParam)
+ + Pol(degree,iFirstParam+4)
+ + "return pol(x,p) + normal(x,p); }" );
+ }
+
+ /// lambda function for hyperbol background only
+ inline std::string HyperbolFitFunc( const unsigned int iFirstParam = 0 ) const noexcept
{
return ( "[&](double *x, double *p) { "
+ Hyperbolic(iFirstParam)
- + "return hyperbol; }" );
+ + "return hyperbol(x,p); }" );
+ }
+
+ /// lambda function for pol background only
+ inline std::string PolFitFunc( const int degree,
+ const unsigned int iFirstParam = 0 ) const noexcept
+ {
+ return ( "[&](double *x, double *p) { "
+ + Pol(degree,iFirstParam)
+ + "return pol(x,p); }" );
+ }
+
+ /// lambda function for skewed normal signal
+ inline std::string SkewedNormalFitFunc( const unsigned int iFirstParam = 0 ) const noexcept
+ {
+ return ( "[&](double *x, double *p) { "
+ + SkewedNormal(iFirstParam)
+ + "return normal(x,p); }" );
+ }
+
+ /// lambda function for skewed normal signal + polynominal background
+ inline std::string SkewedNormalPolFitFunc( const unsigned int degree,
+ const unsigned int iFirstParam = 0 ) const noexcept
+ {
+ return ( "[&](double *x, double *p) { "
+ + SkewedNormal(iFirstParam)
+ + Pol(degree,iFirstParam+4)
+ + "return pol(x,p) + normal(x,p); }" );
+ }
+
+ /// lambda function for skewed normal signal + hyperbolic background
+ inline std::string SkewedNormalHyperbolFitFunc( const unsigned int iFirstParam = 0 ) const noexcept
+ {
+ return ( "[&](double *x, double *p) { "
+ + SkewedNormal(iFirstParam)
+ + Hyperbolic(iFirstParam+4)
+ + "return hyperbol(x,p) + normal(x,p); }" );
}
public: // Fit Quality options
/// Maximum abs fitted mean for OK fit
- std::array<double,2> RichMaxMean { { 0.006, 0.003 } };
+ std::array<double,NRads> RichMaxMean { { 0.006, 0.003 } };
/// Maximum abs fitted sigma for OK fit
- std::array<double,2> RichMaxSigma { { 0.010, 0.005 } };
+ std::array<double,NRads> RichMaxSigma { { 0.010, 0.005 } };
/// Maximum shift error for OK fit
- std::array<double,2> MaxShiftError { { 0.001, 0.001 } };
+ std::array<double,NRads> MaxShiftError { { 0.001, 0.001 } };
/// Maximum resolution error for OK fit
- std::array<double,2> MaxSigmaError { { 0.0002, 0.0001 } };
+ std::array<double,NRads> MaxSigmaError { { 0.0002, 0.0001 } };
};
@@ -200,7 +340,7 @@ namespace Rich
}
/// Get the RICH array index from enum
- inline int rIndex( const Rich::RadiatorType rad ) const noexcept
+ inline std::uint8_t rIndex( const Rich::RadiatorType rad ) const noexcept
{
return ( Rich::Rich1Gas == rad ? 0 : 1 );
}
diff --git a/Rich/RichRecUtils/src/RichCKResolutionFitter.cpp b/Rich/RichRecUtils/src/RichCKResolutionFitter.cpp
index 346744edfb29d2a35851b24286c29d6570bc20a4..df294860b2910371b9f6047e512de93d4f107622 100644
--- a/Rich/RichRecUtils/src/RichCKResolutionFitter.cpp
+++ b/Rich/RichRecUtils/src/RichCKResolutionFitter.cpp
@@ -5,6 +5,9 @@
#include <TMath.h>
#include <TFitResult.h>
+// boost
+#include "boost/algorithm/string.hpp"
+
using namespace Rich::Rec;
//=============================================================================
@@ -47,6 +50,18 @@ CKResolutionFitter::fitImp( TH1& hist,
// which RICH ?
const auto rIn = rIndex(rad);
+ // Extract the signal and background forms from the fit type
+ std::vector<std::string> fitOps;
+ boost::split( fitOps, fitType, boost::is_any_of(":") );
+ // Must have exactly two splits signal:background
+ if ( fitOps.size() != 2 )
+ {
+ std::cerr << "Unknown fit type " << fitType << std::endl;
+ return fitRes;
+ }
+ const auto sigType = fitOps[0];
+ const auto bkgType = fitOps[1];
+
// Firstly do a prefit
// delta for the pre-fit
@@ -65,7 +80,8 @@ CKResolutionFitter::fitImp( TH1& hist,
fitMax = xPeak + delta;
// Do the pre-fit
- const std::string preFitFName = "Rad" + std::to_string(rad) + "PreFitF";
+ //std::cout << rad << " Prefit" << std::endl;
+ const auto preFitFName = "Rad" + std::to_string(rad) + "PreFitF";
TF1 preFitF( preFitFName.c_str(), "gaus", fitMin, fitMax );
// Estimate overall scale from max bin
preFitF.SetParameter( 0, yPeak );
@@ -75,6 +91,7 @@ CKResolutionFitter::fitImp( TH1& hist,
preFitF.SetParameter( 2, params().RichStartRes[rIn] );
// run the pre-fit
const auto preFitR = hist.Fit( &preFitF, "MQRS0" );
+ //const auto preFitR = hist.Fit( &preFitF, "MRS0" );
// Fit result status
auto & fitOK = fitRes.fitOK;
@@ -93,8 +110,13 @@ CKResolutionFitter::fitImp( TH1& hist,
fitMin = params().RichFitMin[rIn];
fitMax = params().RichFitMax[rIn];
+ // Number of signal parameters
+ const bool isSkewed = "SkewedNormal" == sigType;
+ const bool isAsym = "AsymNormal" == sigType;
+ const unsigned int nSigParams = 4; // all forms have 4 parms (1 sometimes fixed).
+
// Fit Type
- if ( "FreeNPol" == fitType )
+ if ( "FreeNPol" == bkgType )
{
// Full free polynominal form for background
@@ -114,33 +136,65 @@ CKResolutionFitter::fitImp( TH1& hist,
// Loop over poly orders
for ( unsigned int nPol = 1; nPol < nPolFull+1; ++nPol )
{
+
+ //std::cout << rad << " Pol " << nPol << " Fit" << std::endl;
const auto fFitFName = "Rad" + std::to_string(rad) + "fFitF" + std::to_string(nPol);
- const auto fFitFType = "gaus(0)+pol" + std::to_string(nPol) + "(3)" ;
- TF1 fFitF( fFitFName.c_str(), fFitFType.c_str(), fitMin, fitMax );
-
- const unsigned int nParamsToSet = ( nPol > 1 ? 3 + nPol : 3 );
- for ( unsigned int p = 0; p < nParamsToSet; ++p )
- {
- fFitF.SetParameter(p,lastFitF.GetParameter(p));
+ const auto fFitFType = ( isSkewed ?
+ params().SkewedNormalPolFitFunc(nPol,0) :
+ params().AsymNormalPolFitFunc(nPol,0) );
+ TF1 fFitF( fFitFName.c_str(), fFitFType.c_str(), fitMin, fitMax, nSigParams+nPol+1 );
+
+ // set parameters from last fit
+ // If first iteration, just 3 Gaussian params otherwise number from last full fit
+ const auto nPToSet = ( 1 == nPol ? 3 : nSigParams + nPol );
+ for ( unsigned int i = 0; i < nPToSet; ++i )
+ {
+ const auto p = lastFitF.GetParameter(i);
+ //std::cout << " -> Setting param " << i << " to " << p << std::endl;
+ fFitF.SetParameter(i,p);
}
// set the gaussian parameter names
fFitF.SetParName(0,"Gaus Const");
fFitF.SetParName(1,"Gaus Mean");
fFitF.SetParName(2,"Gaus Sigma");
+ fFitF.SetParName(3,"Gaus Asym");
// Set the background parameter names
- for ( int i = 3; i < fFitF.GetNpar(); ++i )
- { fFitF.SetParName( i, ("Bkg Par "+std::to_string(i-3)).c_str() ); }
-
- // run the fit
- const auto fitR = hist.Fit( &fFitF, minuitOpts.c_str() );
- const auto validFit = fitIsValid( fitR );
-
+ for ( int i = nSigParams; i < fFitF.GetNpar(); ++i )
+ { fFitF.SetParName( i, ("Bkg Par "+std::to_string(i-nSigParams)).c_str() ); }
+
+ // First time, fix par 3
+ fFitF.FixParameter( 3, params().RichStartAsym[rIn] );
+
+ // run the fit with no asymmetry
+ auto fitR = hist.Fit( &fFitF, minuitOpts.c_str() );
+ auto validFit = fitIsValid( fitR );
+ fitOK = validFit && fitIsOK(fFitF,rad);
+
// save the last fit
lastFitF = fFitF;
-
- // Fit OK?
- fitOK = validFit && fitIsOK(fFitF,rad);
+
+ // If doing Asym Gauss fit, run again with that param released
+ if ( fitOK && ( isAsym || isSkewed ) )
+ {
+ fFitF.ReleaseParameter(3);
+ fFitF.SetParameter( 3, params().RichStartAsym[rIn] );
+ // limit asymmetry to at most +-20%
+ if ( isAsym ) { fFitF.SetParLimits(3,-0.2,0.2); }
+ // run the fit
+ fitR = hist.Fit( &fFitF, minuitOpts.c_str() );
+ const auto vFit = fitIsValid( fitR );
+ const auto ok = vFit && fitIsOK(fFitF,rad);
+ if ( ok )
+ {
+ validFit = vFit;
+ fitOK = ok;
+ // save the last fit
+ lastFitF = fFitF;
+ }
+ }
+
+ // Save end fit result ?
if ( fitOK )
{
bestFitF = fFitF;
@@ -165,7 +219,7 @@ CKResolutionFitter::fitImp( TH1& hist,
fitMin, fitMax );
}
- else if ( "FixedNPol" == fitType )
+ else if ( "FixedNPol" == bkgType )
{
// Fixed free polynominal form for background
@@ -175,7 +229,7 @@ CKResolutionFitter::fitImp( TH1& hist,
//const std::string minuitOpts = "MRS0";
const auto fFitFName = "Rad" + std::to_string(rad) + "fFitF";
- const std::string fFitFType = "gaus(0)+[3]*"+params().RichFixedNPolFunc(rad,4);
+ const auto fFitFType = "gaus(0)+[3]*"+params().RichFixedNPolFunc(rad,4);
TF1 fFitF( fFitFName.c_str(), fFitFType.c_str(), fitMin, fitMax );
// Set starting gauss parameters
@@ -207,7 +261,7 @@ CKResolutionFitter::fitImp( TH1& hist,
fitMin, fitMax );
}
- else if ( "Hyperbolic" == fitType )
+ else if ( "Hyperbolic" == bkgType )
{
// Background : y = ( A - sqrt( B + C*(x-D)^2 ) )
@@ -219,28 +273,34 @@ CKResolutionFitter::fitImp( TH1& hist,
// Full signal + Hyperbolic background function
const auto fFitFName = "Rad" + std::to_string(rad) + "fFitF";
- const auto fFitFType = params().GaussHyperbolFitFunc(0);
- TF1 fFitF( fFitFName.c_str(), fFitFType.c_str(), fitMin, fitMax, 8 );
+ const auto fFitFType = ( isSkewed ?
+ params().SkewedNormalHyperbolFitFunc(0) :
+ params().AsymNormalHyperbolFitFunc(0) );
+ TF1 fFitF( fFitFName.c_str(), fFitFType.c_str(), fitMin, fitMax, 9 );
// Set starting gauss parameters
for ( const int i : {0,1,2} )
{ fFitF.SetParameter( i, preFitF.GetParameter(i) ); }
+ // Fix the asym for now
+ fFitF.FixParameter( 3, params().RichStartAsym[rIn] );
+
// Start background level at 0
- for ( const int i : {3,4,5,6,7} )
+ for ( const int i : {4,5,6,7,8} )
{ fFitF.SetParameter( i, 0 ); }
// set the signal parameter names
fFitF.SetParName(0,"Gaus Const");
fFitF.SetParName(1,"Gaus Mean");
fFitF.SetParName(2,"Gaus Sigma");
+ fFitF.SetParName(3,"Gaus Asym");
// background scale factor
- fFitF.SetParName(3,"Hyper A");
- fFitF.SetParName(4,"Hyper B");
- fFitF.SetParName(5,"Hyper C -ve");
- fFitF.SetParName(6,"Hyper C +ve");
- fFitF.SetParName(7,"Hyper D");
+ fFitF.SetParName(4,"Hyper A");
+ fFitF.SetParName(5,"Hyper B");
+ fFitF.SetParName(6,"Hyper C -ve");
+ fFitF.SetParName(7,"Hyper C +ve");
+ fFitF.SetParName(8,"Hyper D");
// Side band fits to get C slope parameters
TF1 sideBmF( (fFitFName+"SideBm").c_str(), "pol1",
@@ -255,10 +315,10 @@ CKResolutionFitter::fitImp( TH1& hist,
const auto CpInit = std::pow( sideBpF.GetParameter(1), 2 );
// Fix some of the background terms
- fFitF.FixParameter(4,params().HyperInitB[rIn]); // B
- fFitF.FixParameter(5,CmInit); // C -ve
- fFitF.FixParameter(6,CpInit); // C +ve
- fFitF.FixParameter(7,params().HyperInitD[rIn]); // D
+ fFitF.FixParameter(5,params().HyperInitB[rIn]); // B
+ fFitF.FixParameter(6,CmInit); // C -ve
+ fFitF.FixParameter(7,CpInit); // C +ve
+ fFitF.FixParameter(8,params().HyperInitD[rIn]); // D
// run the fit with all but the overall background scale set
auto fitR = hist.Fit( &fFitF, minuitOpts.c_str() );
@@ -266,8 +326,8 @@ CKResolutionFitter::fitImp( TH1& hist,
if ( validFitA ) { bestFitF = fFitF; }
// Release C and refit
- fFitF.ReleaseParameter(5);
fFitF.ReleaseParameter(6);
+ fFitF.ReleaseParameter(7);
fitR = hist.Fit( &fFitF, minuitOpts.c_str() );
const bool validFitB = fitIsValid(fitR);
if ( validFitB ) { bestFitF = fFitF; }
@@ -276,7 +336,7 @@ CKResolutionFitter::fitImp( TH1& hist,
bool validFitC = false;
//if ( validFitB )
{
- fFitF.ReleaseParameter(7);
+ fFitF.ReleaseParameter(8);
fitR = hist.Fit( &fFitF, minuitOpts.c_str() );
validFitC = fitIsValid(fitR) && fitIsOK(fFitF,rad);
if ( validFitC ) { bestFitF = fFitF; }
@@ -286,14 +346,33 @@ CKResolutionFitter::fitImp( TH1& hist,
bool validFitD = false;
//if ( validFitC )
{
- fFitF.ReleaseParameter(4);
+ fFitF.ReleaseParameter(5);
fitR = hist.Fit( &fFitF, minuitOpts.c_str() );
validFitD = fitIsValid(fitR) && fitIsOK(fFitF,rad);
if ( validFitD ) { bestFitF = fFitF; }
}
+ // Release asym/skew param
+ if ( validFitD && ( isAsym || isSkewed ) )
+ {
+ fFitF.ReleaseParameter(3);
+ fitR = hist.Fit( &fFitF, minuitOpts.c_str() );
+ validFitD = fitIsValid(fitR) && fitIsOK(fFitF,rad);
+ if ( validFitD ) { bestFitF = fFitF; }
+ }
+
// final OK status based on best fit
- fitOK = ( validFitD || validFitC );
+ fitOK = ( "AsymNormal" == sigType ? validFitD : // for asym fit all must be ok
+ validFitD || validFitC ); // otherwise use either of last 2
+
+ // Check B is not too small compared to C
+ if ( fitOK )
+ {
+ const auto cm = fabs( fFitF.GetParameter(6) );
+ const auto cp = fabs( fFitF.GetParameter(7) );
+ const auto b = fabs( fFitF.GetParameter(5) );
+ fitOK = ( cm < b*1e7 && cp < b*1e7 );
+ }
// Setup the background function
fitRes.bkgFitFunc = TF1( "Background",
@@ -303,17 +382,17 @@ CKResolutionFitter::fitImp( TH1& hist,
}
else
{
- std::cerr << "Unknown fit type " << fitType << std::endl;
+ std::cerr << "Unknown background type " << bkgType << std::endl;
}
// Set the background function parameters
- for ( int i = 3; i < bestFitF.GetNpar(); ++i )
+ for ( int i = nSigParams; i < bestFitF.GetNpar(); ++i )
{
- fitRes.bkgFitFunc.SetParameter( i-3, bestFitF.GetParameter(i) );
- fitRes.bkgFitFunc.SetParError ( i-3, bestFitF.GetParError(i) );
- fitRes.bkgFitFunc.SetParName ( i-3, bestFitF.GetParName(i) );
+ fitRes.bkgFitFunc.SetParameter( i-nSigParams, bestFitF.GetParameter(i) );
+ fitRes.bkgFitFunc.SetParError ( i-nSigParams, bestFitF.GetParError(i) );
+ fitRes.bkgFitFunc.SetParName ( i-nSigParams, bestFitF.GetParName(i) );
}
-
+
// Final sanity check on the fit results
// limits on fit sigma and mean, to detect complete failures such as when there
// is no Cherenkov signal at all