Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
CBfitter_1state.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file CBfitter_1state.C
1 #include <TF1.h>
2 #include <TMath.h>
3 #include <TFile.h>
4 #include <TCanvas.h>
5 #include <TH1.h>
6 #include <TGraph.h>
7 #include <TLegend.h>
8 #include <TStyle.h>
9 #include <TROOT.h>
10 #include <TLatex.h>
11 
12 Double_t CBcalc(Double_t *xx, Double_t *par)
13 {
14  double f;
15  double x = xx[0];
16 
17  // The four parameters (alpha, n, x_mean, sigma) plus normalization (N) are:
18  double alpha = par[0];
19  double n = par[1];
20  double x_mean = par[2];
21  double sigma = par[3];
22  double N = par[4];
23 
24  // we need:
25  double A = pow( (n/TMath::Abs(alpha)),n) * exp(-pow(alpha,2)/2.0);
26  double B = n/TMath::Abs(alpha) - TMath::Abs(alpha);
27 
28  // The Crystal Ball function is:
29  if( (x-x_mean)/sigma > -alpha)
30  {
31  f = N * exp( -pow(x-x_mean,2) / (2.0*pow(sigma,2)));
32  }
33  else
34  {
35  f = N * A * pow(B - (x-x_mean)/sigma, -n);
36  }
37 
38  double Nexp = par[5];
39  double slope = par[6];
40 
41  TF1 *fexp = new TF1("fexp","[0]*exp([1]*x)",7,11);
42  fexp->SetParameter(0,Nexp);
43  fexp->SetParameter(1,slope);
44 
45  f = f + fexp->Eval(x);
46  return f;
47 }
48 
49 Double_t Upscalc(Double_t *xx, Double_t *par)
50 {
51  double x = xx[0];
52 
53  // The input parameters are:
54  double alpha1s = par[0];
55  double n1s = par[1];
56  double sigma1s = par[2];
57  double m1s = par[3];
58  double m2s = par[4];
59  double m3s = par[5];
60  double N1s = par[6];
61  double N2s = par[7];
62  double N3s = par[8];
63  double Nexp = par[9];
64  double slope = par[10];
65 
66  // Construct the Y(1S) CB shape
67 
68  TF1 *f1 = new TF1("f1",CBcalc,7,11,5);
69  f1->SetParameter(0,alpha1s);
70  f1->SetParameter(1,n1s);
71  f1->SetParameter(2,m1s);
72  f1->SetParameter(3,sigma1s);
73  f1->SetParameter(4,N1s);
74 
75  TF1 *f2 = new TF1("f2",CBcalc,7,11,5);
76  f2->SetParameter(0,alpha1s);
77  f2->SetParameter(1,n1s);
78  f2->SetParameter(2,m2s);
79  f2->SetParameter(3,sigma1s);
80  f2->SetParameter(4,N2s);
81 
82  TF1 *f3 = new TF1("f3",CBcalc,7,11,5);
83  f3->SetParameter(0,alpha1s);
84  f3->SetParameter(1,n1s);
85  f3->SetParameter(2,m3s);
86  f3->SetParameter(3,sigma1s);
87  f3->SetParameter(4,N3s);
88 
89  TF1 *fexp = new TF1("fexp","[0]*exp([1]*x)",7,11);
90  fexp->SetParameter(0,Nexp);
91  fexp->SetParameter(1,slope);
92 
93 
94 
95  double mass = f1->Eval(x) + f2->Eval(x) + f3->Eval(x) + fexp->Eval(x);
96 
97  return mass;
98 }
99 
100 
101 void CBfitter_1state()
102 {
103  gROOT->SetStyle("Plain");
104  gStyle->SetOptStat(0);
105  gStyle->SetOptFit(0);
106  gStyle->SetOptTitle(0);
107 
108  TFile *file1S, *file2S, *file3S;
109  TH1 *recomass1S, *recomass2S, *recomass3S;
110 
111  bool draw_gauss = false;
112 
113  // if all are false, do AuAu 0-10%
114  bool pp = false;
115  bool pAu = false;
116  bool AuAu_20pc = false;
117  bool AuAu_10pc = true;
118 
119  //bool do_subtracted = true;
120 
121  file1S = new TFile("root_files/ntp_quarkonium_out.root");
122 
123  if(!file1S)
124  {
125  cout << "Failed to open file " << endl;
126  exit(1);
127  }
128 
129  recomass1S = (TH1 *)file1S->Get("recomass0");
130 
131  TH1 *recomass = (TH1*)recomass1S->Clone("recomass1");
132  if(!recomass)
133  {
134  cout << "Failed to get recomass histogram " << endl;
135  exit(1);
136  }
137 
138  recomass->GetXaxis()->SetTitle("invariant mass (GeV/c^{2})");
139  recomass->Sumw2();
140 
141  int nrebin = 1; // set to 2 to match background histo binning, but this worsens resolution slightly. Use 1 for signal fit
142  //int nrebin = 2; // set to 2 to match background histo binning, but this worsens resolution slightly. Use 1 for signal fit
143  recomass->Rebin(nrebin);
144 
145  TCanvas *cups = new TCanvas("cups","cups",5,5,800,800);
146  recomass->SetTitle("Y(1S,2S,3S) #rightarrow e^{+}e^{-}");
147  recomass->SetMarkerStyle(20);
148  recomass->SetMarkerSize(1);
149  recomass->SetLineStyle(kSolid);
150  recomass->SetLineWidth(2);
151  // recomass->SetMaximum(700);
152  recomass->DrawCopy("p");
153 
154 
155  TF1 *f1S = new TF1("f1S",CBcalc,7,11,7);
156  f1S->SetParameter(0, 1.0); // alpha
157  f1S->SetParameter(1, 1.0); // n
158  f1S->SetParameter(2, 9.46); // xmean
159  f1S->SetParameter(3, 0.08); // sigma
160  f1S->SetParameter(4, 2000.0); // N
161  //f1S->SetParameter(4, 200.0); // N
162  f1S->SetParameter(5,3.5);
163  f1S->SetParameter(6,0.05);
164 
165  f1S->SetParNames("alpha1S","n1S","m1S","sigma1S","N1S");
166  f1S->SetLineColor(kBlue);
167  f1S->SetLineWidth(3);
168  f1S->SetLineStyle(kDashed);
169 
170  recomass->Fit(f1S);
171  f1S->Draw("same");
172  cout << "f1S pars " << f1S->GetParameter(3) << " " << f1S->GetParError(3) << endl;
173 
174  char resstr[500];
175  sprintf(resstr,"#sigma_{1S} = %.1f #pm %.1f MeV", f1S->GetParameter(3)*1000, f1S->GetParError(3)*1000);
176  TLatex *res = new TLatex(0.13,0.55,resstr);
177  res->SetNDC();
178  res->SetTextSize(0.05);
179  res->Draw();
180 
181 
182  double binw = recomass->GetBinWidth(1);
183  double renorm = 1.0/binw; // (1 / (bin_width of data in GeV) )
184  cout << "renorm = " << renorm << endl;
185 
186  cout << "Area of f1S is " << renorm * f1S->Integral(7,11) << endl;
187 
188  // Extract ratio of yield in central gaussian to total
189 
190  TF1 *fgauss = new TF1("fgauss","gaus(0)",7,11);
191  fgauss->SetParameter(0, f1S->GetParameter(4));
192  fgauss->SetParameter(1, f1S->GetParameter(2));
193  fgauss->SetParameter(2, f1S->GetParameter(3));
194  fgauss->SetLineColor(kRed);
195  if(draw_gauss) fgauss->Draw("same");
196 
197  // calculate fraction of yield in gaussian
198  double area_fgauss = fgauss->Integral(7,11) * renorm;
199  double area_f1S = f1S->Integral(7,11) * renorm;
200  double fraction = area_fgauss / area_f1S;
201 
202 
203  cout << "Parameters of fgauss = " << fgauss->GetParameter(0) << " " << fgauss->GetParameter(1) << " " << fgauss->GetParameter(2) << " Area of fgauss is " << renorm * fgauss->Integral(7,11) << " fraction in fgauss " << area_fgauss / area_f1S << endl;
204 
205  char labfrac[500];
206  sprintf(labfrac, "Gauss fraction %.2f", fraction);
207  TLatex *lab = new TLatex(0.13,0.75,labfrac);
208  lab->SetNDC();
209  lab->SetTextSize(0.05);
210  if(draw_gauss) lab->Draw();
211 
212 
213  /*
214  TLatex *lab;
215 
216  if(pp)
217  lab = new TLatex(0.15,0.75,"p+p, 10 weeks");
218  else if(pAu)
219  lab = new TLatex(0.20,0.75,"#splitline{p+Au 0-20%}{10 weeks}");
220  else if(AuAu_20pc)
221  lab = new TLatex(0.20,0.75,"#splitline{Au+Au 0-20%}{20B events}");
222  else
223  lab = new TLatex(0.20,0.75,"#splitline{Au+Au 0-10%}{10B events}");
224 
225  lab->SetNDC();
226  lab->SetTextSize(0.05);
227  lab->Draw();
228 
229  char resstr[500];
230  sprintf(resstr,"#sigma_{1S} = %.0f #pm %.0f MeV",sigma1s,sigma1s_err);
231  TLatex *res = new TLatex(0.13,0.55,resstr);
232  res->SetNDC();
233  res->SetTextSize(0.05);
234  res->Draw();
235 
236  */
237  /*
238  if (do_subtracted)
239  {
240  //========================================================
241  // Now make the background_subtracted spectrum and show the
242  // CB fits on it
243 
244  //=======================================
245  // Read in the subtracted mass spectrum
246 
247  // There is only one file to read, it was written by upsilon_mass_plus_background(_pp).C
248  TFile *file_all;
249  if(pp)
250  file_all = new TFile("SplusBtotal_pp_to80.root");
251  else
252  file_all = new TFile("SplusBtotal_10Bevts_0_10pc_rej90_eff70.root");
253 
254  TH1 *hsplusb, *hsub, *comb_bgd, *corr_bgd;
255  hsplusb = (TH1 *)file_all->Get("hSplusB");
256  hsub = (TH1 *)file_all->Get("hsub");
257  corr_bgd = (TH1 *)file_all->Get("hdy_new");
258  if(!pp)
259  comb_bgd = (TH1 *)file_all->Get("hfake_new_bgd");
260 
261  // Fit the correlated background with a function
262 
263  //TCanvas *c = new TCanvas("c","c",30,30,900,600);
264  TCanvas *c = new TCanvas("c","c",30,30,800,800);
265 
266  TF1 *fbgd = new TF1("fbgd","[0]*exp([1]+[2]*x)",7,11);
267  fbgd->SetLineColor(kRed);
268  fbgd->SetLineWidth(3);
269  fbgd->SetLineStyle(kDashed);
270  corr_bgd->Fit(fbgd);
271 
272  hsub->SetMarkerStyle(20);
273  hsub->SetMarkerSize(1);
274  hsub->GetXaxis()->SetRangeUser(7.8,11.0);
275  hsub->Draw("p");
276  f1S->Draw("same");
277  f2S->Draw("same");
278  f3S->Draw("same");
279  //fexp->Draw("same");
280  fbgd->Draw("same");
281 
282  // Now want to make the total of all of the fit functions
283 
284  static const int NSTEP = 200;
285  double xm[NSTEP];
286  double ym[NSTEP];
287 
288  double step = 4.0 / (double) NSTEP;
289  for(int i=0;i<NSTEP;i++)
290  {
291  double x = 7.0 + (double) i * step;
292  double y = f1S->Eval(x);
293  y += f2S->Eval(x);
294  y += f3S->Eval(x);
295  y += fbgd->Eval(x);
296  y += fexp->Eval(x);
297 
298  xm[i] = x;
299  ym[i] = y;
300  }
301  TGraph *gtot = new TGraph(NSTEP,xm,ym);
302  gtot->SetLineWidth(3);
303  gtot->SetLineStyle(kSolid);
304  gtot->SetLineColor(kBlue);
305  gtot->Draw("l");
306 
307  TLegend *leg = new TLegend(0.62,0.58,0.9,0.88,"");
308  leg->SetBorderSize(0);
309  leg->SetFillColor(0);
310  leg->SetFillStyle(0);
311  leg->AddEntry(hsub,"++,-- subtracted","p");
312  leg->AddEntry(f1S,"Y(1S)","l");
313  leg->AddEntry(f2S,"Y(2S)","l");
314  leg->AddEntry(f3S,"Y(3S)","l");
315  leg->AddEntry(fbgd,"correlated bkg","l");
316 
317  leg->Draw();
318 
319  TLatex *lab2;
320 
321  if(pp)
322  lab2 = new TLatex(0.15,0.75,"p+p, 10 weeks");
323  else
324  lab2 = new TLatex(0.20,0.75,"#splitline{Au+Au 0-10%}{10B events}");
325 
326  lab2->SetNDC();
327  lab2->SetTextSize(0.05);
328  lab2->Draw();
329 
330  }
331  */
332 }