ResponseFit

Description

Performs a double-sided Crystal-Ball fit of the given response step by step. Three methods are provided to fit the gaussian core and each of the power- law tails.

Inheritance diagram for ResponseFit:

Collaboration diagram for ResponseFit:

Public Types
enum	{   N = 0, MU, SIGMA, KL,   NL, KR, NR, NPARS }
enum	Status { failed = 0b0000, core = 0b0001, LHStail = 0b0010, RHStail = 0b0100 }

Public Member Functions
	ResponseFit (const unique_ptr< TH1 > &h, ostream &cout)
void	gausCore ()
void	SSCB ()
void	DSCB ()
bool	good () const override

Private Member Functions
void	Write (const char *) override

Additional Inherited Members
Public Attributes inherited from AbstractFit
const std::unique_ptr< TH1 > &	h
std::uint32_t	status
std::ostream &	cout
double *	p
double *	e
std::unique_ptr< TF1 >	f
std::optional< double >	chi2ndf
std::optional< double >	chi2ndfErr
std::pair< float, float >	interval
Protected Member Functions inherited from AbstractFit
	AbstractFit (const unique_ptr< TH1 > &h, ostream &cout, int npars)
virtual	~AbstractFit ()
void	fit (std::pair< float, float >, const char *)

Member Enumeration Documentation

anonymous enum

anonymous enum

Enumerator
N
MU
SIGMA
KL
NL
KR
NR
NPARS

         {
        N = 0,
        MU, SIGMA,
        KL, NL, KR, NR,
        NPARS
    }; 

Status

enum Status

Enumerator
failed
core
LHStail
RHStail

                {
        failed     = 0b0000,
        core       = 0b0001,
        LHStail    = 0b0010,
        RHStail    = 0b0100,
    };

Constructor & Destructor Documentation

ResponseFit()

ResponseFit ( const unique_ptr< TH1 > &  h, ostream &  cout  )

inline

Preparing for fit, setting ranges. Assuming response is centered around 1

Todo:

tune range?

                                                          :
        AbstractFit(h, cout, NPARS)
    {
        // Calculate 1st & 2nd derivatives of the log(h)
        unique_ptr<TH1> logd  = DerivativeFivePointStencil(h, safelog), // 1st derivative of log(h)
                        logdd = DerivativeFivePointStencil(logd      ); // 2nd derivative of log(h)
        logd->Smooth();
        logd->Write("logd");
        logdd->Write("logdd");
 
        double mu      = h->GetMean(),
               sigma   = h->GetStdDev();
 
        float rmin = 0.65, rmax = 1.5; 
        interval = findDLogExtrema(logd, mu, {rmin,rmax});
 
        // Initial estimate on the tail parameters
        static auto LHS = make_unique<TF1>("tail", "[0]*exp( [1]*x)", 0., rmin),
                    RHS = make_unique<TF1>("tail", "[0]*exp(-[1]*x)", rmax, 2.);
        h->Fit(LHS.get(), "N0QR");
        h->Fit(RHS.get(), "N0QR");
 
        f = make_unique<TF1>("f", DAS::DoubleCrystalBall::Distribution(),
                                interval.first, interval.second, NPARS);
        p[N    ] = 1.;                      e[N    ] = 0.001;
        p[MU   ] = mu;                      e[MU   ] = h->GetMeanError();
        p[SIGMA] = sigma;                   e[SIGMA] = h->GetStdDevError();
        p[KL   ] = interval.first;          e[KL   ] = sigma;
        p[NL   ] = LHS->GetParameter(1);    e[NL   ] = LHS->GetParError(1);
        p[KR   ] = interval.second;         e[KR   ] = sigma;
        p[NR   ] = RHS->GetParameter(1);    e[NR   ] = RHS->GetParError(1);
        f->SetParameters(p);
        f->SetParErrors (e);
 
        f->SetParName(N    ,  "N"     );
        f->SetParName(MU   ,  "#mu"   );
        f->SetParName(SIGMA,  "#sigma");
        f->SetParName(KL   ,  "k_{L}" );
        f->SetParName(NL   ,  "N_{L}" );
        f->SetParName(KR   ,  "k_{R}" );
        f->SetParName(NR   ,  "N_{R}" );
 
        f->SetNpx(2000); // Increase points used to draw TF1 object when saving
 
        // The Integral of a DSCB function using the DSCB parameters from fit.
        auto IDSCB = make_unique<TF1>("DSCB_Integral", DAS::DoubleCrystalBall::Integral(), 0., 2., NPARS-1);
        IDSCB->SetParameters(p);
        IDSCB->SetLineColor(f->GetLineColor());
        IDSCB->SetLineStyle(f->GetLineStyle());
 
        IDSCB->Write("DSCB_Integral");
    }

Member Function Documentation

DSCB()

void DSCB

(

)

Crystal-Ball fit with tail on the RHS

The fit performance is scanned for all possible RHS boundaries.

