Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
plotembed_vscent.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file plotembed_vscent.C
1 #include <TNtuple.h>
2 #include <TF1.h>
3 #include <TLine.h>
4 #include "sPHENIXStyle/sPhenixStyle.C"
5 
6 #include <iostream>
7 #include <fstream>
8 
9 //----- MY STUFF -------------------------------------------------
10 double CBFunction(double *x, double *p)
11 {
12  double norm = p[0];
13  double alpha = p[1]; // tail start
14  double n = p[2]; // tail shape
15  double mu = p[3]; // upsilon mass
16  double sigma = p[4]; // upsilon width
17 
18  double A = pow(n/fabs(alpha),n)*TMath::Exp(-pow(fabs(alpha),2)/2.);
19  double B = n/fabs(alpha) - fabs(alpha);
20  double k = (x[0]-mu)/sigma;
21 
22  double val;
23  if( k > -alpha )
24  val = norm*TMath::Exp(-0.5*pow(k,2));
25  else
26  val = norm*A*pow(B-k,-n);
27 
28  if( TMath::IsNaN(val) )
29  val = 0.0;
30 
31  return val;
32 }
33 
34 double CBFunction_withBG(double *x, double *p)
35 {
36  double norm = p[0];
37  double alpha = p[1]; // tail start
38  double n = p[2]; // tail shape
39  double mu = p[3]; // upsilon mass
40  double sigma = p[4]; // upsilon width
41 
42  double A = pow(n/fabs(alpha),n)*TMath::Exp(-pow(fabs(alpha),2)/2.);
43  double B = n/fabs(alpha) - fabs(alpha);
44  double k = (x[0]-mu)/sigma;
45 
46  double val;
47  if( k > -alpha )
48  val = norm*TMath::Exp(-0.5*pow(k,2));
49  else
50  val = norm*A*pow(B-k,-n);
51 
52  double bgnorm1 = p[5];
53  double bgslope1 = p[6];
54  double bg = exp(bgnorm1+x[0]*bgslope1);
55 
56  val = val + bg;
57 
58  if( TMath::IsNaN(val) )
59  val = 0.0;
60 
61  return val;
62 }
63 //---- ALICE STUFF ------------------------------------------------------------------
64 
65 Double_t CrystallBall2( Double_t x, Double_t mean, Double_t sigma, Double_t alpha1, Double_t n1, Double_t alpha2, Double_t n2 )
66 {
67  Double_t t = (x-mean)/sigma;
68  if( t < -alpha1 )
69  {
70  Double_t a = TMath::Power( n1/alpha1, n1 )*TMath::Exp( -TMath::Power( alpha1, 2 )/2 );
71  Double_t b = n1/alpha1 - alpha1;
72  return a/TMath::Power( b - t, n1 );
73  } else if( t > alpha2 ) {
74 
75  Double_t a = TMath::Power( n2/alpha2, n2 )*TMath::Exp( -TMath::Power( alpha2,2 )/2 );
76  Double_t b = n2/alpha2 - alpha2;
77  return a/TMath::Power( b + t, n2 );
78  } else return TMath::Exp( -TMath::Power( t, 2 )/2 );
79 }
80 
81 Double_t CrystallBall2Integral( Double_t sigma, Double_t alpha1, Double_t n1, Double_t alpha2, Double_t n2 )
82 {
83 // get corresponding integral
84  alpha1 = fabs( alpha1 );
85  alpha2 = fabs( alpha2 );
86  return sigma*(
87  n1/(alpha1*(n1-1))*TMath::Exp( -pow( alpha1, 2 )/2 ) +
88  n2/(alpha2*(n2-1))*TMath::Exp( -pow( alpha2,2 )/2 ) +
89  TMath::Sqrt( TMath::Pi()/2)*TMath::Erfc( -alpha1/TMath::Sqrt(2) ) -
90  TMath::Sqrt( TMath::Pi()/2)*TMath::Erfc( alpha2/TMath::Sqrt(2) ) );
91 }
92 
93 Double_t CrystallBall2( Double_t *x, Double_t *par )
94 {
95  // get normalized Crystal ball
96  Double_t result = CrystallBall2( x[0], par[1], par[2], par[3], par[4], par[5], par[6] );
97  // get integral
98  Double_t integral = CrystallBall2Integral( par[2], par[3], par[4], par[5], par[6] );
99  // return scaled Crystalball so that par[0] corresponds to integral
100  return par[0] * result/integral;
101 }
102 
103 //======================================================================
104 
105 void plotembed_vscent()
106 {
107 SetsPhenixStyle();
108 gStyle->SetOptStat(0);
109 gStyle->SetOptFit(1);
110 
111 double RR = 15.;
112 const int nbins = 10;
113 double bins[nbins+1];
114 double AA = RR*RR/double(nbins);
115 bins[0] = 0.;
116 bins[nbins] = 16.;
117 for(int i=1; i<nbins; i++) { bins[i] = sqrt(i*AA); }
118 for(int i=0; i<=nbins; i++) { cout << "bin edge: " << i << " " << bins[i] << endl; }
119 
120 TF1* fCB = new TF1("fCB",CBFunction,6.,12.,5);
121 TF1* fCBbg = new TF1("fCBbg",CBFunction_withBG,6.,12.,7);
122 
123 float mass,type,pt,eta,rap,pt1,pt2,eta1,eta2,bimp;
124 float chisq1,chisq2,dca2d1,dca2d2,eop3x3_1,eop3x3_2;
125 float nmvtx1,nmvtx2,ntpc1,ntpc2;
126 
127 char hname[99];
128 TH1D* hm[nbins];
129 TH1D* hmss[nbins];
130 TH1D* hmnbg[nbins];
131 TH1D* hbimp[nbins];
132 for(int i=0; i<nbins; i++) {
133  sprintf(hname,"hm_%d",i);
134  hm[i] = new TH1D(hname,hname,120,6.,12.);
135  sprintf(hname,"hmss_%d",i);
136  hmss[i] = new TH1D(hname,hname,120,6.,12.);
137  sprintf(hname,"hmnbg_%d",i);
138  hmnbg[i] = new TH1D(hname,hname,120,6.,12.);
139  sprintf(hname,"hbimp_%d",i);
140  hbimp[i] = new TH1D(hname,hname,1000,0.,20.);
141  hm[i]->Sumw2();
142  hmss[i]->Sumw2();
143  hmnbg[i]->Sumw2();
144 }
145 
146 TH1D* hmass = new TH1D("hmass","",120,6.,12.);
147 TH1D* hmass_ss = new TH1D("hmass_ss","",120,6.,12.);
148 TH1D* hmass_nobg = new TH1D("hmass_nobg","",120,6.,12.);
149 hmass->Sumw2();
150 hmass_ss->Sumw2();
151 hmass_nobg->Sumw2();
152 
153 TLine* l1 = new TLine(7.,0.,11.,0.);
154 l1->SetLineStyle(2);
155 
156 //TNtuple* ntp2;
157 TChain* cntp2 = new TChain("ntp2");
158 TChain* cntp2c = new TChain("ntp2");
159 
160 //-------------------------------------------------------------------
161 
162 // central ---------------------------------------
163 
164 string infname_central;
165 ifstream ifs_central;
166 ifs_central.open("centrallist.txt");
167 if ( ifs_central.fail() ) { cout << "FAIL TO READ INPUT FILE 1" << endl; ifs_central.close(); return; }
168 while(!ifs_central.eof()) {
169  ifs_central >> infname_central;
170  cntp2c->AddFile(infname_central.c_str());
171 }
172 ifs_central.close();
173 
174 cout << "number of CENTRAL entries = " << cntp2c->GetEntries() << endl;
175 
176 cntp2c->SetBranchAddress("type",&type);
177 cntp2c->SetBranchAddress("mass",&mass);
178 cntp2c->SetBranchAddress("pt",&pt);
179 cntp2c->SetBranchAddress("eta",&eta);
180 cntp2c->SetBranchAddress("rap",&rap);
181 cntp2c->SetBranchAddress("pt1",&pt1);
182 cntp2c->SetBranchAddress("pt2",&pt2);
183 cntp2c->SetBranchAddress("eta1",&eta1);
184 cntp2c->SetBranchAddress("eta2",&eta2);
185 cntp2c->SetBranchAddress("chisq1",&chisq1);
186 cntp2c->SetBranchAddress("dca2d1",&dca2d1);
187 cntp2c->SetBranchAddress("eop3x3_1",&eop3x3_1);
188 cntp2c->SetBranchAddress("chisq2",&chisq2);
189 cntp2c->SetBranchAddress("dca2d2",&dca2d2);
190 cntp2c->SetBranchAddress("eop3x3_2",&eop3x3_2);
191 cntp2c->SetBranchAddress("nmvtx1",&nmvtx1);
192 cntp2c->SetBranchAddress("nmvtx2",&nmvtx2);
193 cntp2c->SetBranchAddress("ntpc1",&ntpc1);
194 cntp2c->SetBranchAddress("ntpc2",&ntpc2);
195 
196  for(int j=0; j<cntp2c->GetEntries(); j++) {
197  cntp2c->GetEvent(j);
198  if(chisq1>3 || chisq2>3) continue;
199  if(nmvtx1<2 || nmvtx2<2) continue;
200  if(ntpc1<26 || ntpc2<26) continue;
201  if(eop3x3_1<0.7 || eop3x3_2<0.7) continue;
202  //if(type==1) { hmass->Fill(mass); hmass_nobg->Fill(mass); }
203  //if(type==2 || type==3) { hmass_ss->Fill(mass); }
204  if(type==1) { hm[0]->Fill(mass); hmnbg[0]->Fill(mass); }
205  if(type==2 || type==3) { hmss[0]->Fill(mass); }
206  }
207 // end central -----------------------------------
208 
209 // minbias ---------------------------------------
210 string infname;
211 ifstream ifs;
212 ifs.open("mblist.txt");
213 if ( ifs.fail() ) { cout << "FAIL TO READ INPUT FILE 1" << endl; ifs.close(); return; }
214 while(!ifs.eof()) {
215  ifs >> infname;
216  cntp2->AddFile(infname.c_str());
217 }
218 ifs.close();
219 
220 cout << "number of MB entries = " << cntp2->GetEntries() << endl;
221 
222 cntp2->SetBranchAddress("b",&bimp);
223 cntp2->SetBranchAddress("type",&type);
224 cntp2->SetBranchAddress("mass",&mass);
225 cntp2->SetBranchAddress("pt",&pt);
226 cntp2->SetBranchAddress("eta",&eta);
227 cntp2->SetBranchAddress("rap",&rap);
228 cntp2->SetBranchAddress("pt1",&pt1);
229 cntp2->SetBranchAddress("pt2",&pt2);
230 cntp2->SetBranchAddress("eta1",&eta1);
231 cntp2->SetBranchAddress("eta2",&eta2);
232 cntp2->SetBranchAddress("chisq1",&chisq1);
233 cntp2->SetBranchAddress("dca2d1",&dca2d1);
234 cntp2->SetBranchAddress("eop3x3_1",&eop3x3_1);
235 cntp2->SetBranchAddress("chisq2",&chisq2);
236 cntp2->SetBranchAddress("dca2d2",&dca2d2);
237 cntp2->SetBranchAddress("eop3x3_2",&eop3x3_2);
238 cntp2->SetBranchAddress("nmvtx1",&nmvtx1);
239 cntp2->SetBranchAddress("nmvtx2",&nmvtx2);
240 cntp2->SetBranchAddress("ntpc1",&ntpc1);
241 cntp2->SetBranchAddress("ntpc2",&ntpc2);
242 
243  for(int j=0; j<cntp2->GetEntries(); j++) {
244  cntp2->GetEvent(j);
245  if(chisq1>3 || chisq2>3) continue;
246  if(nmvtx1<2 || nmvtx2<2) continue;
247  if(ntpc1<26 || ntpc2<26) continue;
248  if(eop3x3_1<0.7 || eop3x3_2<0.7) continue;
249  if(type==1) { hmass->Fill(mass); hmass_nobg->Fill(mass); }
250  if(type==2 || type==3) { hmass_ss->Fill(mass); }
251 
252  for(int j=0; j<nbins; j++) {
253  if(bimp>bins[j] && bimp<bins[j+1]) {
254  if(type==1) { hm[j]->Fill(mass); hmnbg[j]->Fill(mass); hbimp[j]->Fill(bimp); }
255  if(type==2 || type==3) { hmss[j]->Fill(mass); }
256  }
257  }
258 
259  }
260 // end minbias -------------------------------------------------------
261 
262 for(int j=0; j<nbins; j++) { hmnbg[j]->Add(hmss[j],-1.0); }
263 hmass_nobg->Add(hmass_ss,-1.0);
264 
265 double massres[nbins];
266 double massreserr[nbins];
267 double x[nbins];
268 for(int i=0; i<nbins; i++) { x[i] = hbimp[i]->GetMean(); cout << "mean b = " << x[i] << endl; }
269 
270 //----------------------------------------------------------------
271 
272 // ALICE Double Crystal Ball function
273 TF1 *f2 = new TF1("f2",CrystallBall2,7,11,7);
274  f2->SetParameter(0, 2000. );
275  f2->SetParameter(1, 9.46 );
276  f2->SetParameter(2, 0.13 );
277  f2->SetParameter(3, 1);
278  f2->SetParameter(4, 3);
279  f2->SetParameter(5, 1 );
280  f2->SetParameter(6, 5 );
281  f2->SetParNames("normalization", "mean", "sigma","alpha1","n1","alpha2","n2");
282  f2->SetParLimits(1, 9.40, 9.55);
283  f2->SetParLimits(2, 0.06, 0.20);
284  f2->SetParLimits(3, 0.120, 10);
285  f2->SetParLimits(4, 1.01, 10);
286  f2->SetParLimits(5, 0.1, 10);
287  f2->SetParLimits(6, 1.01, 10);
288 
289 
290 TCanvas* c3 = new TCanvas("c3","inv. mass allpT",20,20,600,500);
291 
292 hmass->Draw();
293 hmass_ss->SetLineColor(kRed);
294 hmass_ss->Scale(1.2);
295 hmass_ss->Draw("same");
296 
297 TCanvas* c33 = new TCanvas("c33","inv. mass nobg",20,20,600,500);
298 //c33->Divide(3,2);
299 /*
300  fCB->SetParameter(0,1000.);
301  fCB->SetParameter(1,1.0);
302  fCB->SetParameter(2,1.0);
303  fCB->SetParLimits(2,0.1,99.);
304  fCB->SetParameter(3,9.0);
305  fCB->SetParameter(4,0.100);
306 */
307 
308 // c33->cd(1);
309  hmass_nobg->SetAxisRange(7.,11.);
310  hmass_nobg->Fit("f2","rlq","",7.5,10.0);
311 
312  for(int j=0; j<nbins; j++) {
313 // c33->cd(j+2);
314  fCB->SetParameter(0,hmnbg[j]->GetMaximum());
315  hmnbg[j]->SetTitle(";Invariant mass, GeV ; Counts");
316  hmnbg[j]->Fit("f2","rlq","",7.5,10.0);
317  hmnbg[j]->SetAxisRange(7.,11.);
318 // l1->Draw();
319  cout << "mass resolution = " << j << " " << f2->GetParameter(2) << " +- " << f2->GetParError(2) << endl;
320  massres[j] = 1000.*f2->GetParameter(2);
321  massreserr[j] = 1000.*f2->GetParError(2);
322  }
323  hmass_nobg->SetTitle(";Invariant mass, GeV ; Counts");
324  hmass_nobg->Draw();
325  TLatex* lat1 = new TLatex(7.2,hmass_nobg->GetMaximum()*0.4,"Minimum Bias"); lat1->Draw();
326  l1->Draw();
327 
328 
329 // double fitndf = hmass_nobg->GetFunction("fCB")->GetNDF();
330 // double fitchi2 = hmass_nobg->GetFunction("fCB")->GetChisquare();
331 // cout << "chi2 = " << fitchi2/fitndf << endl;
332 
333 double scale = 1.4;
334 TF1* fws = new TF1("fws","[0]/(1+exp((x-[1])/[2]))",0.,20.);
335 fws->SetParameter(0,8.68398e+01*scale);
336 fws->SetParameter(1,1.51845e+01);
337 fws->SetParameter(2,6.52512e-01);
338 TF1* fws_sum = new TF1("fws_sum","(x*[0])/(1+exp((x-[1])/[2]))",0.,20.);
339 fws_sum->SetParameter(0,8.68398e+01*scale);
340 fws_sum->SetParameter(1,1.51845e+01);
341 fws_sum->SetParameter(2,6.52512e-01);
342 double avg[nbins];
343 for(int j=0; j<nbins; j++) {
344  avg[j] = fws_sum->Integral(bins[j],bins[j+1]);
345 // cout << "expected = " << avg[j] << endl;
346 }
347 avg[0] = avg[0]*0.9;
348 
349 double sumstart = 8.8;
350 double sumstop = 9.8;
351 int fbin = hmass_nobg->FindBin(sumstart + 0.001);
352 int lbin = hmass_nobg->FindBin(sumstop - 0.001);
353 cout << "counting bins from " << fbin << " to " << lbin << endl;
354 cout << " in mass range from " << hmass_nobg->GetBinLowEdge(fbin) << " to " << hmass_nobg->GetBinLowEdge(lbin)+hmass_nobg->GetBinWidth(lbin) << endl;
355 
356  double sum =0.;
357  double sumsum[nbins]; for(int i=0; i<nbins; i++) { sumsum[i]=0.; }
358  double y[nbins],ey[nbins];
359  for(int k=fbin; k<=lbin; k++) {
360  sum += hmass_nobg->GetBinContent(k);
361  for(int i=0; i<nbins; i++) {
362  sumsum[i] += hmnbg[i]->GetBinContent(k);
363  y[i] = sumsum[i]/avg[i];
364  ey[i] = y[i]*0.05;
365  }
366  }
367  cout << "Number of Upsilons = " << sum << endl;
368  for(int i=0; i<nbins; i++) { cout << sumsum[i] << endl; }
369 
370  cout << "Number of Upsilons from fit = " << fCB->Integral(sumstart,sumstop)/hmass_nobg->GetBinWidth(lbin) << " " << fCB->Integral(5.0,11.0)/hmass_nobg->GetBinWidth(lbin) << endl;
371 
372 TCanvas* c4 = new TCanvas("c4","mass resolution vs Bimp",20,20,600,500);
373 TH2D* hh = new TH2D("hh","",10,0.,18.,10,60.,140.);
374 hh->SetTitle(";Impact parameter, fm ; Mass resolution, MeV");
375 hh->Draw();
376 
377 TGraphErrors* gr1 = new TGraphErrors(nbins,x,massres,0,massreserr);
378 gr1->SetMarkerStyle(20);
379 gr1->SetMarkerColor(kBlack);
380 gr1->SetMarkerSize(1.5);
381 gr1->SetLineColor(kBlack);
382 gr1->SetLineWidth(2);
383 gr1->Draw("p");
384 
385 /*
386 TCanvas* c44 = new TCanvas("c44","mass resolution vs Bimp",120,120,600,500);
387 TH2D* hhh = new TH2D("hhh","",10,0.,18.,10,0.,1.5);
388 hhh->Draw();
389 
390 TGraphErrors* gr2 = new TGraphErrors(nbins,x,y,0,ey);
391 gr2->SetMarkerStyle(20);
392 gr2->SetMarkerColor(kBlue);
393 gr2->SetMarkerSize(1.5);
394 gr2->SetLineColor(kBlue);
395 gr2->Draw("p");
396 */
397 
398 TCanvas* c5 = new TCanvas("c5","",60,60,600,500);
399 hmnbg[0]->Draw();
400 TLatex* lat2 = new TLatex(7.2,hmnbg[0]->GetMaximum()*0.4,"10% central Au+Au"); lat2->Draw();
401 l1->Draw();
402 
403 }
404