FairRoot/PandaRoot
analyse_phi2.C
Go to the documentation of this file.
1 // --------------------------------------------------------------------
2 
3 void confgraph(TGraph *g, TString tit, int col=1, int marker=20)
4 {
5  g->GetHistogram()->SetTitle(tit);
6  g->GetHistogram()->SetMinimum(0);
7  g->SetLineColor(col);
8  g->SetMarkerColor(col);
9  g->SetMarkerStyle(marker);
10  g->SetMarkerSize(1);
11  g->SetLineWidth(2);
12 }
13 
14 // --------------------------------------------------------------------
15 
16 void confhist(TH1* h)
17 {
18  h->GetXaxis()->SetLabelSize(0.05);
19  h->GetXaxis()->SetTitleSize(0.05);
20  h->GetYaxis()->SetLabelSize(0.05);
21  h->GetYaxis()->SetTitleSize(0.05);
22  h->GetYaxis()->SetTitleOffset(1.4);
23 
24  h->GetXaxis()->SetNdivisions(505);
25  h->SetMinimum(0);
26  h->SetStats(0);
27 }
28 
29 // --------------------------------------------------------------------
30 
31 TH1F* createHistoGraph(TGraph *g, TString tit="", double xmin=0, double xmax=0, double ymaxf=1.05)
32 {
33  static int cnt=0;
34  if (tit=="") tit=g->GetTitle();
35 
36  double ymax = TMath::MaxElement(g->GetN(), g->GetY());
37  double dymax = g->GetErrorY(TMath::LocMax(g->GetN(), g->GetEY()));
38 
39  if (xmin==xmax)
40  {
41  xmin = TMath::MinElement(g->GetN(), g->GetX());
42  xmax = TMath::MaxElement(g->GetN(), g->GetX());
43  }
44  TH1F *h=new TH1F(Form("h%03d",cnt++),tit,g->GetN(),xmin-0.5*(xmax-xmin)/g->GetN(), xmax+0.5*(xmax-xmin)/g->GetN());
45 
46  confhist(h);
47  //h->GetXaxis()->SetLabelSize(0.05);
48  //h->GetXaxis()->SetTitleSize(0.05);
49  //h->GetYaxis()->SetLabelSize(0.05);
50  //h->GetYaxis()->SetTitleSize(0.05);
51  //h->GetYaxis()->SetTitleOffset(1.4);
52 
53  //h->GetXaxis()->SetNdivisions(505);
54  //h->SetMinimum(0);
55  //h->SetStats(0);
56 
57  h->SetMaximum((ymax+dymax)*ymaxf);
58 
59  return h;
60 }
61 // --------------------------------------------------------------------
62 // --------------------------------------------------------------------
63 
64 
65 void analyse_phi2(double t=2, double sigsmin=500, double sigsmax=500, int Nrep=1000, int steps = 25, double emin = 2.25, double emax = 2.35, double m=2.3, double G=0.015)
66 {
67  //gStyle->SetStatX(0.5);
68  //gStyle->SetStatY(0.40);
69  gStyle->SetStatW(0.20);
70  gStyle->SetStatH(0.12);
71 
72  double sigs=(sigsmax+sigsmin)*0.5;
73 
74  // C. Evangelista et al., Phys. Rev. D57, 5370 (1998).
75  double x[25] = {2.219, 2.221, 2.221, 2.222, 2.224, 2.226, 2.226, 2.228, 2.229, 2.231, 2.231, 2.233, 2.235, 2.236, 2.242, 2.247, 2.254, 2.255, 2.272, 2.289, 2.307, 2.333, 2.360, 2.404, 2.430};
76  double dx[25] = {0.0};
77 
78  double y[25] = {4.15, 4.16, 3.87, 3.94, 4.40, 3.90, 3.68, 3.00, 3.20, 3.73, 2.97, 4.20, 3.52, 2.48, 3.38, 2.68, 2.63, 3.01, 2.25, 2.15, 1.97, 2.04, 2.04, 1.61, 1.93 };
79  double dy[25] = {0.34, 0.36, 0.54, 0.35, 0.36, 0.34, 0.34, 0.30, 0.32, 0.35, 0.33, 0.35, 0.34, 0.30, 0.45, 0.39, 0.35, 0.16, 0.17, 0.21, 0.24, 0.14, 0.12, 0.18, 0.11 };
80 
81 
82  TGraphErrors *gphi = new TGraphErrors(25,x,y,dx,dy);
83 
84  confgraph(gphi,"#phi#phi cross section [JetSet];E_{cm} [GeV];#sigma [#mub]");
85  TH1F *hphi = createHistoGraph(gphi,"#phi#phi cross section [JetSet];E_{cm} [GeV];#sigma [#mub]",2.21,2.435);
86 
87  //TF1 *f1=new TF1("f1","expo",2.21,2.435);
88  TF1 *f1=new TF1("f1","0.5*[0]*(1.0-TMath::Erf((x-[1])/[2]))+[3]",2.21,2.435);
89  f1->SetParameters(3,2.23,0.01,1);
90 
91  // inelasic cross section in mb
92  TF1 *finel=new TF1("finel","[0]+[1]/x^[2]-[3]/x^[4]");
93  finel->SetParameters(32,82,0.67,45,0.85);
94  finel->SetRange(1.,100.);
95 
96  TCanvas *c1=new TCanvas("c1","c1",20,20,1000,1000);
97  c1->Divide(2,2,0.0001,0.0001);
98 
99  c1->cd(1);
100  hphi->Draw();
101  gphi->Fit("f1","R");
102  gphi->Draw("Psame");
103 
104  TLine l;
105  l.SetLineColor(4);
106  l.SetLineStyle(9);
107  l.SetLineWidth(2);
108  l.DrawLine(emin,0,emin,hphi->GetMaximum()*0.5);
109  l.DrawLine(emax,0,emax,hphi->GetMaximum()*0.5);
110 
111 
112  double effs = 22463./100000.; // signal efficiency: cut = abs(summ4-2.15)<0.15&&abs(f4cxd0m-1.02)<0.012&&abs(f4cxd1m-1.02)<0.012&chi24c<50
113  //double sigs = 100.; // signal cross section [nb]
114  double fBR = 0.489*0.489; // BR(phi->K+ K-)^2
115 
116  double effg = 77./840.5e6; // generic background efficiency
117 
118  double sigphi = f1->Eval(2.3)*1000; // phi phi cross section [nb]
119  double sigg = 60e6; // generic cross section [nb]
120 
121  //double t = 0.5; // [d]
122  double Lint = 788; // [1./(nb*d)]
123 
124  double S = sigs*effs*fBR*t*Lint; // expected number of signals
125  double Bphi = sigphi*effs*fBR*t*Lint; // expected number of non-res phi phi
126  double Bg = sigg*effg*t*Lint; // expected number of
127 
128  cout <<"S="<<S<<" B_phi(sig ="<<sigphi<<" nb)="<<Bphi<<" B_g="<<Bg<<endl;
129 
130 
131  // signal cross section (BW)
132  TF1 *fbw = new TF1("fbw","[0]*[3]*TMath::BreitWigner(x,[1],[2])");
133  fbw->SetParNames("#sigma_{S}","m","#Gamma","f_{scale}");
134 
135  // background function
136  TF1 *fbg=new TF1("fbg","pol3",2.21,2.435);
137 
138  // combined fit function
139  TF1 *ffit = new TF1("ffit","[0]*TMath::BreitWigner(x,[1],[2]) + pol3(3)");
140  ffit->SetParNames("A","m","#Gamma","a_{0}","a_{1}","a_{2}","a_{3}");
141 
142 
143 
144 
145  double stsize = (emax-emin)/(steps-1);
146 
147  TGraphErrors *gscan = new TGraphErrors(steps);
148  TGraphErrors *gscansig = new TGraphErrors(steps);
149  TGraphErrors *gscanbg = new TGraphErrors(steps);
150  confgraph(gscan,Form("#phi#phi scan (t=%.2fd);E_{cm} [GeV];counts",t));
151  confgraph(gscanbg,Form("#phi#phi scan (t=%.2fd);E_{cm} [GeV];counts",t),16);
152  confgraph(gscansig,Form("#phi#phi scan (t=%.2fd);E_{cm} [GeV];counts",t));
153 
154  double st_t = t/steps;
155 
156  TH1F *hsigs = new TH1F("hsigs",Form("#sigma_{S}, %d fits",Nrep),100,0,sigs*3);
157  confhist(hsigs);
158  hsigs->SetLineWidth(2);
159 
160  TH1F *hscan=0, *hscansig=0;
161 
162  for (int j=0;j<Nrep;++j)
163  {
164  fbw->SetParameters(1,m,G,1);
165  fbw->FixParameter(3,1);
166 
167  for (int i=0;i<steps;++i)
168  {
169  double ecm = emin + i*stsize;
170 
171  double st_sigs = fbw->Eval(ecm)/fbw->Eval(m)*sigs;
172  double st_sigphi = f1->Eval(ecm)*1000; // mub -> nb
173  double st_sigg = finel->Eval(ecm)*1e6; // mub -> nb
174 
175  double yieldS = gRandom->Poisson(st_sigs * effs * fBR * st_t * Lint);
176  double yieldBphi = gRandom->Poisson(st_sigphi * effs * fBR * st_t * Lint);
177  double yieldBg = gRandom->Poisson(st_sigg * effg * st_t * Lint);
178 
179  //printf("(%2d) Ecm = %5.3f GeV : sig_S = %.2f nb (S = %6.0f) sig_phi = %.2f nb (B = %6.0f) sig_gen = %.2f mub (Bg = %5.0f)\n", i, ecm, st_sigs, yieldS, st_sigphi, yieldBphi, st_sigg/1000, yieldBg);
180 
181  double yield = (yieldS + yieldBg + yieldBphi);
182 
183  gscan->SetPoint(i,ecm, yield);
184  gscan->SetPointError(i, 0, sqrt(yield));
185 
186  gscanbg->SetPoint(i,ecm, gRandom->Poisson(yieldBg));
187  gscanbg->SetPointError(i, 0, sqrt(yieldBg));
188  }
189 
190  ffit->SetParameters(sigs,m,G,1,1,1);
191  ffit->FixParameter(1,m);
192  ffit->FixParameter(2,G);
193 
194  gscan->Fit("ffit","0q");
195  ffit->SetParLimits(1,m*0.7,m*1.3);
196  ffit->SetParLimits(2,G*0.5,G*1.5);
197  gscan->Fit("ffit","qm");
198 
199  c1->cd(2);
200  if (hscan==0) hscan = createHistoGraph(gscan,Form("#phi#phi scan (t = %.2fd, #sigma = %.0fnb);E_{cm} [GeV];counts",t,sigs),0,0,1.3);
201  hscan->Draw();
202  gscan->Draw("Psame");
203  gscanbg->Draw("Psame");
204 
205  fbg->SetParameters(ffit->GetParameter(3),ffit->GetParameter(4),ffit->GetParameter(5),ffit->GetParameter(6));
206 
207 
208  for (int i=0;i<steps;++i)
209  {
210  double x,y,dx,dy;
211  gscan->GetPoint(i,x,y);
212  gscansig->SetPoint(i,x,y-fbg->Eval(x));
213  gscansig->SetPointError(i,gscan->GetErrorX(i),gscan->GetErrorY(i));
214 
215  //cout <<x<<" "<<fbg->Eval(x)<<endl;
216  }
217 
218  fbw->FixParameter(3,1);
219  fbw->SetParLimits(1,m*0.7,m*1.3);
220  fbw->SetParLimits(2,G*0.6,G*1.5);
221  gscansig->Fit("fbw","q");
222  fbw->SetParameter(0,1);
223  //cout<<fbw->Eval(fbw->GetParameter(1))<<" "<<fbw->GetParameter(0)<<" "<<fbw->GetParameter(1)<<" "<<fbw->GetParameter(2)<<" "<<fbw->GetParameter(3)<<endl;
224  fbw->FixParameter(3, 1./fbw->Eval(fbw->GetParameter(1))*(effs*fBR*st_t*Lint) );
225  fbw->SetParameter(0,TMath::MaxElement(gscansig->GetN(),gscansig->GetY())/(effs*fBR*st_t*Lint));
226  gscansig->Fit("fbw","q");
227 
228  c1->cd(3);
229  if (hscansig==0) hscansig = createHistoGraph(gscansig,Form("#phi#phi scan (t = %.2fd, #sigma = %.0fnb);E_{cm} [GeV];counts",t,sigs),0,0,1.3);
230  hscansig->SetMinimum(-100);
231  hscansig->Draw();
232  gscansig->Draw("Psame");
233 
234  //fbw->Print();
235  double ss = fbw->GetParameter(0);
236  hsigs->Fill(ss);
237 
238 
239  //cout <<"sig_in = "<<sigs<<" sig_out = "<< ss<<" "<<TMath::MaxElement(gscansig->GetN(),gscansig->GetY())<<" "<<fbw->GetParameter(0)<<" "<<fbw->GetParameter(1)<<" "<<fbw->GetParameter(2)<<" "<<fbw->GetParameter(3)<<endl;
240 
241  c1->cd(4);
242  hsigs->Draw();
243 
244  if (j%100==0) c1->Update();
245 
246  //if (ss>880) return;
247  //sleep(1);
248  }
249 
250  //c1->cd(2);
251  //TH1F *hscan = createHistoGraph(gscan,Form("#phi#phi scan (t = %.2fd, #sigma = %.0fnb);E_{cm} [GeV];counts",t,sigs),0,0,1.3);
252  //hscan->Draw();
253  //gscan->Draw("Psame");
254  //gscanbg->Draw("Psame");
255 
256  //c1->cd(3);
257  //TH1F *hscansig = createHistoGraph(gscansig,Form("#phi#phi scan (t = %.2fd, #sigma = %.0fnb);E_{cm} [GeV];counts",t,sigs),0,0,1.3);
258  //hscansig->SetMinimum(-100);
259  //hscansig->Draw();
260  //gscansig->Draw("Psame");
261 
262  c1->cd(4);
263  hsigs->Draw("e");
264  hsigs->SetStats(1);
265  hsigs->Fit("gaus");
266 
267  //fbw->FixParameter(3,1);
268  //fbw->SetParLimits(1,m*0.7,m*1.3);
269  //fbw->SetParLimits(2,G*0.5,G*1.5);
270  //gscansig->Fit("fbw");
271  //fbw->SetParameter(0,1);
272  //fbw->FixParameter(3, 1./fbw->GetMaximum(emin,emax) );
273  //fbw->SetParameter(0,TMath::MaxElement(gscansig->GetN(),gscansig->GetY()));
274  //gscansig->Fit("fbw");
275 
276 
277  c1->Update();
278  TPaveStats *st = (TPaveStats *)gscan->GetListOfFunctions()->FindObject("stats");
279  if (st)
280  {
281  cout <<"good"<<endl;
282  st->SetX1NDC(0.5); //new x start position
283  st->SetX2NDC(0.5); //new x end position
284  c1->Modified();
285  c1->Update();
286  }
287  else cout <<"***"<<endl;
288 
289 }
290 
dy
double dy
Definition: runGeaneTestSimBox.cxx:45
i
Int_t i
Definition: run_full.C:25
m
__m128 m
Definition: P4_F32vec4.h:28
f1
TF1 * f1
Definition: reco_analys2.C:50
sqrt
friend F32vec4 sqrt(const F32vec4 &a)
Definition: P4_F32vec4.h:29
xmax
Double_t xmax
Definition: convolutionAnalysis.C:11
col
int col
Definition: anaLmdDigi.C:67
g
TFile * g
Definition: GammaSpectraAnalysis.C:22
confhist
void confhist(TH1 *h)
Definition: analyse_phi2.C:16
createHistoGraph
TH1F * createHistoGraph(TGraph *g, TString tit="", double xmin=0, double xmax=0)
Definition: analyse_etac1.C:60
confgraph
void confgraph(TGraph *g, TString tit, int col=1, int marker=20)
Definition: analyse_etac1.C:47
analyse_phi2
void analyse_phi2(double t=2, double sigsmin=500, double sigsmax=500, int Nrep=1000, int steps=25, double emin=2.25, double emax=2.35, double m=2.3, double G=0.015)
Definition: analyse_phi2.C:65
c1
c1
Definition: plot_dirc.C:35
dx
double dx
Definition: runGeaneTestSimBox.cxx:44
effs
std::map< int, double > effs
Definition: simubg.C:29
x
Double_t x
Definition: createRootFscGeometryFile.C:159
cnt
Int_t cnt
Definition: hist-t7.C:106
t
TTree * t
Definition: bump_analys.C:13
xmin
Double_t xmin
Definition: convolutionAnalysis.C:10
y
Double_t y
Definition: createRootFscGeometryFile.C:159
Generated on Wed Apr 3 2019 08:43:48 for FairRoot/PandaRoot by   doxygen 1.8.5