Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
CBfitter.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file CBfitter.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  return f;
39 }
40 
41 Double_t Upscalc(Double_t *xx, Double_t *par)
42 {
43  double x = xx[0];
44 
45  // The input parameters are:
46  double alpha1s = par[0];
47  double n1s = par[1];
48  double sigma1s = par[2];
49  double m1s = par[3];
50  double m2s = par[4];
51  double m3s = par[5];
52  double N1s = par[6];
53  double N2s = par[7];
54  double N3s = par[8];
55  double Nexp = par[9];
56  double slope = par[10];
57 
58  // Construct the Y(1S) CB shape
59 
60  TF1 *f1 = new TF1("f1",CBcalc,7,11,5);
61  f1->SetParameter(0,alpha1s);
62  f1->SetParameter(1,n1s);
63  f1->SetParameter(2,m1s);
64  f1->SetParameter(3,sigma1s);
65  f1->SetParameter(4,N1s);
66 
67  TF1 *f2 = new TF1("f2",CBcalc,7,11,5);
68  f2->SetParameter(0,alpha1s);
69  f2->SetParameter(1,n1s);
70  f2->SetParameter(2,m2s);
71  f2->SetParameter(3,sigma1s);
72  f2->SetParameter(4,N2s);
73 
74  TF1 *f3 = new TF1("f3",CBcalc,7,11,5);
75  f3->SetParameter(0,alpha1s);
76  f3->SetParameter(1,n1s);
77  f3->SetParameter(2,m3s);
78  f3->SetParameter(3,sigma1s);
79  f3->SetParameter(4,N3s);
80 
81  TF1 *fexp = new TF1("fexp","[0]*exp([1]*x)",7,11);
82  fexp->SetParameter(0,Nexp);
83  fexp->SetParameter(1,slope);
84 
85  double mass = f1->Eval(x) + f2->Eval(x) + f3->Eval(x) + fexp->Eval(x);
86 
87  return mass;
88 }
89 
90 
91 void CBfitter()
92 {
93  gROOT->SetStyle("Plain");
94  gStyle->SetOptStat(0);
95  gStyle->SetOptTitle(1);
96 
97  bool ups1s = true;
98  bool ups2s = true;
99  bool ups3s = true;
100 
101  TFile *file1S, *file2S, *file3S;
102  TH1 *recomass1S, *recomass2S, *recomass3S;
103 
104  if(ups1s)
105  {
106  file1S = new TFile("root_files/unsuppressed_auau_24B_ntp_quarkonium_out.root");
107  recomass1S = (TH1 *)file1S->Get("recomass0");
108  }
109 
110  if(ups2s)
111  {
112  //file2S = new TFile(" root_files/ups2s_80ns_100pions_pp.root ");
113  recomass2S = (TH1 *)file1S->Get("recomass1");
114  }
115  if(ups3s)
116  {
117  //file3S = new TFile(" root_files/ups3s_80ns_100pions_pp.root");
118  recomass3S = (TH1 *)file1S->Get("recomass2");
119  }
120 
121  TFile *massout = new TFile("upsilon_mass_histos_out.root","recreate");
122 
123  TH1 *recomass = (TH1*)recomass1S->Clone("recomass");
124  recomass->Sumw2();
125  if(ups2s)
126  recomass->Add(recomass2S);
127  if(ups3s)
128  recomass->Add(recomass3S);
129  int nrebin = 1; // a rebin of 2 is needed to match background histo binning, use 1 if want best resolution
130  recomass->Rebin(nrebin);
131 
132  //TCanvas *cups = new TCanvas("cups","cups",5,5,800,600);
133  TCanvas *cups = new TCanvas("cups","cups",5,5,800,800);
134  recomass->SetTitle("Y(1S,2S,3S) #rightarrow e^{+}e^{-}");
135  recomass->GetYaxis()->SetNdivisions(205);
136  recomass->GetXaxis()->SetTitle("Invariant mass (GeV/c^{2})");
137  recomass->SetMarkerStyle(20);
138  recomass->SetMarkerSize(1);
139  recomass->SetLineStyle(kSolid);
140  recomass->SetLineWidth(2);
141  //recomass->SetMaximum(810);
142  recomass->DrawCopy("p");
143 
144 
145  TF1 *fitups = new TF1("fitups",Upscalc,7,11,11);
146  fitups->SetParameter(0, 1.05);
147  fitups->SetParameter(1, 1.5);
148  fitups->SetParameter(2, 0.09);
149  fitups->SetParameter(3, 9.44);
150  fitups->SetParameter(4, 10.0);
151  fitups->SetParameter(5, 10.3);
152  fitups->SetParameter(6, 900.0);
153  fitups->SetParameter(7, 200.0);
154  fitups->SetParameter(8, 100.0);
155  //fitups->SetParameter(9, 10.0);
156  fitups->FixParameter(9, 0.0);
157  fitups->FixParameter(10, -0.2);
158  fitups->SetLineColor(kBlack);
159  fitups->SetParNames("alpha1S","n1S","sigma1S","m1S","m2S","m3S","N1S","N2S","N3S","Nexp","expSlope");
160 
161  recomass->Fit(fitups);
162  //fitups->Draw("same");
163 
164  double sigma1s = 1000.0 * fitups->GetParameter(2);
165  double sigma1s_err = 1000.0 * fitups->GetParError(2);
166  cout << " sigma1s " << sigma1s << " par2 " << fitups->GetParameter(2) << endl;
167  // Now draw the individual states
168 
169  TF1 *f1S = new TF1("f1S",CBcalc,7,11,5);
170  f1S->SetParameter(0,fitups->GetParameter(0));
171  f1S->SetParameter(1,fitups->GetParameter(1));
172  f1S->SetParameter(2,fitups->GetParameter(3));
173  f1S->SetParameter(3,fitups->GetParameter(2));
174  f1S->SetParameter(4,fitups->GetParameter(6));
175  f1S->SetParNames("alpha1S","n1S","m1S","sigma1S","N1S");
176  f1S->SetLineColor(kBlue);
177  f1S->SetLineWidth(3);
178  f1S->SetLineStyle(kDashed);
179 
180  f1S->Draw("same");
181 
182  TF1 *f2S = new TF1("f2S",CBcalc,7,11,5);
183  f2S->SetParameter(0,fitups->GetParameter(0));
184  f2S->SetParameter(1,fitups->GetParameter(1));
185  f2S->SetParameter(2,fitups->GetParameter(4));
186  f2S->SetParameter(3,fitups->GetParameter(2));
187  f2S->SetParameter(4,fitups->GetParameter(7));
188  f2S->SetParNames("alpha1S","n1S","m2S","sigma1S","N2S");
189  f2S->SetLineColor(kRed);
190  f2S->SetLineWidth(3);
191  f2S->SetLineStyle(kDashed);
192 
193  f2S->Draw("same");
194 
195  TF1 *f3S = new TF1("f3S",CBcalc,7,11,5);
196  f3S->SetParameter(0,fitups->GetParameter(0));
197  f3S->SetParameter(1,fitups->GetParameter(1));
198  f3S->SetParameter(2,fitups->GetParameter(5));
199  f3S->SetParameter(3,fitups->GetParameter(2));
200  f3S->SetParameter(4,fitups->GetParameter(8));
201  f3S->SetParNames("alpha1S","n1S","m3S","sigma1S","N3S");
202  f3S->SetLineColor(kMagenta);
203  f3S->SetLineWidth(3);
204  f3S->SetLineStyle(kDashed);
205 
206  f3S->Draw("same");
207 
208  TF1 *fexp = new TF1("fexp","[0]*exp([1]*x)",7,11);
209  fexp->SetParameter(0,fitups->GetParameter(9));
210  fexp->SetParameter(1,fitups->GetParameter(10));
211 
212  fexp->SetLineColor(kGreen);
213  fexp->SetLineWidth(3);
214  fexp->SetLineStyle(kDashed);
215  fexp->Draw("same");
216 
217  double binw = recomass->GetBinWidth(1);
218  double renorm = 1.0/binw; // (1 / (bin_width of data in GeV) )
219  //double renorm = 25.0; // (1 / (bin_width of data in GeV) )
220  cout << "renorm = " << renorm << endl;
221 
222  cout << "Area of f1S is " << renorm * f1S->Integral(7,11) << endl;
223  cout << "Area of f2S is " << renorm * f2S->Integral(7,11) << endl;
224  cout << "Area of f3S is " << renorm * f3S->Integral(7,11) << endl;
225  cout << "Area of fexp is " << renorm * fexp->Integral(7,11) << endl;
226 
227  bool pp = true;
228  bool AuAu_20pc = false;
229 
230  bool do_subtracted = false;
231 
232  TLatex *lab;
233 
234  if(pp)
235  lab = new TLatex(0.15,0.75,"#splitline{p+p, 197 pb^{-1}}{signal only}");
236  else if(AuAu_20pc)
237  lab = new TLatex(0.20,0.75,"#splitline{Au+Au 0-20%}{10B events}");
238  else
239  lab = new TLatex(0.20,0.75,"#splitline{Au+Au MB}{50B events}");
240 
241  lab->SetNDC();
242  lab->SetTextSize(0.05);
243  //lab->Draw();
244 
245  char resstr[500];
246  sprintf(resstr,"#sigma_{1S} = %.0f #pm %.1f MeV",sigma1s,sigma1s_err);
247  TLatex *res = new TLatex(0.13,0.55,resstr);
248  res->SetNDC();
249  res->SetTextSize(0.05);
250  //res->Draw();
251 
252  recomass->Write();
253  f1S->Write();
254  f2S->Write();
255  f3S->Write();
256 
257  if (do_subtracted)
258  {
259  //========================================================
260  // Now make the background_subtracted spectrum and show the
261  // CB fits on it
262 
263  //=======================================
264  // Read in the subtracted mass spectrum
265 
266  // There is only one file to read, it was written by upsilon_mass_plus_background(_pp).C
267  TFile *file_all;
268  if(pp)
269  file_all = new TFile("SplusBtotal_pp_to80.root");
270  else
271  file_all = new TFile("SplusBtotal_rej90_to80.root");
272 
273  TH1 *hsplusb, *hsub, *comb_bgd, *corr_bgd;
274  hsplusb = (TH1 *)file_all->Get("hSplusB");
275  hsub = (TH1 *)file_all->Get("hsub");
276  corr_bgd = (TH1 *)file_all->Get("hdy_new");
277  if(!pp)
278  comb_bgd = (TH1 *)file_all->Get("hfake_new_bgd");
279 
280  // Fit the correlated background with a function
281 
282  //TCanvas *c = new TCanvas("c","c",30,30,900,600);
283  TCanvas *c = new TCanvas("c","c",30,30,800,800);
284 
285  TF1 *fbgd = new TF1("fbgd","[0]*exp([1]+[2]*x)",7,11);
286  fbgd->SetLineColor(kRed);
287  fbgd->SetLineWidth(3);
288  fbgd->SetLineStyle(kDashed);
289  corr_bgd->Fit(fbgd);
290 
291  hsub->SetMarkerStyle(20);
292  hsub->SetMarkerSize(1);
293  hsub->Draw("p");
294  f1S->Draw("same");
295  f2S->Draw("same");
296  f3S->Draw("same");
297  //fexp->Draw("same");
298  fbgd->Draw("same");
299 
300  // Now want to make the total of all of the fit functions
301 
302  static const int NSTEP = 200;
303  double xm[NSTEP];
304  double ym[NSTEP];
305 
306  double step = 4.0 / (double) NSTEP;
307  for(int i=0;i<NSTEP;i++)
308  {
309  double x = 7.0 + (double) i * step;
310  double y = f1S->Eval(x);
311  y += f2S->Eval(x);
312  y += f3S->Eval(x);
313  y += fbgd->Eval(x);
314  y += fexp->Eval(x);
315 
316  xm[i] = x;
317  ym[i] = y;
318  }
319  TGraph *gtot = new TGraph(NSTEP,xm,ym);
320  gtot->SetLineWidth(3);
321  gtot->SetLineStyle(kSolid);
322  gtot->SetLineColor(kBlue);
323  gtot->Draw("l");
324 
325  TLegend *leg = new TLegend(0.62,0.58,0.9,0.88,"");
326  leg->SetBorderSize(0);
327  leg->SetFillColor(0);
328  leg->SetFillStyle(0);
329  leg->AddEntry(hsub,"++,-- subtracted","p");
330  leg->AddEntry(f1S,"Y(1S)","l");
331  leg->AddEntry(f2S,"Y(2S)","l");
332  leg->AddEntry(f3S,"Y(3S)","l");
333  leg->AddEntry(fbgd,"correlated bkg","l");
334 
335  leg->Draw();
336 
337  TLatex *lab2;
338 
339  if(pp)
340  lab2 = new TLatex(0.15,0.75,"p+p, 10 weeks");
341  else
342  lab2 = new TLatex(0.20,0.75,"#splitline{Au+Au 0-20%}{10B events}");
343 
344  lab2->SetNDC();
345  lab2->SetTextSize(0.05);
346  lab2->Draw();
347 
348  }
349 }