CbmRoot
PairAnalysisFunction.cxx
Go to the documentation of this file.
1 // PairAnalysis 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, PowGaussPow, ExpGaussExp)
13  TF1 *fS = new TF1("fitSign",PairAnalysisFunction::PeakFunCB,minFit,maxFit,5); // has 5 parameters
14  // sig->SetMCSignalShape(hMCsign);
15  // TF1 *fS = new TF1("fitSign",PairAnalysisFunction::PeakFunMC,minFit,maxFit,1); // requires a MC shape
16 
17  OR
18 
19  2.2) one of the other predefined function (Boltzmann, PtExp, Hagedorn, Levi)
20 
21 
22  3) combined fit of bgrd+signal
23 
24  3.1) combine the two functions
25  sig->CombineFunc(fS,fB);
26  3.2) apply fitting ranges and the fit options
27  sig->SetFitRange(minFit, maxFit);
28  sig->SetFitOption("NR");
29 
30 
31  6) access the spectra and values created
32 
33  6.1) fit function
34  TF1 *fFitSign = sig->GetCombinedFunction(); // combined fit function
35  TF1 *fFitExtr = sig->GetSignalFunction(); // signal function
36  TF1 *fFitBgrd = sig->GetBackgroundFunction(); // background function
37 
38 */
39 // //
41 
42 #include <TF1.h>
43 #include <TFitResult.h>
44 #include <TGraph.h>
45 #include <TH1.h>
46 #include <TList.h>
47 #include <TMath.h>
48 #include <TNamed.h>
49 #include <TPaveText.h>
50 #include <TString.h>
51 //#include <../hist/hist/src/TF1Helper.h>
52 
53 #include "PairAnalysisFunction.h"
54 
55 
57 
59 
61  : PairAnalysisFunction("PairAnalysisFunction", "fcttitle") {
62  //
63  // Default Constructor
64  //
65 }
66 
67 //______________________________________________
68 PairAnalysisFunction::PairAnalysisFunction(const char* name, const char* title)
69  : TNamed(name, title) {
70  //
71  // Named Constructor
72  //
73 }
74 
75 //______________________________________________
77  : TNamed(c.GetName(), c.GetTitle())
78  , fFitMin(c.GetFitMin())
79  , fFitMax(c.GetFitMax())
80  , fPOIpdg(c.GetParticleOfInterest()) {
81  //
82  // Copy Constructor
83  //
84 }
85 
86 
87 //______________________________________________
89  //
90  // Default Destructor
91  //
92  if (fFuncSigBack) delete fFuncSigBack;
93  if (fFuncSignal) delete fFuncSignal;
94  if (fFuncBackground) delete fFuncBackground;
95 }
96 
97 
98 //______________________________________________________________________________
99 Double_t PairAnalysisFunction::PeakFunMC(const Double_t* x,
100  const Double_t* par) {
101  // Fit MC signal shape
102  // parameters
103  // [0]: scale for simpeak
104 
105  Double_t xx = x[0];
106 
107  TH1F* hPeak = fgHistSimPM;
108  if (!hPeak) {
109  printf(
110  "E-PairAnalysisFunction::PeakFun: No histogram for peak fit defined!\n");
111  return 0.0;
112  }
113 
114  Int_t idx = hPeak->FindBin(xx);
115  if ((idx <= 1) || (idx >= hPeak->GetNbinsX())) { return 0.0; }
116 
117  return (par[0] * hPeak->GetBinContent(idx));
118 }
119 
120 //______________________________________________________________________________
121 Double_t PairAnalysisFunction::PeakFunCB(const Double_t* x,
122  const Double_t* par) {
123  // Crystal Ball fit function
124 
125  Double_t n = par[0];
126  Double_t alpha = par[1];
127  Double_t meanx = par[2];
128  Double_t sigma = par[3];
129  Double_t nn = par[4];
130 
131  Double_t a =
132  TMath::Power((n / TMath::Abs(alpha)), n) * TMath::Exp(-.5 * alpha * alpha);
133  Double_t b = n / TMath::Abs(alpha) - TMath::Abs(alpha);
134 
135  Double_t arg = (x[0] - meanx) / sigma;
136  Double_t fitval = 0;
137 
138  if (arg > -1. * alpha) {
139  fitval = nn * TMath::Exp(-.5 * arg * arg); //gaussian part
140  } else {
141  fitval = nn * a * TMath::Power((b - arg), (-1 * n));
142  }
143 
144  return fitval;
145 }
146 
147 //______________________________________________________________________________
149  const Double_t* par) {
150  // PowGaussPow function fit function (both sided Crystall Ball)
151 
152  Double_t n = par[0];
153  Double_t alpha = par[1];
154  Double_t nn = par[4];
155  Double_t meanx = par[2];
156  Double_t sigma = par[3];
157 
158  Double_t a =
159  TMath::Power((n / TMath::Abs(alpha)), n) * TMath::Exp(-.5 * alpha * alpha);
160  Double_t b = n / TMath::Abs(alpha) - TMath::Abs(alpha);
161 
162  Double_t arg = (x[0] - meanx) / sigma;
163  Double_t fitval = 0;
164 
165  if (arg > alpha) {
166  fitval = nn * a * TMath::Power((b + arg), (-1 * n));
167  } else if (arg < -alpha) {
168  fitval = nn * a * TMath::Power((b - arg), (-1 * n));
169  } else {
170  fitval = nn * TMath::Exp(-0.5 * arg * arg); //gaussian part
171  }
172 
173  return fitval;
174 }
175 
176 //______________________________________________________________________________
178  const Double_t* par) {
179  // ExpGaussExp function fit function
180 
181  Double_t n = par[0];
182  Double_t alpha = par[1];
183 
184  Double_t meanx = par[2];
185  Double_t sigma = par[3];
186 
187  Double_t arg = (x[0] - meanx) / sigma;
188  Double_t fitval = 0;
189 
190  if (arg > alpha) {
191  fitval = n * TMath::Exp(-0.5 * alpha * alpha - alpha * arg);
192  } else if (arg < -alpha) {
193  fitval = n * TMath::Exp(-0.5 * alpha * alpha + alpha * arg);
194  } else {
195  fitval = n * TMath::Exp(-0.5 * arg * arg); //gaussian part
196  }
197 
198  return fitval;
199 }
200 
201 //______________________________________________________________________________
202 Double_t PairAnalysisFunction::PeakFunGauss(const Double_t* x,
203  const Double_t* par) {
204  // Gaussian fit function
205 
206  //printf("fNparBgrd %d \n",fNparBgnd);
207  Double_t n = par[0];
208  Double_t mean = par[1];
209  Double_t sigma = par[2];
210  Double_t xx = x[0];
211 
212  return (n * TMath::Exp(-0.5 * TMath::Power((xx - mean) / sigma, 2)));
213 }
214 
215 //______________________________________________
216 void PairAnalysisFunction::SetFunctions(TF1* const combined,
217  TF1* const sig,
218  TF1* const back,
219  Int_t parM,
220  Int_t parMres) {
221  //
222  // Set the signal, background functions and combined fit function
223  // Note: The process method assumes that the first n parameters in the
224  // combined fit function correspond to the n parameters of the signal function
225  // and the n+1 to n+m parameters to the m parameters of the background function!!!
226 
227  if (!sig || !back || !combined) {
228  Error("SetFunctions",
229  "Both, signal and background function need to be set!");
230  return;
231  }
232  fFuncSignal = sig;
233  fFuncBackground = back;
234  fFuncSigBack = combined;
235  fParMass = parM;
236  fParMassWidth = parMres;
237 }
238 
239 //______________________________________________
244  switch (predefinedFunc) {
245  case kBoltzmann: GetBoltzmann(); break;
246  case kPtExp: GetPtExp(); break;
247  case kHagedorn: GetHagedorn(); break;
248  case kLevi: GetLevi(); break;
249  default: Error("SetDefault", "predefined function not yet implemented");
250  }
251 }
252 
253 //______________________________________________
255  //
256  // Setup some default functions:
257  // type = 0: gaus signal + linear background in 2.5 - 4 GeV inv. mass
258  // type = 1: gaus signal + exponential background in 2.5 - 4 GeV inv. mass
259  // type = 2: half gaussian, half exponential signal function
260  // type = 3: Crystal-Ball function
261  // type = 4: Crystal-Ball signal + exponential background
262  //
263  // TODO: use PDG mass and width + fPOI for parameter settings
264 
265  if (type == 0) {
266  fFuncSignal = new TF1("DieleSignal", "gaus", 2.5, 4);
267  fFuncBackground = new TF1("DieleBackground", "pol1", 2.5, 4);
268  fFuncSigBack = new TF1("DieleCombined", "gaus+pol1(3)", 2.5, 4);
269 
270  fFuncSigBack->SetParameters(1, 3.1, .05, 2.5, 1);
271  fFuncSigBack->SetParLimits(0, 0, 10000000);
272  fFuncSigBack->SetParLimits(1, 3.05, 3.15);
273  fFuncSigBack->SetParLimits(2, .02, .1);
274  } else if (type == 1) {
275  fFuncSignal = new TF1("DieleSignal", "gaus", 2.5, 4);
277  new TF1("DieleBackground", "[0]*exp(-(x-[1])/[2])", 2.5, 4);
278  fFuncSigBack =
279  new TF1("DieleCombined", "gaus+[3]*exp(-(x-[4])/[5])", 2.5, 4);
280 
281  fFuncSigBack->SetParameters(1, 3.1, .05, 1, 2.5, 1);
282  fFuncSigBack->SetParLimits(0, 0, 10000000);
283  fFuncSigBack->SetParLimits(1, 3.05, 3.15);
284  fFuncSigBack->SetParLimits(2, .02, .1);
285  } else if (type == 2) {
286  // half gaussian, half exponential signal function
287  // exponential background
288  fFuncSignal = new TF1("DieleSignal",
289  "(x<[1])*([0]*(exp(-0.5*((x-[1])/[2])^2)+exp((x-[1])/"
290  "[3])*(1-exp(-0.5*((x-[1])/"
291  "[2])^2))))+(x>=[1])*([0]*exp(-0.5*((x-[1])/[2])^2))",
292  2.5,
293  4);
295  new TF1("DieleBackground", "[0]*exp(-(x-[1])/[2])+[3]", 2.5, 4);
296  fFuncSigBack = new TF1(
297  "DieleCombined",
298  "(x<[1])*([0]*(exp(-0.5*((x-[1])/[2])^2)+exp((x-[1])/"
299  "[3])*(1-exp(-0.5*((x-[1])/[2])^2))))+(x>=[1])*([0]*exp(-0.5*((x-[1])/"
300  "[2])^2))+[4]*exp(-(x-[5])/[6])+[7]",
301  2.5,
302  4);
303  fFuncSigBack->SetParameters(1., 3.1, .05, .1, 1, 2.5, 1, 0);
304 
305  fFuncSigBack->SetParLimits(0, 0, 10000000);
306  fFuncSigBack->SetParLimits(1, 3.05, 3.15);
307  fFuncSigBack->SetParLimits(2, .02, .1);
308  fFuncSigBack->FixParameter(6, 2.5);
309  fFuncSigBack->FixParameter(7, 0);
310  }
311 }
312 
313 //______________________________________________________________________________
314 void PairAnalysisFunction::CombineFunc(TF1* const peak, TF1* const bgnd) {
315  //
316  // combine the bgnd and the peak function
317  //
318 
319  if (!peak) {
320  Error("CombineFunc", "signal function need to be set!");
321  return;
322  }
323  fFuncSignal = peak;
324  fFuncBackground = bgnd;
325 
326  fNparPeak = fFuncSignal->GetNpar();
327  fNparBgnd = (bgnd ? fFuncBackground->GetNpar() : 0);
328 
329  fFuncSigBack = new TF1("BgndPeak",
330  this,
332  fFitMin,
333  fFitMax,
334  fNparPeak + fNparBgnd);
335  return;
336 }
337 
338 //______________________________________________________________________________
339 Double_t PairAnalysisFunction::PeakBgndFun(const Double_t* x,
340  const Double_t* par) {
341  //
342  // merge peak and bgnd functions
343  //
344  return (
345  fFuncSignal->EvalPar(x, par)
346  + (fFuncBackground ? fFuncBackground->EvalPar(x, par + fNparPeak) : 0.));
347 }
348 
349 //______________________________________________
350 //void PairAnalysisFunction::Print(Option_t */*option*/) const
351 //{
352 //
353 // Print the statistics
354 //
355 // printf("Fit int. : %.5g - %.5g \n",fFitMin,fFitMax);
356 // printf("Fit chi/%d: %.5g \n",fDof,fChi2Dof);
357 //}
358 
359 
360 //______________________________________________
362  // Boltzmann (exp in 1/mt*dNdmT times mt) as a function of dNdpt
363  fFuncSigBack = new TF1("Boltzmann",
364  "[1]*x*sqrt(x*x+[0]*[0])*exp(-sqrt(x*x+[0]*[0])/[2])",
365  0.,
366  10.);
367  // fFuncSigBack->SetParameters(fPOI->Mass(), norm, temp);
368  if (fPOI) fFuncSigBack->FixParameter(0, fPOI->Mass());
369  fFuncSigBack->SetParLimits(2, 0.01, 10);
370  fFuncSigBack->SetParNames("mass", "norm", "T");
371  return fFuncSigBack;
372 }
373 
374 //______________________________________________
376  // Simple exponential in 1/pt*dNdpT, as a function of dNdpt
377  fFuncSigBack = new TF1("Exponential", "[0]*x*exp(-x/[1])", 0., 10.);
378  // fFuncSigBack->SetParameters(norm, temp);
379  fFuncSigBack->SetParLimits(1, 0.01, 10);
380  fFuncSigBack->SetParNames("norm", "T");
381  return fFuncSigBack;
382 }
383 
384 //______________________________________________
386  // PowerLaw function, dNdpt
387  // power law Nuclear Physics B, Vol. 335, No. 2. (7 May 1990), pp. 261-287.
388  // This is sometimes also called Hagedorn or modified Hagedorn
389  fFuncSigBack = new TF1("Hagedorn", "x*[0]*( 1 + x/[1] )^(-[2])", 0., 10.);
390  // fFuncSigBack->SetParameters(norm, pt0, n);
391  fFuncSigBack->SetParLimits(1, 0.01, 10);
392  //fFuncSigBack->SetParLimits(2, 0.01, 50);
393  fFuncSigBack->SetParNames("norm", "pt0", "n");
394  return fFuncSigBack;
395 }
396 
397 //______________________________________________
399  // Levi function (aka Tsallis), dNdpt
400  fFuncSigBack =
401  new TF1("Levi-Tsallis",
402  "( x*[0]*([1]-1)*([1]-2) )/( [1]*[2]*( [1]*[2]+[3]*([1]-2) ) ) * "
403  "( 1 + (sqrt([3]*[3]+x*x) -[3])/([1]*[2]) )^(-[1])",
404  0.,
405  10.);
406  // fFuncSigBack->SetParameters(norm, n, temp,mass);
407  if (fPOI) fFuncSigBack->FixParameter(3, fPOI->Mass());
408  fFuncSigBack->SetParLimits(2, 0.01, 10);
409  fFuncSigBack->SetParNames("norm (dN/dy)", "n", "T", "mass");
410  return fFuncSigBack;
411 }
PairAnalysisFunction
Definition: PairAnalysisFunction.h:22
PairAnalysisFunction::EFunction
EFunction
Definition: PairAnalysisFunction.h:25
PairAnalysisFunction::GetPtExp
TF1 * GetPtExp()
Definition: PairAnalysisFunction.cxx:375
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::GetBoltzmann
TF1 * GetBoltzmann()
Definition: PairAnalysisFunction.cxx:361
PairAnalysisFunction::kBoltzmann
@ kBoltzmann
Definition: PairAnalysisFunction.h:25
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::SetDefault
void SetDefault(EFunction predefinedFunc)
Definition: PairAnalysisFunction.cxx:240
PairAnalysisFunction::PeakFunGauss
Double_t PeakFunGauss(const Double_t *x, const Double_t *par)
Definition: PairAnalysisFunction.cxx:202
PairAnalysisFunction::fFuncSignal
TF1 * fFuncSignal
Definition: PairAnalysisFunction.h:102
ClassImp
ClassImp(CbmConverterManager) InitStatus CbmConverterManager
Definition: CbmConverterManager.cxx:12
PairAnalysisFunction::kPtExp
@ kPtExp
Definition: PairAnalysisFunction.h:25
PairAnalysisFunction::fPOI
TParticlePDG * fPOI
Definition: PairAnalysisFunction.h:109
PairAnalysisFunction::fNparPeak
Int_t fNparPeak
Definition: PairAnalysisFunction.h:123
PairAnalysisFunction::kHagedorn
@ kHagedorn
Definition: PairAnalysisFunction.h:25
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::kLevi
@ kLevi
Definition: PairAnalysisFunction.h:25
PairAnalysisFunction::GetHagedorn
TF1 * GetHagedorn()
Definition: PairAnalysisFunction.cxx:385
PairAnalysisFunction::fParMassWidth
Int_t fParMassWidth
Definition: PairAnalysisFunction.h:113
PairAnalysisFunction::~PairAnalysisFunction
virtual ~PairAnalysisFunction()
Definition: PairAnalysisFunction.cxx:88
PairAnalysisFunction::PeakFunPowGaussPow
Double_t PeakFunPowGaussPow(const Double_t *x, const Double_t *par)
Definition: PairAnalysisFunction.cxx:148
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
PairAnalysisFunction::PeakFunExpGaussExp
Double_t PeakFunExpGaussExp(const Double_t *x, const Double_t *par)
Definition: PairAnalysisFunction.cxx:177
PairAnalysisFunction::GetLevi
TF1 * GetLevi()
Definition: PairAnalysisFunction.cxx:398