CbmRoot
PairAnalysisSignalFunc.cxx
Go to the documentation of this file.
1 // Dielectron Function //
3 // //
4 // //
5 /*
6 
7  1) add any background function you like
8  TF1 *fB = new TF1("fitBgrd","pol3",minFit,maxFit);
9 
10  2) configure the signal extraction
11 
12  2.1) chose one of the signal functions (MCshape, CrystalBall, Gauss)
13  TF1 *fS = new TF1("fitSign",PairAnalysisFunction::PeakFunCB,minFit,maxFit,5); // has 5 parameters
14  // TF1 *fS = new TF1("fitSign",PairAnalysisFunction::PeakFunGaus,minFit,maxFit,3); // has 3 parameters
15  // sig->SetMCSignalShape(hMCsign);
16  // TF1 *fS = new TF1("fitSign",PairAnalysisFunction::PeakFunMC,minFit,maxFit,1); // requires a MC shape
17 
18 
19  3) combined fit of bgrd+signal
20 
21  3.1) combine the two functions
22  sig->CombineFunc(fS,fB);
23  3.2) apply fitting ranges and the fit options
24  sig->SetFitRange(minFit, maxFit);
25  sig->SetFitOption("NR");
26 
27 
28  6) access the spectra and values created
29 
30  6.1) fit function
31  TF1 *fFitSign = sig->GetCombinedFunction(); // combined fit function
32  TF1 *fFitExtr = sig->GetSignalFunction(); // signal function
33  TF1 *fFitBgrd = sig->GetBackgroundFunction(); // background function
34 
35 */
36 // //
38 
39 #include <TF1.h>
40 #include <TFitResult.h>
41 #include <TGraph.h>
42 #include <TH1.h>
43 #include <TList.h>
44 #include <TMath.h>
45 #include <TPaveText.h>
46 #include <TString.h>
47 //#include <../hist/hist/src/TF1Helper.h>
48 
49 #include "PairAnalysisFunction.h"
50 
52 
53 
55  : TNamed("function", "function") {
56  //
57  // Default Constructor
58  //
59 }
60 
61 //______________________________________________
62 PairAnalysisFunction::PairAnalysisFunction(const char* name, const char* title)
63  : TNamed(name, title) {
64  //
65  // Named Constructor
66  //
67 }
68 
69 //______________________________________________
71  //
72  // Default Destructor
73  //
74  if (fFuncSigBack) delete fFuncSigBack;
75  if (fFuncSignal) delete fFuncSignal;
76  if (fFuncBackground) delete fFuncBackground;
77 }
78 
79 
80 //______________________________________________________________________________
81 Double_t PairAnalysisFunction::PeakFunMC(const Double_t* x,
82  const Double_t* par) {
83  // Fit MC signal shape
84  // parameters
85  // [0]: scale for simpeak
86 
87  Double_t xx = x[0];
88 
89  TH1F* hPeak = fgHistSimPM;
90  if (!hPeak) {
91  printf(
92  "E-PairAnalysisFunction::PeakFun: No histogram for peak fit defined!\n");
93  return 0.0;
94  }
95 
96  Int_t idx = hPeak->FindBin(xx);
97  if ((idx <= 1) || (idx >= hPeak->GetNbinsX())) { return 0.0; }
98 
99  return (par[0] * hPeak->GetBinContent(idx));
100 }
101 
102 //______________________________________________________________________________
103 Double_t PairAnalysisFunction::PeakFunCB(const Double_t* x,
104  const Double_t* par) {
105  // Crystal Ball fit function
106 
107  Double_t n = par[0];
108  Double_t alpha = par[1];
109  Double_t meanx = par[2];
110  Double_t sigma = par[3];
111  Double_t nn = par[4];
112 
113  Double_t a =
114  TMath::Power((n / TMath::Abs(alpha)), n) * TMath::Exp(-.5 * alpha * alpha);
115  Double_t b = n / TMath::Abs(alpha) - TMath::Abs(alpha);
116 
117  Double_t arg = (x[0] - meanx) / sigma;
118  Double_t fitval = 0;
119 
120  if (arg > -1. * alpha) {
121  fitval = nn * TMath::Exp(-.5 * arg * arg);
122  } else {
123  fitval = nn * a * TMath::Power((b - arg), (-1 * n));
124  }
125 
126  return fitval;
127 }
128 
129 //______________________________________________________________________________
130 Double_t PairAnalysisFunction::PeakFunGaus(const Double_t* x,
131  const Double_t* par) {
132  // Gaussian fit function
133  //printf("fNparBgrd %d \n",fNparBgnd);
134  Double_t n = par[0];
135  Double_t mean = par[1];
136  Double_t sigma = par[2];
137  Double_t xx = x[0];
138 
139  return (n * TMath::Exp(-0.5 * TMath::Power((xx - mean) / sigma, 2)));
140 }
141 
142 //______________________________________________
143 void PairAnalysisFunction::SetFunctions(TF1* const combined,
144  TF1* const sig,
145  TF1* const back,
146  Int_t parM,
147  Int_t parMres) {
148  //
149  // Set the signal, background functions and combined fit function
150  // Note: The process method assumes that the first n parameters in the
151  // combined fit function correspond to the n parameters of the signal function
152  // and the n+1 to n+m parameters to the m parameters of the background function!!!
153 
154  if (!sig || !back || !combined) {
155  Error("SetFunctions",
156  "Both, signal and background function need to be set!");
157  return;
158  }
159  fFuncSignal = sig;
160  fFuncBackground = back;
161  fFuncSigBack = combined;
162  fParMass = parM;
163  fParMassWidth = parMres;
164 }
165 
166 //______________________________________________
167 void PairAnalysisFunction::SetDefaults(Int_t type) {
168  //
169  // Setup some default functions:
170  // type = 0: gaus signal + linear background in 2.5 - 4 GeV inv. mass
171  // type = 1: gaus signal + exponential background in 2.5 - 4 GeV inv. mass
172  // type = 2: half gaussian, half exponential signal function
173  // type = 3: Crystal-Ball function
174  // type = 4: Crystal-Ball signal + exponential background
175  //
176  // TODO: use PDG mass and width + fPOI for parameter settings
177 
178  if (type == 0) {
179  fFuncSignal = new TF1("DieleSignal", "gaus", 2.5, 4);
180  fFuncBackground = new TF1("DieleBackground", "pol1", 2.5, 4);
181  fFuncSigBack = new TF1("DieleCombined", "gaus+pol1(3)", 2.5, 4);
182 
183  fFuncSigBack->SetParameters(1, 3.1, .05, 2.5, 1);
184  fFuncSigBack->SetParLimits(0, 0, 10000000);
185  fFuncSigBack->SetParLimits(1, 3.05, 3.15);
186  fFuncSigBack->SetParLimits(2, .02, .1);
187  } else if (type == 1) {
188  fFuncSignal = new TF1("DieleSignal", "gaus", 2.5, 4);
190  new TF1("DieleBackground", "[0]*exp(-(x-[1])/[2])", 2.5, 4);
191  fFuncSigBack =
192  new TF1("DieleCombined", "gaus+[3]*exp(-(x-[4])/[5])", 2.5, 4);
193 
194  fFuncSigBack->SetParameters(1, 3.1, .05, 1, 2.5, 1);
195  fFuncSigBack->SetParLimits(0, 0, 10000000);
196  fFuncSigBack->SetParLimits(1, 3.05, 3.15);
197  fFuncSigBack->SetParLimits(2, .02, .1);
198  } else if (type == 2) {
199  // half gaussian, half exponential signal function
200  // exponential background
201  fFuncSignal = new TF1("DieleSignal",
202  "(x<[1])*([0]*(exp(-0.5*((x-[1])/[2])^2)+exp((x-[1])/"
203  "[3])*(1-exp(-0.5*((x-[1])/"
204  "[2])^2))))+(x>=[1])*([0]*exp(-0.5*((x-[1])/[2])^2))",
205  2.5,
206  4);
208  new TF1("DieleBackground", "[0]*exp(-(x-[1])/[2])+[3]", 2.5, 4);
209  fFuncSigBack = new TF1(
210  "DieleCombined",
211  "(x<[1])*([0]*(exp(-0.5*((x-[1])/[2])^2)+exp((x-[1])/"
212  "[3])*(1-exp(-0.5*((x-[1])/[2])^2))))+(x>=[1])*([0]*exp(-0.5*((x-[1])/"
213  "[2])^2))+[4]*exp(-(x-[5])/[6])+[7]",
214  2.5,
215  4);
216  fFuncSigBack->SetParameters(1., 3.1, .05, .1, 1, 2.5, 1, 0);
217 
218  fFuncSigBack->SetParLimits(0, 0, 10000000);
219  fFuncSigBack->SetParLimits(1, 3.05, 3.15);
220  fFuncSigBack->SetParLimits(2, .02, .1);
221  fFuncSigBack->FixParameter(6, 2.5);
222  fFuncSigBack->FixParameter(7, 0);
223  }
224 }
225 
226 //______________________________________________________________________________
227 void PairAnalysisFunction::CombineFunc(TF1* const peak, TF1* const bgnd) {
228  //
229  // combine the bgnd and the peak function
230  //
231 
232  if (!peak || !bgnd) {
233  Error("CombineFunc",
234  "Both, signal and background function need to be set!");
235  return;
236  }
237  fFuncSignal = peak;
238  fFuncBackground = bgnd;
239 
240  fNparPeak = fFuncSignal->GetNpar();
241  fNparBgnd = fFuncBackground->GetNpar();
242 
243  fFuncSigBack = new TF1("BgndPeak",
244  this,
246  fFitMin,
247  fFitMax,
248  fNparPeak + fNparBgnd);
249  return;
250 }
251 
252 //______________________________________________________________________________
253 Double_t PairAnalysisFunction::PeakBgndFun(const Double_t* x,
254  const Double_t* par) {
255  //
256  // merge peak and bgnd functions
257  //
258  return (
259  fFuncSignal->EvalPar(x, par)
260  + (fFuncBackground ? fFuncBackground->EvalPar(x, par + fNparPeak) : 0.));
261 }
262 
263 //______________________________________________
264 void PairAnalysisFunction::Print(Option_t* /*option*/) const {
265  //
266  // Print the statistics
267  //
268  printf("Fit int. : %.5g - %.5g \n", fFitMin, fFitMax);
269  printf("Fit chi/%d: %.5g \n", fDof, fChi2Dof);
270 }
PairAnalysisFunction
Definition: PairAnalysisFunction.h:22
PairAnalysisFunction::fDof
Int_t fDof
Definition: PairAnalysisFunction.h:120
PairAnalysisFunction::CombineFunc
void CombineFunc(TF1 *const peak=0, TF1 *const bgnd=0)
Definition: PairAnalysisFunction.cxx:314
PairAnalysisFunction::SetFunctions
void SetFunctions(TF1 *const combined, TF1 *const sig=0, TF1 *const back=0, Int_t parM=1, Int_t parMres=2)
Definition: PairAnalysisFunction.cxx:216
PairAnalysisFunction::fFitMax
Double_t fFitMax
Definition: PairAnalysisFunction.h:107
PairAnalysisFunction::Print
virtual void Print(Option_t *option="") const
Definition: PairAnalysisSignalFunc.cxx:264
PairAnalysisFunction::fChi2Dof
Double_t fChi2Dof
Definition: PairAnalysisFunction.h:121
PairAnalysisFunction::SetDefaults
void SetDefaults(Int_t type)
Definition: PairAnalysisFunction.cxx:254
PairAnalysisFunction::fParMass
Int_t fParMass
Definition: PairAnalysisFunction.h:111
PairAnalysisFunction::fFitMin
Double_t fFitMin
Definition: PairAnalysisFunction.h:106
PairAnalysisFunction::fFuncBackground
TF1 * fFuncBackground
Definition: PairAnalysisFunction.h:103
PairAnalysisFunction::fFuncSigBack
TF1 * fFuncSigBack
Definition: PairAnalysisFunction.h:104
PairAnalysisFunction::PairAnalysisFunction
PairAnalysisFunction()
PairAnalysisFunction.h
PairAnalysisFunction::fFuncSignal
TF1 * fFuncSignal
Definition: PairAnalysisFunction.h:102
PairAnalysisFunction::PeakFunGaus
Double_t PeakFunGaus(const Double_t *x, const Double_t *par)
Definition: PairAnalysisSignalFunc.cxx:130
ClassImp
ClassImp(PairAnalysisFunction) PairAnalysisFunction
Definition: PairAnalysisSignalFunc.cxx:51
PairAnalysisFunction::fNparPeak
Int_t fNparPeak
Definition: PairAnalysisFunction.h:123
PairAnalysisFunction::PeakBgndFun
Double_t PeakBgndFun(const Double_t *x, const Double_t *par)
Definition: PairAnalysisFunction.cxx:339
PairAnalysisFunction::PeakFunMC
Double_t PeakFunMC(const Double_t *x, const Double_t *par)
Definition: PairAnalysisFunction.cxx:99
PairAnalysisFunction::fgHistSimPM
static TH1F * fgHistSimPM
Definition: PairAnalysisFunction.h:96
x
Double_t x
Definition: CbmMvdSensorDigiToHitTask.cxx:68
PairAnalysisFunction::fParMassWidth
Int_t fParMassWidth
Definition: PairAnalysisFunction.h:113
PairAnalysisFunction::~PairAnalysisFunction
virtual ~PairAnalysisFunction()
Definition: PairAnalysisFunction.cxx:88
PairAnalysisFunction::fNparBgnd
Int_t fNparBgnd
Definition: PairAnalysisFunction.h:124
PairAnalysisFunction::PeakFunCB
Double_t PeakFunCB(const Double_t *x, const Double_t *par)
Definition: PairAnalysisFunction.cxx:121