Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
CBfitter_jpsi_LL_single.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file CBfitter_jpsi_LL_single.C
1 #include <TF1.h>
2 #include <TMath.h>
3 #include <TFile.h>
4 #include <TCanvas.h>
5 #include <TH1.h>
6 #include <TH2.h>
7 #include <TGraph.h>
8 #include <TLegend.h>
9 #include <TStyle.h>
10 #include <TROOT.h>
11 #include <TLatex.h>
12 
13 TH1D *recomass;
14 TH1D *recomassFG;
15 TH1D *recomassBG;
16 
17 Double_t CBcalc(Double_t *xx, Double_t *par)
18 {
19  // Crystal Ball lineshape plus exponential background
20  // Used for chisq fit to mixed event background data
21 
22  double f;
23  double x = xx[0];
24 
25  // The four CB parameters (alpha, n, x_mean, sigma) plus normalization (N) are:
26  double alpha = par[0];
27  double n = par[1];
28  double x_mean = par[2];
29  double sigma = par[3];
30  double N = par[4];
31 
32  // we need:
33  double A = pow( (n/TMath::Abs(alpha)),n) * exp(-pow(alpha,2)/2.0);
34  double B = n/TMath::Abs(alpha) - TMath::Abs(alpha);
35 
36  // The Crystal Ball function is:
37  if( (x-x_mean)/sigma > -alpha)
38  {
39  f = N * exp( -pow(x-x_mean,2) / (2.0*pow(sigma,2)));
40  }
41  else
42  {
43  f = N * A * pow(B - (x-x_mean)/sigma, -n);
44  }
45 
46  return f;
47 }
48 
49 Double_t CBcalc_exp(Double_t *xx, Double_t *par)
50 {
51  // Crystal Ball lineshape plus exponential background
52  // Used for chisq fit to mixed event background data
53 
54  double f;
55  double x = xx[0];
56 
57  // The four CB parameters (alpha, n, x_mean, sigma) plus normalization (N) are:
58  double alpha = par[0];
59  double n = par[1];
60  double x_mean = par[2];
61  double sigma = par[3];
62  double N = par[4];
63  // and we add an exponential background
64  double Nexp = par[5];
65  double slope = par[6];
66 
67  // we need:
68  double A = pow( (n/TMath::Abs(alpha)),n) * exp(-pow(alpha,2)/2.0);
69  double B = n/TMath::Abs(alpha) - TMath::Abs(alpha);
70 
71  // The Crystal Ball function is:
72  if( (x-x_mean)/sigma > -alpha)
73  {
74  f = N * exp( -pow(x-x_mean,2) / (2.0*pow(sigma,2)));
75  }
76  else
77  {
78  f = N * A * pow(B - (x-x_mean)/sigma, -n);
79  }
80 
81  f += Nexp * exp(slope * x);
82 
83  return f;
84 }
85 
86 Double_t CBcalc_LL(Double_t *xx, Double_t *par)
87 {
88  // Crystal Ball lineshape plus exponential background plus mixed event combinatorial background
89  // Used for LL fit to foreground data
90  double f;
91  double x = xx[0];
92 
93  // The four CB parameters (alpha, n, x_mean, sigma) plus normalization (N) are:
94  double alpha = par[0];
95  double n = par[1];
96  double x_mean = par[2];
97  double sigma = par[3];
98  double N = par[4];
99  // and we add an exponential background
100  double Nexp = par[5];
101  double slope = par[6];
102 
103  /*
104  TF1 *fcbexp = new TF1("fcbexp",CBcalc_exp,2.0,4.0,7);
105  fcbexp->SetParameter(0,alpha);
106  fcbexp->SetParameter(1,n);
107  fcbexp->SetParameter(2,x_mean);
108  fcbexp->SetParameter(3,sigma);
109  fcbexp->SetParameter(4, N);
110  fcbexp->SetParameter(5,Nexp);
111  fcbexp->SetParameter(6,slope);
112  f = fcbexp->Eval(x);
113  */
114 
115  // we need:
116  double A = pow( (n/TMath::Abs(alpha)),n) * exp(-pow(alpha,2)/2.0);
117  double B = n/TMath::Abs(alpha) - TMath::Abs(alpha);
118 
119  // The Crystal Ball function is:
120  if( (x-x_mean)/sigma > -alpha)
121  {
122  f = N * exp( -pow(x-x_mean,2) / (2.0*pow(sigma,2)));
123  }
124  else
125  {
126  f = N * A * pow(B - (x-x_mean)/sigma, -n);
127  }
128 
129  f += Nexp * exp(slope * x);
130 
131  // Add the mixed event background to the fit function
132  int ibin = recomassBG->FindBin(x);
133 
134  double bg12 = recomassBG->GetBinContent(ibin);
135  f += bg12;
136 
137 
138  return f;
139 }
140 
141 Double_t Upscalc(Double_t *xx, Double_t *par)
142 {
143  double x = xx[0];
144 
145  // The input parameters are:
146  double alpha1s = par[0];
147  double n1s = par[1];
148  double sigma1s = par[2];
149  double m1s = par[3];
150  double m2s = par[4];
151  double m3s = par[5];
152  double N1s = par[6];
153  double N2s = par[7];
154  double N3s = par[8];
155  double Nexp = par[9];
156  double slope = par[10];
157 
158  // Construct the Y(1S) CB shape
159 
160  TF1 *f1 = new TF1("f1",CBcalc,7,11,5);
161  f1->SetParameter(0,alpha1s);
162  f1->SetParameter(1,n1s);
163  f1->SetParameter(2,m1s);
164  f1->SetParameter(3,sigma1s);
165  f1->SetParameter(4,N1s);
166 
167  TF1 *f2 = new TF1("f2",CBcalc,7,11,5);
168  f2->SetParameter(0,alpha1s);
169  f2->SetParameter(1,n1s);
170  f2->SetParameter(2,m2s);
171  f2->SetParameter(3,sigma1s);
172  f2->SetParameter(4,N2s);
173 
174  TF1 *f3 = new TF1("f3",CBcalc,7,11,5);
175  f3->SetParameter(0,alpha1s);
176  f3->SetParameter(1,n1s);
177  f3->SetParameter(2,m3s);
178  f3->SetParameter(3,sigma1s);
179  f3->SetParameter(4,N3s);
180 
181  TF1 *fexp = new TF1("fexp","[0]*exp([1]*x)",7,11);
182  fexp->SetParameter(0,Nexp);
183  fexp->SetParameter(1,slope);
184 
185 
186  double mass = f1->Eval(x) + f2->Eval(x) + f3->Eval(x) + fexp->Eval(x);
187 
188  return mass;
189 }
190 
191 void CBfitter_jpsi_LL_single()
192 {
193  gROOT->SetStyle("Plain");
194  gStyle->SetOptStat(0);
195  gStyle->SetOptFit(0);
196  gStyle->SetOptTitle(0);
197 
198  bool LL_fit = true;
199 
200  double maxmass = 3.5;
201 
202  TFile *file1S;
203  file1S = new TFile("Final_Hists_mee_meeBG_Type_VetoB0.root");
204 
205  if(!file1S)
206  {
207  cout << " Failed to open input root file" << endl;
208  return;
209  }
210 
211  recomassFG = (TH1D *)file1S->Get("Hunlk");
212  if(!recomassFG)
213  cout << "Failed to get Hunlk" << endl;
214 
215  recomassBG = (TH1D *)file1S->Get("HmxNorm");
216  if(!recomassBG)
217  cout << "Failed to get HmxNorm" << endl;
218 
219  TCanvas *cups = new TCanvas("cups","cups",5,5,1000,800);
220  recomassFG->Draw();
221  //recomassBG->Draw("same");
222 
223  TF1 *f1S;
224 
225  if(LL_fit)
226  {
227  // fit FG with CB+exp+BG12 using LL
228  f1S = new TF1("f1SLL",CBcalc_LL,1.3,maxmass,7);
229  }
230  else
231  {
232  // fit FG - BG with CB+exp using chisq
233  f1S = new TF1("f1S",CBcalc_exp,1.3,maxmass,7);
234  }
235 
236  // These are an error weighted average of the four free fits below in the comments
237  // Use these values to describe the J/psi line shape in fits at all cenralities
238 
239  f1S->SetParameter(0, -0.5); // alpha
240  f1S->SetParameter(1, 2.0); // n
241  f1S->SetParameter(2, 3.07); // mass
242  f1S->SetParameter(3, 0.11); // sigma
243  f1S->SetParameter(4, 850); // N
244  f1S->FixParameter(5, 0.0); // Nexp
245  f1S->FixParameter(6, -1.00); // slope
246  f1S->SetParNames("alpha1S","n1S","m1S","sigma1S","N1S","Nexp","slope");
247 
248  f1S->SetLineColor(kBlue);
249  f1S->SetLineWidth(3);
250  f1S->SetLineStyle(kDashed);
251 
252  TCanvas *csig = new TCanvas("csig","csig",5,5,1000,800);
253 
254  TF1 *fitResult= 0;
255  if(LL_fit)
256  {
257  // log likelihood fit with integer data option is "L"
258  //recomassFG->Fit(f1S,"LQ");
259  recomassFG->Fit(f1S,"L R");
260  recomassFG->DrawCopy();
261  fitResult = recomassFG->GetFunction("f1SLL");
262  }
263  else
264  {
265  recomassFG->Fit(f1S);
266  recomassFG->DrawCopy();
267  fitResult = recomassFG->GetFunction("CBcalc_exp");
268  }
269  if(fitResult)
270  {
271  cout << " chisquare " << fitResult->GetChisquare() << " NDF " << fitResult->GetNDF() << endl;
272  }
273  else
274  {
275  cout << "Failed to get fitResult" << endl;
276  return;
277  }
278 
279  // Plot the fit results
280  //======================
281  // Make a function that is only the CB lineshape and plot it
282  // We want CBcalc for this
283  TF1 *fcbonly = new TF1("fcbonly",CBcalc,1.3,4.0,5);
284  fcbonly->SetParameter(0,f1S->GetParameter(0));
285  fcbonly->SetParameter(1,f1S->GetParameter(1));
286  fcbonly->SetParameter(2,f1S->GetParameter(2));
287  fcbonly->SetParameter(3,f1S->GetParameter(3));
288  fcbonly->SetParameter(4,f1S->GetParameter(4));
289  fcbonly->SetLineColor(kMagenta);
290  fcbonly->Draw("same");
291 
292  // make a function that is only the exponential lineshape and plot it
293  TF1 *fexp = new TF1("fexp","[0]*exp([1]*x)",1.3,4.0);
294  fexp->SetParameter(0,f1S->GetParameter(5));
295  fexp->SetParameter(1,f1S->GetParameter(6));
296 
297  fexp->SetLineColor(kGreen);
298  fexp->SetLineWidth(3);
299  fexp->SetLineStyle(kDashed);
300  fexp->Draw("same");
301 
302  // Finally, plot the mixed event background
303  if(LL_fit)
304  {
305  recomassBG->SetLineColor(kRed);
306  recomassBG->DrawCopy("same");
307  }
308 
309  // calculate yield uncertainty from fit
310  double frac_yerror = f1S->GetParError(4) / f1S->GetParameter(4);
311 
312  /*
313  CB_yield[icent][ipt] = renorm * fcbonly->Integral(2.8,3.4) ;
314  CB_yield_err[icent][ipt] = CB_yield[icent][ipt] * frac_yerror;
315 
316  best_CB_yield[ialpha][in1][im1][isig] = CB_yield[icent][ipt];
317  best_CB_yield_err[ialpha][in1][im1][isig] = CB_yield_err[icent][ipt];
318  */
319 
320  // On a separate canvas, plot the FG - BG and the CB+exp
321  //========================================
322 
323  recomassFG->DrawCopy("pe");
324  recomassBG->DrawCopy("same");
325  fexp->Draw("same");
326  fcbonly->Draw("same");
327 
328 
329 }
330 
331 
332 
333