FairRoot/PandaRoot
analyse_etac1_slc.C
Go to the documentation of this file.
1 // --------------------------------------------------------------------
2 // counts different events for a cut
3 long int cntEvt(TTree *t1, TString cut, TString var="ev")
4 {
5  float fbev;
6  int ibev;
7 
8  TBranch *br = t1->GetBranch(var);
9  if (br==0) return 0;
10 
11  TString tit = br->GetTitle(); // should end with '/I' or '/F'
12 
13  t1->SetEventList(0);
14  t1->Draw(">>el",cut);
15  TEventList *el=(TEventList*)gROOT->FindObject("el");
16 
17  t1->SetBranchStatus("*",0);
18  t1->SetBranchStatus(var,1);
19 
20  // int or float branch
21  bool isint=tit.EndsWith("/I");
22 
23  if (isint)
24  t1->SetBranchAddress(var,&ibev);
25  else
26  t1->SetBranchAddress(var,&fbev);
27 
28  // last and current event number; counter
29  int lev = -1, cev=-1, cnt=0;
30  for (int i=0;i<el->GetN();++i)
31  {
32  t1->GetEvent(el->GetEntry(i));
33 
34  cev = isint ? ibev : fbev;
35 
36  if (lev!=cev) cnt++;
37  lev = cev;
38  }
39 
40  t1->SetBranchStatus("*",1);
41 
42  return cnt;
43 }
44 
45 // --------------------------------------------------------------------
46 
47 void confgraph(TGraph *g, TString tit, int col=1, int marker=20)
48 {
49  g->GetHistogram()->SetTitle(tit);
50  g->GetHistogram()->SetMinimum(0);
51  g->SetLineColor(col);
52  g->SetMarkerColor(col);
53  g->SetMarkerStyle(marker);
54  g->SetMarkerSize(1);
55  g->SetLineWidth(2);
56 }
57 
58 // --------------------------------------------------------------------
59 
60 TH1F* createHistoGraph(TGraph *g, TString tit="", double xmin=0, double xmax=0)
61 {
62  static int cnt=0;
63  if (tit=="") tit=g->GetTitle();
64 
65  double ymax = TMath::MaxElement(g->GetN(), g->GetY());
66  double dymax = g->GetErrorY(TMath::LocMax(g->GetN(), g->GetEY()));
67 
68  if (xmin>=xmax)
69  {
70  xmin = TMath::MinElement(g->GetN(), g->GetX());
71  xmax = TMath::MaxElement(g->GetN(), g->GetX());
72  }
73  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());
74 
75  h->GetXaxis()->SetLabelSize(0.05);
76  h->GetXaxis()->SetTitleSize(0.05);
77  h->GetYaxis()->SetLabelSize(0.05);
78  h->GetYaxis()->SetTitleSize(0.05);
79  h->GetYaxis()->SetTitleOffset(1.4);
80 
81  h->SetMaximum((ymax+dymax)*1.05);
82  h->SetMinimum(0);
83 
84  h->SetStats(0);
85 
86  return h;
87 }
88 
89 // --------------------------------------------------------------------
90 
91 void analyse_etac1_slc(TString fname="ntp5_slc_etac1_A.root", TString cut="chi24c<50", int sigmode=300, TString tit="")
92 {
93  //---------------------------------------------------------------------------------------------------------------------------------
94  //Processing analyse_etac1_slc.C("ntp5_slc_etac1_A.root","chi24c<50",300, "Setup A, #bar{p}p#rightarrow#tilde{#eta_{c1}}#eta")...
95  //---------------------------------------------------------------------------------------------------------------------------------
96  //i | S0 | B0 | S | B | fS | fB | S* | B* |
97  //----+----------+------------+----------+--------+----------+-----------+---------+---------+
98  //0 | 100000 | 2.691 G | 17743 | 17 | 0.0098 | 95.2150 | 173 | 1619 |
99  //8 | 100000 | 2.690 G | 142 | 4 | 0.0098 | 95.2268 | 1 | 381 |
100  //---------------------------------------------------------------------------------------------------------------------------------
101  //Processing analyse_etac1_slc.C("ntp5_slc_etac1_B.root","chi24c<50",320, "Setup B, #bar{p}p#rightarrow#tilde{#eta_{c1}}#eta")...
102  //----+----------+------------+----------+--------+----------+-----------+---------+---------+
103  //0 | 100000 | 2.675 G | 3503 | 11 | 0.0098 | 95.7917 | 34 | 1054 |
104  //8 | 100000 | 2.819 G | 2 | 1 | 0.0098 | 90.8956 | 0 | 91 |
105  //######################################################################
106  //---------------------------------------------------------------------------------------------------------------------------------
107  //Processing analyse_etac1.C("ntp5_etac1_A.root","chi24c<50",300, "Setup A, #bar{p}p#rightarrow#tilde{#eta_{c1}}#eta")...
108  //---------------------------------------------------------------------------------------------------------------------------------
109  //i | S0 | B0 | S | B | fS | fB | S* | B* |
110  //----+----------+------------+----------+--------+----------+-----------+---------+---------+
111  //0 | 100000 | 2.866 G | 17744 | 12 | 0.0098 | 89.3904 | 173 | 1073 |
112  //7 | 100000 | 2.706 G | 131 | 1 | 0.0098 | 94.6622 | 1 | 95 |
113  //---------------------------------------------------------------------------------------------------------------------------------
114  //Processing analyse_etac1.C("ntp5_etac1_B.root","chi24c<50",320, "Setup B, #bar{p}p#rightarrow#tilde{#eta_{c1}}#eta")...
115  //----+----------+------------+----------+--------+----------+-----------+---------+---------+
116  //0 | 100000 | 2.931 G | 3617 | 14 | 0.0098 | 87.4188 | 35 | 1224 |
117  //7 | 100000 | 2.739 G | 5 | 1 | 0.0098 | 93.5508 | 0 | 94 |
118 
119  //data/DPMetac1_A_SLC1 ( 168/ 200 files) = 2690752186
120  //data/DPMetac1_A_SLC2 ( 182/ 200 files) = 2914674237
121  //data/DPMetac1_A_SLC3 ( 174/ 200 files) = 2786406995
122  //data/DPMetac1_A_SLC4 ( 171/ 200 files) = 2738298947
123  //data/DPMetac1_A_SLC5 ( 176/ 200 files) = 2818739807
124  //data/DPMetac1_A_SLC6 ( 179/ 200 files) = 2866193988
125  //data/DPMetac1_A_SLC7 ( 169/ 200 files) = 2706488227
126  //data/DPMetac1_A_SLC8 ( 175/ 200 files) = 2801957793
127  //data/DPMetac1_A_SLC9 ( 168/ 200 files) = 2690419363
128  //--------------------
129  //data/DPMetac1_B_SLC1 ( 167/ 200 files) = 2674554579
130  //data/DPMetac1_B_SLC2 ( 172/ 200 files) = 2754101731
131  //data/DPMetac1_B_SLC3 ( 172/ 200 files) = 2754649841
132  //data/DPMetac1_B_SLC4 ( 171/ 200 files) = 2738459876
133  //data/DPMetac1_B_SLC5 ( 177/ 200 files) = 2834829596
134  //data/DPMetac1_B_SLC6 ( 173/ 200 files) = 2770416159
135  //data/DPMetac1_B_SLC7 ( 175/ 200 files) = 2802568964
136  //data/DPMetac1_B_SLC8 ( 174/ 200 files) = 2786802048
137  //data/DPMetac1_B_SLC9 ( 176/ 200 files) = 2818616809
138 
139  std::map<long int,long int> evcnts =
140  { {1300, 2690752186}, {1301, 2914674237}, {1302, 2786406995}, {1303, 2738298947}, {1304, 2818739807}, {1305, 2866193988}, {1306, 2706488227}, {1307, 2801957793}, {1308, 2690419363}, // Setup A, pbp -> J/psi (-> e+ e-) pi+ pi-
141  {1320, 2674554579}, {1321, 2754101731}, {1322, 2754649841}, {1323, 2738459876}, {1324,2834829596 }, {1325, 2770416159}, {1326, 2802568964}, {1327, 2786802048}, {1328, 2818616809}
142  };
143 
144 
145  // we merge the values at 0% and 100% from both studies by hand...
146  // **** from theta study
147  //data/DPMetac1_A_EMC1 ( 179/ 199 files) = 2866079091
148  //data/DPMetac1_A_EMC8 ( 169/ 199 files) = 2706466096
149  //--------------------
150  //data/DPMetac1_B_EMC1 ( 183/ 198 files) = 2930718917
151  //data/DPMetac1_B_EMC8 ( 171/ 198 files) = 2738620309
152 
153  evcnts[1300] += 2866079091;
154  evcnts[1320] += 2706466096;
155  evcnts[1308] += 2930718917;
156  evcnts[1328] += 2738620309;
157 
158  std::map<long int,long int> mrgB =
159  { {1300, 17+12}, {1308, 4+1},
160  {1320, 11+14}, {1328, 1+1}
161  };
162 
163  //for ( auto x:evcnts) cout <<"mode "<<x.first<<" : "<<x.second<<endl;
164  //cout <<endl;
165 
166  int bkgmode = sigmode+1000;
167  //if (sigmode%100>19) bkgmode+=1;
168 
169  TFile *f = new TFile(fname);
170  TTree *t = (TTree*)f->Get("ntp5");
171 
172  TFile fana("ana_etac_slc.root","UPDATE");
173 
174  int Nred = 9;
175 
176  TGraphErrors *g[4];
177  g[0] = new TGraphErrors(Nred);
178  g[1] = new TGraphErrors(Nred);
179  g[2] = new TGraphErrors(Nred);
180  g[3] = new TGraphErrors(Nred);
181 
182  confgraph(g[0], "signal to noise");
183  confgraph(g[1], "significance");
184  confgraph(g[2], "signal efficiency");
185  confgraph(g[3], "background efficiency");
186 
187  TF1 *f1=new TF1("f1","0.5*[0]*(1.0-TMath::Erf((x-[1])/[2]))+[3]",0,100.);
188  //TF1 *f1=new TF1("f1","0.5*[0]*(1.0-TMath::Erf((x-[1])/[2]))+pol2(3)",0,100.);
189  //f1->SetParLimits(2,15,300);
190 
191  // -------------------------------------------------------------------
192  // 8x2 phi slices, symmetric
193  // -------------------------------------------------------------------
194 
195  // missing fraction in [%] of EMC = [ 1.0 - (tht_max - 22°)/118° ]* 100
196  double emc_rmv[9] = { 0., 12.5, 25., 37.5, 50., 62.5, 75., 87.5, 100. };
197 
198  //TString lab[4] = { ";supermodules missing;S/B", ";supermodules missing;significance [#sigma]", ";supermodules missing;signal efficiency [%]", ";supermodules missing;background efficiency [%]"};
199  TString lab[4] = { ";EMC missing [%];signal efficiency (#phi) [%]", ";EMC missing (#phi) [%];background efficiency [%]", ";EMC missing (#phi) [%];S/B", ";EMC missing (#phi) [%];significance [#sigma]"};
200 
201 
202  TCanvas *c1 = new TCanvas("c1","c1",1000,800);
203  c1->Divide(2,2,0.0001,0.0001);
204 
205  double sig_S = 20;
206  double sig_B = 42e6;
207 
208  double Lint = 1220*5;
209  double fBR = 0.06*0.339*0.394; // BR_J * BR_chic * BR_eta
210 
211  double S_dat = sig_S * fBR * Lint;
212  double B_dat = sig_B * Lint;
213 
214  cout <<"S:B = "<<S_dat/B_dat<<" mode_S="<<sigmode<<" mode_B="<<bkgmode<<endl;
215 
216  printf (" i | S0 | B0 | S | B | fS | fB | S* | B* |\n");
217  printf ("----+----------+------------+----------+--------+----------+-----------+---------+---------+\n");
218  for (int i=0;i<Nred;++i)
219  {
220  long int S0 = 1e5;
221  long int B0 =evcnts[sigmode+i+1000];
222 
223  //cout <<"S0 = "<<S0<<" B0 = "<<B0<<endl;
224 
225  double fS = S_dat/S0;
226  double fB = B_dat/B0;
227 
228  TString sigcut = Form("%s && mode==%d && xmct", cut.Data(), sigmode+i);
229  TString bkgcut = Form("%s && mode==%d", cut.Data(), bkgmode+i);
230 
231 
232  double S = (double) cntEvt(t, sigcut);
233  double B = (double) (mrgB.find(sigmode+i+1000) == mrgB.end() ? cntEvt(t, bkgcut)+1 : mrgB[sigmode+i+1000]);
234 
235  (double) cntEvt(t, bkgcut)+1;
236  double dS = sqrt(S);
237  double dB = sqrt(B);
238 
239  printf("%2d | %6ld | %6.3f G | %6d | %4d | %6.4f | %6.4f | %5.0f | %5.0f |\n", i, S0, (double)B0/1e9, (int)S, (int)B, fS, fB, S*fS, B*fB);
240  //cout <<i<<": S="<<S<<" B="<<B<<" --> S*"<<fS<<"="<<S*fS<<" B*"<<fB<<"="<<B*fB<<endl;
241 
242  double SN = S*fS/(B*fB);
243  double dSN = SN*sqrt(dS*dS/(S*S) + dB*dB/(B*B));
244 
245  g[2]->SetPoint(i, emc_rmv[i], SN);
246  g[2]->SetPointError(i, 0, dSN);
247 
248 
249  double Z = S*fS/sqrt(S*fS+B*fB);
250  double dZ = 0.5 * sqrt( ((fS*fS*S+2*fB*fS*B)*(fS*fS*S+2*fB*fS*B)*S + fS*fS*fB*fB*S*S*B)/pow(fS*S+fB*B,3)) ;//0.5 * sqrt( (pow(2*B*fB + S*fS,2)*S + fS*fS*S*S*B)/pow(fS*S+fB*B,3));
251 
252  g[3]->SetPoint(i, emc_rmv[i], Z);
253  g[3]->SetPointError(i, 0, dZ);
254 
255 
256  double effS = S/S0;
257  double deffS = effS*sqrt( 1./S + 1./S0 );
258 
259  g[0]->SetPoint(i, emc_rmv[i], effS*100.);
260  g[0]->SetPointError(i, 0, deffS*100.);
261 
262 
263  double effB = B/B0;
264  double deffB = effB*sqrt( 1./B + 1./B0 );
265 
266  g[1]->SetPoint(i, emc_rmv[i], effB*100.);
267  g[1]->SetPointError(i, 0, deffB*100.);
268  }
269 
270  TH1F *h[4];
271 
272  gStyle->SetOptFit(0);
273  gStyle->SetOptStat(0);
274 
275  TString plotnam[4]={"effs","effb","sn","sign"};
276 
277  for (int i=0;i<4;++i)
278  {
279  c1->cd(i+1);
280  gPad->SetTopMargin(0.10);
281  h[i] = createHistoGraph(g[i],tit+lab[i]);
282  h[i]->Draw();
283  g[i]->Draw("P same");
284 
285  f1->SetParameters(TMath::MaxElement(g[i]->GetN(), g[i]->GetY()),70, 10, TMath::MinElement(g[i]->GetN(), g[i]->GetY()));
286  f1->SetParLimits(2,50,200);
287 
288  //if (i==1)
289  //{
290  //f1->SetParameters(TMath::MaxElement(g[i]->GetN(), g[i]->GetY()),120, -10, TMath::MinElement(g[i]->GetN(), g[i]->GetY()));
291  //f1->SetParLimits(1,110,200);
292  //}
293 
294 
295  //if (i<3)
296  //if (i>0)
297  g[i]->Fit("f1","q");
298  //else g[i]->Fit("f1","q","",0,95);
299  //else g[i]->Fit("f2");
300  g[i]->SetName(plotnam[i]+"_"+((TString)fname(9,7)));
301  g[i]->Write();
302  }
303 
304  fname.ReplaceAll(".root","");
305  cut.ReplaceAll("&&","_AND_");
306  cut.ReplaceAll("||","_OR_");
307  cut.ReplaceAll("!","_NOT_");
308  cut.ReplaceAll(">","_lg_");
309  cut.ReplaceAll("<","_sm_");
310  cut.ReplaceAll(".","_");
311 
312  c1->SaveAs(Form("fig/%s__%03d__%s.gif",fname.Data(), sigmode, cut.Data()));
313  c1->SaveAs(Form("fig/%s__%03d__%s.C",fname.Data(), sigmode, cut.Data()));
314 
315  fana.Close();
316 }
317 
318 /*
319 Found 244 files with pattern data/DPMetac1_A_EMC1 containing event info Sum of events = 781589138
320 Found 243 files with pattern data/DPMetac1_A_EMC2 containing event info Sum of events = 778264034
321 Found 247 files with pattern data/DPMetac1_A_EMC3 containing event info Sum of events = 791030134
322 Found 244 files with pattern data/DPMetac1_A_EMC4 containing event info Sum of events = 781662143
323 Found 248 files with pattern data/DPMetac1_A_EMC5 containing event info Sum of events = 794307759
324 Found 243 files with pattern data/DPMetac1_A_EMC6 containing event info Sum of events = 778257467
325 Found 242 files with pattern data/DPMetac1_A_EMC7 containing event info Sum of events = 775037895
326 Found 242 files with pattern data/DPMetac1_A_EMC8 containing event info Sum of events = 775095496
327 --------------------
328 Found 241 files with pattern data/DPMetac1_B_EMC1 containing event info Sum of events = 771776262
329 Found 242 files with pattern data/DPMetac1_B_EMC2 containing event info Sum of events = 775016383
330 Found 244 files with pattern data/DPMetac1_B_EMC3 containing event info Sum of events = 781427165
331 Found 235 files with pattern data/DPMetac1_B_EMC4 containing event info Sum of events = 752754673
332 Found 247 files with pattern data/DPMetac1_B_EMC5 containing event info Sum of events = 791054235
333 Found 242 files with pattern data/DPMetac1_B_EMC6 containing event info Sum of events = 775064262
334 Found 241 files with pattern data/DPMetac1_B_EMC7 containing event info Sum of events = 771961964
335 Found 244 files with pattern data/DPMetac1_B_EMC8 containing event info Sum of events = 781628321
336 */
337 
338  //data/DPMetac1_A_SLC1 ( 82/ 100 files) = 1313354597
339  //data/DPMetac1_A_SLC2 ( 86/ 100 files) = 1377489004
340  //data/DPMetac1_A_SLC3 ( 85/ 100 files) = 1361147492
341  //data/DPMetac1_A_SLC4 ( 87/ 100 files) = 1393104477
342  //data/DPMetac1_A_SLC5 ( 86/ 100 files) = 1377270184
343  //data/DPMetac1_A_SLC6 ( 86/ 98 files) = 1376947715
344  //data/DPMetac1_A_SLC7 ( 84/ 96 files) = 1345215200
345  //data/DPMetac1_A_SLC8 ( 83/ 95 files) = 1328843629
346  //data/DPMetac1_A_SLC9 ( 83/ 100 files) = 1329393334
347  //--------------------
348  //data/DPMetac1_B_SLC1 ( 87/ 100 files) = 1393415089
349  //data/DPMetac1_B_SLC2 ( 89/ 100 files) = 1425022249
350  //data/DPMetac1_B_SLC3 ( 84/ 100 files) = 1345470462
351  //data/DPMetac1_B_SLC4 ( 82/ 99 files) = 1313202603
352  //data/DPMetac1_B_SLC5 ( 86/ 100 files) = 1377424566
353  //data/DPMetac1_B_SLC6 ( 81/ 97 files) = 1297113544
354  //data/DPMetac1_B_SLC7 ( 87/ 98 files) = 1393303719
355  //data/DPMetac1_B_SLC8 ( 84/ 100 files) = 1345397603
356  //data/DPMetac1_B_SLC9 ( 86/ 98 files) = 1377182415
357 
358  //std::map<long int,long int> evcnts =
359  //{ {1300, 1313354597}, {1301, 1377489004}, {1302, 1361147492}, {1303, 1393104477}, {1304, 1377270184}, {1305, 1376947715}, {1306, 1345215200}, {1307, 1328843629}, {1308, 1329393334}, // Setup A, pbp -> J/psi (-> e+ e-) pi+ pi-
360  //{1320, 1393415089}, {1321, 1425022249}, {1322, 1345470462}, {1323, 1313202603}, {1324, 1377424566}, {1325, 1297113544}, {1326, 1393303719}, {1327, 1345397603}, {1328, 1377182415}
361  //};
362 
printf
printf("RealTime=%f seconds, CpuTime=%f seconds\n", rtime, ctime)
t1
Int_t t1
Definition: hist-t7.C:106
createHistoGraph
TH1F * createHistoGraph(TGraph *g, TString tit="", double xmin=0, double xmax=0)
Definition: analyse_etac1_slc.C:60
i
Int_t i
Definition: run_full.C:25
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
confgraph
void confgraph(TGraph *g, TString tit, int col=1, int marker=20)
Definition: analyse_etac1_slc.C:47
cntEvt
long int cntEvt(TTree *t1, TString cut, TString var="ev")
Definition: analyse_etac1_slc.C:3
cut
double cut[MAX]
Definition: autocutx.C:36
analyse_etac1_slc
void analyse_etac1_slc(TString fname="ntp5_slc_etac1_A.root", TString cut="chi24c<50", int sigmode=300, TString tit="")
Definition: analyse_etac1_slc.C:91
f
TFile * f
Definition: bump_analys.C:12
c1
c1
Definition: plot_dirc.C:35
cnt
Int_t cnt
Definition: hist-t7.C:106
Z
double Z
Definition: anaLmdDigi.C:68
t
TTree * t
Definition: bump_analys.C:13
xmin
Double_t xmin
Definition: convolutionAnalysis.C:10
Generated on Sun Aug 18 2019 03:01:05 for FairRoot/PandaRoot by   doxygen 1.8.8