Belle II Software development
SliceFit.h
1/**************************************************************************
2 * basf2 (Belle II Analysis Software Framework) *
3 * Author: The Belle II Collaboration *
4 * *
5 * See git log for contributors and copyright holders. *
6 * This file is licensed under LGPL-3.0, see LICENSE.md. *
7 **************************************************************************/
8#include "TH2D.h"
9#include "TH2.h"
10#include "TH1D.h"
11#include "TF1.h"
12#include "TDirectory.h"
13#include "TError.h"
14#include "TROOT.h"
15#include "TString.h"
16namespace Belle2 {
25 class SliceFit {
26 public:
30 static TH1D* doSliceFitY(TH2D* h2, int minHitCut = 0)
31 {
32 gPrintViaErrorHandler = kTRUE;
33 gErrorIgnoreLevel = 3001;
34 TString hist_name = h2->GetName();
35 double ub = h2->GetYaxis()->GetXmax();
36 double lb = h2->GetYaxis()->GetXmin();
37 B2DEBUG(199, "Axis: " << lb << " " << ub);
38 if ((h2->GetEntries() / h2->GetNbinsX()) < 30) {
39 B2WARNING("Low statictic: " << h2->GetEntries() << " Hits");
40 h2->Rebin2D(2, 2, hist_name);
41 }
42
43 B2DEBUG(199, "Slice fit for histo " << hist_name);
44 B2DEBUG(199, "Number of entries: " << h2->GetEntries());
45 TF1* g1 = new TF1("g1", "gaus", lb, ub);
46 h2->FitSlicesY(0, 0, -1, minHitCut);
47
48
49 TString m_name = hist_name + "_1";
50 TH1D* hm = (TH1D*)gDirectory->Get(m_name)->Clone("hm");
51 if (!hm) return 0;
52
53 B2DEBUG(199, "Number of entries: " << hm->GetEntries());
54 TH1D* hlast = (TH1D*)hm->Clone("hlast");
55 hlast->Reset();
56 hlast->SetName(m_name);
57 for (int i = 1; i < h2->GetNbinsX(); ++i) {
58 double sum = 0;
59 double err = 0;
60 double mean = -99;
61 TH1D* h1d = h2->ProjectionY("h1d", i, i);
62 if (!h1d) continue;
63 sum = h1d->GetEntries();
64 if (sum < minHitCut) continue; //skip low data slice
65 mean = h1d->GetMean();
66 double sg = h1d->GetRMS();
67 double max = h1d->GetMaximum();
68 g1->SetParameters(max, mean, sg);
69 h1d->Fit("g1", "QNR", "");
70 // TF1 *f1=h1d->GetFunction("gaus");
71 mean = g1->GetParameter(1);
72 err = g1->GetParError(1);
73 if (sum > 50) {
74 double sg2 = g1->GetParameter(2);
75 h1d->Fit("g1", "Q0", "", mean - 1.1 * sg2, mean + 1.1 * sg2);
76 mean = g1->GetParameter(1);
77 // err=g1->GetParError(1);
78 }
79 hlast->SetBinContent(i, mean);
80 hlast->SetBinError(i, err);
81 h1d->Delete();
82 }
83 return hlast;
84 }
85 };
87}
Class to do the slice fit.
Definition: SliceFit.h:25
static TH1D * doSliceFitY(TH2D *h2, int minHitCut=0)
Do the slice fit for the 2d histogram.
Definition: SliceFit.h:30
Abstract base class for different kinds of events.