Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
CBfitter_jpsi_1state.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file CBfitter_jpsi_1state.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 <TGraphErrors.h>
9 #include <TLegend.h>
10 #include <TStyle.h>
11 #include <TROOT.h>
12 #include <TLatex.h>
13 
14 TH1D *recomass;
15 TH1D *recomassFG;
16 TH1D *recomassBG;
17 
18 Double_t CBcalc(Double_t *xx, Double_t *par)
19 {
20  // Crystal Ball lineshape plus exponential background
21  // Used for chisq fit to mixed event background data
22 
23  double f;
24  double x = xx[0];
25 
26  // The four CB parameters (alpha, n, x_mean, sigma) plus normalization (N) are:
27  double alpha = par[0];
28  double n = par[1];
29  double x_mean = par[2];
30  double sigma = par[3];
31  double N = par[4];
32 
33  // we need:
34  double A = pow( (n/TMath::Abs(alpha)),n) * exp(-pow(alpha,2)/2.0);
35  double B = n/TMath::Abs(alpha) - TMath::Abs(alpha);
36 
37  // The Crystal Ball function is:
38  if( (x-x_mean)/sigma > -alpha)
39  {
40  f = N * exp( -pow(x-x_mean,2) / (2.0*pow(sigma,2)));
41  }
42  else
43  {
44  f = N * A * pow(B - (x-x_mean)/sigma, -n);
45  }
46 
47  return f;
48 }
49 
50 Double_t CBcalc_exp(Double_t *xx, Double_t *par)
51 {
52  // Crystal Ball lineshape plus exponential background
53  // Used for chisq fit to mixed event background data
54 
55  double f;
56  double x = xx[0];
57 
58  // The four CB parameters (alpha, n, x_mean, sigma) plus normalization (N) are:
59  double alpha = par[0];
60  double n = par[1];
61  double x_mean = par[2];
62  double sigma = par[3];
63  double N = par[4];
64  // and we add an exponential background
65  double Nexp = par[5];
66  double slope = par[6];
67 
68  // we need:
69  double A = pow( (n/TMath::Abs(alpha)),n) * exp(-pow(alpha,2)/2.0);
70  double B = n/TMath::Abs(alpha) - TMath::Abs(alpha);
71 
72  // The Crystal Ball function is:
73  if( (x-x_mean)/sigma > -alpha)
74  {
75  f = N * exp( -pow(x-x_mean,2) / (2.0*pow(sigma,2)));
76  }
77  else
78  {
79  f = N * A * pow(B - (x-x_mean)/sigma, -n);
80  }
81 
82  f += Nexp * exp(slope * x);
83 
84  return f;
85 }
86 
87 Double_t CBcalc_LL(Double_t *xx, Double_t *par)
88 {
89  // Crystal Ball lineshape plus exponential background plus mixed event combinatorial background
90  // Used for LL fit to foreground data
91  double f;
92  double x = xx[0];
93 
94  // The four CB parameters (alpha, n, x_mean, sigma) plus normalization (N) are:
95  double alpha = par[0];
96  double n = par[1];
97  double x_mean = par[2];
98  double sigma = par[3];
99  double N = par[4];
100  // and we add an exponential background
101  double Nexp = par[5];
102  double slope = par[6];
103 
104  /*
105  TF1 *fcbexp = new TF1("fcbexp",CBcalc_exp,2.0,4.0,7);
106  fcbexp->SetParameter(0,alpha);
107  fcbexp->SetParameter(1,n);
108  fcbexp->SetParameter(2,x_mean);
109  fcbexp->SetParameter(3,sigma);
110  fcbexp->SetParameter(4, N);
111  fcbexp->SetParameter(5,Nexp);
112  fcbexp->SetParameter(6,slope);
113  f = fcbexp->Eval(x);
114  */
115 
116  // we need:
117  double A = pow( (n/TMath::Abs(alpha)),n) * exp(-pow(alpha,2)/2.0);
118  double B = n/TMath::Abs(alpha) - TMath::Abs(alpha);
119 
120  // The Crystal Ball function is:
121  if( (x-x_mean)/sigma > -alpha)
122  {
123  f = N * exp( -pow(x-x_mean,2) / (2.0*pow(sigma,2)));
124  }
125  else
126  {
127  f = N * A * pow(B - (x-x_mean)/sigma, -n);
128  }
129 
130  f += Nexp * exp(slope * x);
131 
132  // Add the mixed event background to the fit function
133  int ibin = recomassBG->FindBin(x);
134 
135  double bg12 = recomassBG->GetBinContent(ibin);
136  f += bg12;
137 
138 
139  return f;
140 }
141 
142 Double_t Upscalc(Double_t *xx, Double_t *par)
143 {
144  double x = xx[0];
145 
146  // The input parameters are:
147  double alpha1s = par[0];
148  double n1s = par[1];
149  double sigma1s = par[2];
150  double m1s = par[3];
151  double m2s = par[4];
152  double m3s = par[5];
153  double N1s = par[6];
154  double N2s = par[7];
155  double N3s = par[8];
156  double Nexp = par[9];
157  double slope = par[10];
158 
159  // Construct the Y(1S) CB shape
160 
161  TF1 *f1 = new TF1("f1",CBcalc,7,11,5);
162  f1->SetParameter(0,alpha1s);
163  f1->SetParameter(1,n1s);
164  f1->SetParameter(2,m1s);
165  f1->SetParameter(3,sigma1s);
166  f1->SetParameter(4,N1s);
167 
168  TF1 *f2 = new TF1("f2",CBcalc,7,11,5);
169  f2->SetParameter(0,alpha1s);
170  f2->SetParameter(1,n1s);
171  f2->SetParameter(2,m2s);
172  f2->SetParameter(3,sigma1s);
173  f2->SetParameter(4,N2s);
174 
175  TF1 *f3 = new TF1("f3",CBcalc,7,11,5);
176  f3->SetParameter(0,alpha1s);
177  f3->SetParameter(1,n1s);
178  f3->SetParameter(2,m3s);
179  f3->SetParameter(3,sigma1s);
180  f3->SetParameter(4,N3s);
181 
182  TF1 *fexp = new TF1("fexp","[0]*exp([1]*x)",7,11);
183  fexp->SetParameter(0,Nexp);
184  fexp->SetParameter(1,slope);
185 
186 
187  double mass = f1->Eval(x) + f2->Eval(x) + f3->Eval(x) + fexp->Eval(x);
188 
189  return mass;
190 }
191 
192 void CBfitter_jpsi_1state()
193 {
194  gROOT->SetStyle("Plain");
195  gStyle->SetOptStat(0);
196  gStyle->SetOptFit(0);
197  gStyle->SetOptTitle(0);
198 
199  bool LL_fit = true;
200 
201  double maxmass = 3.7;
202  double minmass = 2.4;
203 
204  TFile *file1S;
205  file1S = new TFile("Runcuts_DIMU_MUTRONLY_Run15pp200_JPSISIM_NOEMBED_C_PTBIN0_0_NONE.root");
206  if(!file1S)
207  {
208  cout << " Failed to open input root file" << endl;
209  return;
210  }
211 
212  TH2D *recomass2D[2];
213  TH1D *recomass[2][4];
214  for(int iarm = 0; iarm < 2; ++iarm)
215  {
216  // Get the 2D histos
217  char hname[500];
218  sprintf(hname, "mass_y_same_arm%i_chg0", iarm);
219  recomass2D[iarm] = (TH2D *)file1S->Get(hname);
220 
221  if(!recomass2D[iarm])
222  cout << "Failed to get 2D histogram " << endl;
223 
224  // and project the four rapidity bins
225  for(int irap =0; irap < 4; ++irap)
226  {
227  char h1dname[500];
228  sprintf(h1dname, "recomass%i%i", iarm, irap);
229  recomass[iarm][irap] = recomass2D[iarm]->ProjectionX(h1dname, irap+1, irap+1);
230  }
231  }
232 
233  TF1 *f1S[2][4];
234  TCanvas *cups = new TCanvas("cups","cups",5,5,1600,800);
235  cups->Divide(4,2);
236  for(int iarm = 0; iarm<2; ++iarm)
237  {
238  for(int irap = 0; irap < 4; ++irap)
239  {
240  cout << "Fit for iarm = " << iarm << " irap = " << irap << endl;
241 
242  cups->cd(irap + iarm*4 + 1);
243  recomass[iarm][irap]->Draw();
244 
245  char fname[500];
246  sprintf(fname, "f1S%i%i", iarm, irap);
247  f1S[iarm][irap] = new TF1(fname,CBcalc, minmass, maxmass, 5);
248 
249  f1S[iarm][irap]->SetParameter(0, 1.5); // alpha
250  f1S[iarm][irap]->SetParameter(1, 100.0); // n
251  f1S[iarm][irap]->SetParameter(2, 3.12); // mass
252  f1S[iarm][irap]->SetParameter(3, 0.13); // sigma
253  f1S[iarm][irap]->SetParameter(4, 300); // N
254  f1S[iarm][irap]->SetParNames("alpha1S","n1S","m1S","sigma1S","N1S");
255  f1S[iarm][irap]->SetLineColor(kBlue);
256  f1S[iarm][irap]->SetLineWidth(3);
257  f1S[iarm][irap]->SetLineStyle(kDashed);
258 
259  recomass[iarm][irap]->GetXaxis()->SetRangeUser(2.2,3.8);
260  recomass[iarm][irap]->Fit(f1S[iarm][irap]);
261  //recomass[iarm][irap]->Fit(f1S[iarm][irap],"L");
262  recomass[iarm][irap]->DrawCopy();
263  f1S[iarm][irap]->Draw("same");
264  }
265  }
266 
267  // Get the mass and width so we can plot them
268 
269  TCanvas *cgr = new TCanvas("cgr", "cgr", 5, 5, 1600, 800);
270  cgr->Divide(2,2);
271 
272  TLatex *lab[2];
273  double rap[2][4] = {-2.075, -1.825, -1.575, -1.325, 1.325, 1.575, 1.825, 2.075};
274  double rap_err[4] = {0};
275  //TGraph *grmean[2];
276  //TGraph *grsigma[2];
277  TGraphErrors *grmean[2];
278  TGraphErrors *grsigma[2];
279  double mean[2][4];
280  double sigma[2][4];
281  double mean_err[2][4];
282  double sigma_err[2][4];
283  for(int iarm = 0; iarm < 2; ++iarm)
284  {
285  for(int irap = 0; irap < 4; ++irap)
286  {
287  mean[iarm][irap] = f1S[iarm][irap]->GetParameter(2);
288  mean_err[iarm][irap] = f1S[iarm][irap]->GetParError(2);
289  sigma[iarm][irap] = f1S[iarm][irap]->GetParameter(3);
290  sigma_err[iarm][irap] = f1S[iarm][irap]->GetParError(3);
291  }
292  grmean[iarm] = new TGraphErrors(4, rap[iarm], mean[iarm], rap_err, mean_err[iarm]);
293  grsigma[iarm] = new TGraphErrors(4, rap[iarm], sigma[iarm], rap_err, sigma_err[iarm]);
294 
295  cgr->cd(1+iarm);
296  grmean[iarm]->SetMarkerStyle(20);
297  grmean[iarm]->SetMarkerSize(2);
298  grmean[iarm]->GetYaxis()->SetTitle("mass (GeV/c^{2})");
299  grmean[iarm]->GetYaxis()->SetTitleSize(0.05);
300  grmean[iarm]->GetXaxis()->SetTitle("y");
301  grmean[iarm]->GetXaxis()->SetTitleSize(0.05);
302  grmean[iarm]->Draw();
303 
304  char arm[500];
305  sprintf(arm,"arm %i",iarm);
306  lab[iarm] = new TLatex(0.2, 0.350, arm);
307  lab[iarm]->SetNDC();
308  lab[iarm]->SetTextSize(0.1);
309  lab[iarm]->Draw();
310 
311  cgr->cd(3+iarm);
312  grsigma[iarm]->SetMarkerStyle(20);
313  grsigma[iarm]->SetMarkerSize(2);
314  grsigma[iarm]->GetYaxis()->SetTitle("width (GeV/c^{2})");
315  grsigma[iarm]->GetYaxis()->SetTitleSize(0.05);
316  grsigma[iarm]->GetXaxis()->SetTitle("y");
317  grsigma[iarm]->GetYaxis()->SetTitleSize(0.05);
318  grsigma[iarm]->Draw();
319  lab[iarm]->Draw();
320  }
321 
322 
323 
324 }
325 
326 
327 
328 
329