sigmoid.h

Go to the documentation of this file.

#include <cassert>
#include <cstdlib>
 
#include <iostream>
#include <fstream>
 
#include <TF1.h>
#include <TH1.h>
 
using namespace std;
 
namespace DAS::Normalisation {
 
const float threshold = 0.995;
float GetX (TF1 * f)
{
    double m = 0, M = 7000;
    f->GetRange(m,M);
 
    for (float pt = m; pt < M; ++pt)
        if (f->Eval(pt) > threshold)
            return pt;
 
    cerr << "No threshold was found from " << m << " to " << M << ".\nAborting.\n";
    exit(EXIT_FAILURE);
}
 
template<typename STREAM>
[[ deprecated ]]
TF1 * FitTriggerEfficiency
            (int trigger,                   
             TH1 * h,                       
             STREAM& Stream,                
             int from,                      
             const char * options = "NQERS")
{
    assert(from > 0);
    int mbin = h->GetXaxis()->FindBin(  from),
        Mbin = h->GetXaxis()->FindBin(3*from);
    float m = h->GetBinLowEdge(mbin+1), M = h->GetBinLowEdge(Mbin);
    for (int bin = 0; bin <= mbin; ++bin) {
        h->SetBinContent(bin, 0);
        h->SetBinError  (bin, 0);
    }
    for (int bin = Mbin+1; bin <= h->GetNbinsX()+1; ++bin) {
        h->SetBinContent(bin, 0);
        h->SetBinError  (bin, 0);
    }
 
    cout << "Fitting trigger " << trigger << " from " << m << " to " << M << endl;
 
    // function
    TF1 * f = new TF1(Form("sigmoid%d", trigger),
            "[0]+0.5*(1-[0])*(1+erf((x-[1])/[2]))",
            m, M);
 
    f->SetParNames("a", "#mu", "#sigma");
 
    //f->FixParameter(0,0);
    f->SetParameter(0,0.1);
    f->SetParLimits(0,0,0.9);
 
    f->SetParameter(1,1.1*trigger);
    f->SetParLimits(1,m,M);
 
    f->SetParameter(2,10);
    f->SetParLimits(2,1,150);
 
    // fit
    h->Fit(f, options, "", m, M);
 
    float turnon = GetX(f);
 
    Stream << trigger << '\t' << turnon << '\n';
 
    return f;
}
 
template<typename STREAM> float GetTriggerEfficiency
            (int trigger,    
             TH1 * h,        
             STREAM& Stream, 
             float from,     
             float minEff)   
{
    int mbin = h->GetXaxis()->FindBin(  from),
        Mbin = h->GetXaxis()->FindBin(6*from); 
    for (int bin = 0; bin < mbin; ++bin) {
        h->SetBinContent(bin, 0);
        h->SetBinError  (bin, 0);
    }
    for (int bin = Mbin+1; bin <= h->GetNbinsX()+1; ++bin) {
        h->SetBinContent(bin, 0);
        h->SetBinError  (bin, 0);
    }
    for (int bin = mbin; bin <= Mbin; ++bin) {
        float content = h->GetBinContent(bin);
        if (content < minEff) continue;
        float turnon = h->GetBinLowEdge(bin);
        Stream << trigger << '\t' << turnon << '\n';
        return turnon;
    }
    return -1;
}
 
} // end of DAS::Trigger namespace