{
    f->SetParLimits(N, min(0.9,p[N]), max(1.1,p[N]));
    if ((status & Status::LHStail) == Status::LHStail)
        f->SetParLimits(KL, p[KL]-0.5*p[SIGMA], p[KL]+0.5*p[SIGMA]);
    else {
        f->FixParameter(KL, p[KL]);
        f->FixParameter(NL, p[NL]);
    }
 
    auto lastChi2ndf = chi2ndf;
    float L = interval.first;
    int BinL = h->FindBin(L);
    for (int BinR = h->FindBin(interval.second); BinR <= h->GetNbinsX(); ++BinR) {
        if (h->GetBinContent(BinR) == 0) continue;
        if (BinR - BinL < f->GetNumberFreeParameters()) continue; // Fit range should include more bins than the number of params we want to fit
        float R = h->GetBinLowEdge(BinR+1);
        if (R <= p[MU]) continue; // if true, range doesnt make sense, skip fit
 
        f->SetParLimits(KR, p[MU], R);
        f->SetParLimits(NR, 1.001, min(1e3,p[NR]*2));
        fit({L, R}, __func__); // Fit response distribution in given range
    }
 
    if (!chi2ndf || *chi2ndf == *lastChi2ndf) return;
    status |= Status::RHStail;
    Write(__func__);
}

gausCore()

void gausCore

(

)

Get Gaussian core parameters.

{
    f->SetParLimits(N    , 0.9, 1.1);
    f->SetParLimits(MU   , p[MU]-0.5*p[SIGMA], p[MU]+0.5*p[SIGMA]);
    f->SetParLimits(SIGMA, 0.5*p[SIGMA], 2*p[SIGMA]);
    f->FixParameter(KL, p[KL]);
    f->FixParameter(NL, p[NL]);
    f->FixParameter(KR, p[KR]);
    f->FixParameter(NR, p[NR]);
 
    fit(interval, __func__);
 
    if (!chi2ndf) return;
    status |= Status::core;
    Write(__func__);
}

good()

bool good ( ) const

inlineoverridevirtual

Reimplemented from AbstractFit.

{ return AbstractFit::good() && *chi2ndf < 1000; }

SSCB()

void SSCB

(

)

Crystal-Ball fit with tail on the LHS

The fit performance is scanned for all possible LHS boundaries.

{
    f->SetParLimits(N    , min(0.9,p[N]), max(1.1,p[N]));
    f->SetParLimits(MU   , p[MU]-0.1*p[SIGMA], p[MU]+0.1*p[SIGMA]);
    f->SetParLimits(SIGMA, 0.5*p[SIGMA], 1.5*p[SIGMA]);
 
    auto lastChi2ndf = chi2ndf;
    float R = interval.second;
    int BinR = h->FindBin(R);
    for (int BinL = h->FindBin(interval.first); BinL > 0; --BinL) {
        if (h->GetBinContent(BinL) == 0) continue;
        if (BinR - BinL < f->GetNumberFreeParameters()) continue; // Fit range should include more bins than the number of params we want to fit
        float L = h->GetBinLowEdge(BinL);
        if (L >= p[MU]) continue; // if true, range doesnt make sense, skip fit
 
        f->SetParLimits(KL, L, p[MU]);
        f->SetParLimits(NL, 1.001, min(1e3,p[NL]*2));
        fit({L, R}, __func__); // Fit response distribution in given range
    }
 
    if (!chi2ndf || *chi2ndf >= *lastChi2ndf) return;
    status |= Status::LHStail;
    Write(__func__);
}

Write()

void Write ( const char *  name )

overrideprivatevirtual

Write best current estimate in current directory.

Four color codes are used to categorize fits on their final state (after all 3 fits are performed):

• Red for a failed fit (all fit stages failed, or fit quality did not pass some criteria, e.g., chi2/ndf)
• Magenta a successful Gaussian core fit (fits at stage 2 and 3 failed)
• Orange a successful SSCB fit (stage 1 fit could either succeed or fail, and stage 3 failed)
• Green a successful DSCB fit (stage 1 and 2 fits could have either failed or succeded)

Reimplemented from AbstractFit.

{
    bool goodCore = (status & Status::core   ) == Status::core   ;
    bool goodLHS  = (status & Status::LHStail) == Status::LHStail;
    bool goodRHS  = (status & Status::RHStail) == Status::RHStail;
    if      (goodLHS &&  goodRHS) f->SetLineColor(kGreen+1);
    else if (goodLHS xor goodRHS) f->SetLineColor(kOrange+1);
    else if (goodCore           ) f->SetLineColor(kMagenta+1);
    else                          f->SetLineColor(kRed+1);
 
    AbstractFit::Write(name);
 
    auto dlog = make_unique<TF1>("DLog", DAS::DoubleCrystalBall::DLog(), 0., 2., NPARS-1);
    dlog->SetParameters(p);
    dlog->SetLineColor(f->GetLineColor());
    dlog->SetLineStyle(f->GetLineStyle());
 
    static TString dlogName = "DLog";
    dlog->Write(dlogName + name);
}

---

The documentation for this struct was generated from the following file:

• /builds/cms-analysis/general/DasAnalysisSystem/Core/Installer/Core/JEC/bin/fitJetResponse.cc