FairRoot/PandaRoot
analyse_etac1.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(TString fname="ntp5_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 
120 
121  //data/DPMetac1_A_EMC1 ( 179/ 199 files) = 2866079091
122  //data/DPMetac1_A_EMC2 ( 173/ 198 files) = 2770144980
123  //data/DPMetac1_A_EMC3 ( 166/ 196 files) = 2658679173
124  //data/DPMetac1_A_EMC4 ( 178/ 198 files) = 2850295200
125  //data/DPMetac1_A_EMC5 ( 179/ 200 files) = 2866373874
126  //data/DPMetac1_A_EMC6 ( 169/ 198 files) = 2706157924
127  //data/DPMetac1_A_EMC7 ( 175/ 197 files) = 2801993946
128  //data/DPMetac1_A_EMC8 ( 169/ 199 files) = 2706466096
129  //--------------------
130  //data/DPMetac1_B_EMC1 ( 183/ 198 files) = 2930718917
131  //data/DPMetac1_B_EMC2 ( 168/ 199 files) = 2690367302
132  //data/DPMetac1_B_EMC3 ( 167/ 197 files) = 2674423426
133  //data/DPMetac1_B_EMC4 ( 177/ 200 files) = 2834563180
134  //data/DPMetac1_B_EMC5 ( 174/ 199 files) = 2786239006
135  //data/DPMetac1_B_EMC6 ( 169/ 198 files) = 2706374195
136  //data/DPMetac1_B_EMC7 ( 173/ 198 files) = 2770320696
137  //data/DPMetac1_B_EMC8 ( 171/ 198 files) = 2738620309
138 
139 
140 
141  std::map<long int,long int> evcnts =
142  { {1300, 2866079091}, {1301, 2770144980}, {1302, 2658679173}, {1303, 2850295200}, {1304, 2866373874}, {1305, 2706157924}, {1306, 2801993946}, {1307, 2706466096}, // Setup A, pbp -> J/psi (-> e+ e-) pi+ pi-
143  {1320, 2930718917}, {1321, 2690367302}, {1322, 2674423426}, {1323, 2834563180}, {1324, 2786239006}, {1325, 2706374195}, {1326, 2770320696}, {1327, 2738620309}
144  };
145 
146  // we merge the values at 0% and 100% from both studies by hand...
147  // **** from phi study
148  //data/DPMetac1_A_SLC1 ( 168/ 200 files) = 2690752186
149  //data/DPMetac1_A_SLC9 ( 168/ 200 files) = 2690419363
150  //--------------------
151  //data/DPMetac1_B_SLC1 ( 167/ 200 files) = 2674554579
152  //data/DPMetac1_B_SLC9 ( 176/ 200 files) = 2818616809
153 
154  evcnts[1300] += 2690752186;
155  evcnts[1320] += 2690419363;
156  evcnts[1307] += 2674554579;
157  evcnts[1327] += 2818616809;
158 
159  std::map<long int,long int> mrgB =
160  { {1300, 17+12}, {1307, 4+1},
161  {1320, 11+14}, {1327, 1+1}
162  };
163 
164 
165  //for ( auto x:evcnts) cout <<"mode "<<x.first<<" : "<<x.second<<endl;
166  //cout <<endl;
167 
168  int bkgmode = sigmode+1000;
169  //if (sigmode%100>19) bkgmode+=1;
170 
171  TFile *f = new TFile(fname);
172  TTree *t = (TTree*)f->Get("ntp5");
173 
174  TFile fana("ana_etac.root","UPDATE");
175 
176  TGraphErrors *g[4];
177  g[0] = new TGraphErrors(8);
178  g[1] = new TGraphErrors(8);
179  g[2] = new TGraphErrors(8);
180  g[3] = new TGraphErrors(8);
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  // Supermodul | 1 | 2 | 3 | 4 | 5 | 6 | 7 |
193  // num Alveolen in theta | 1 | 3 | 3 | 3 | 3 | 3 | 2 |
194  // coverage from 22 to |140.0|133.4|113.8| 94.1| 74.4| 54.8| 35.1|
195  // -------------------------------------------------------------------
196 
197  // missing fraction in [%] of EMC = [ 1.0 - (tht_max - 22°)/118° ]* 100
198  double emc_rmv[8] = { 0., 6., 23., 39., 56., 72., 89., 100.};
199 
200  //TString lab[4] = { ";supermodules missing;S/B", ";supermodules missing;significance [#sigma]", ";supermodules missing;signal efficiency [%]", ";supermodules missing;background efficiency [%]"};
201  TString lab[4] = { ";EMC missing (polar #theta) [%];signal efficiency [%]", ";EMC missing (polar #theta) [%];background efficiency [%]", ";EMC missing (polar #theta) [%];S/B", ";EMC missing (polar #theta) [%];significance [#sigma]"};
202 
203 
204  TCanvas *c1 = new TCanvas("c1","c1",1000,800);
205  c1->Divide(2,2,0.0001,0.0001);
206 
207  double sig_S = 20;
208  double sig_B = 42e6;
209 
210  double Lint = 1220*5;
211  double fBR = 0.06*0.339*0.394; // BR_J * BR_chic * BR_eta
212 
213  double S_dat = sig_S * fBR * Lint;
214  double B_dat = sig_B * Lint;
215 
216  cout <<"S:B = "<<S_dat/B_dat<<" mode_S="<<sigmode<<" mode_B="<<bkgmode<<endl;
217 
218  printf (" i | S0 | B0 | S | B | fS | fB | S* | B* |\n");
219  printf ("----+----------+------------+----------+--------+----------+-----------+---------+---------+\n");
220  for (int i=0;i<8;++i)
221  {
222  long int S0 = 1e5;
223  long int B0 = evcnts[sigmode+i+1000];
224 
225 
226  double fS = S_dat/S0;
227  double fB = B_dat/B0;
228 
229  TString sigcut = Form("%s && mode==%d && xmct", cut.Data(), sigmode+i);
230  TString bkgcut = Form("%s && mode==%d", cut.Data(), bkgmode+i);
231 
232 
233  double S = (double) cntEvt(t, sigcut);
234  double B = (double) (mrgB.find(sigmode+i+1000) == mrgB.end() ? cntEvt(t, bkgcut)+1 : mrgB[sigmode+i+1000]);
235  double dS = sqrt(S);
236  double dB = sqrt(B);
237 
238  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);
239  //cout <<i<<" : " <<"S0 = "<<S0<<" B0 = "<<B0<<" : S="<<S<<" B="<<B<<" --> S*"<<fS<<"="<<S*fS<<" B*"<<fB<<"="<<B*fB<<endl;
240 
241  double SN = S*fS/(B*fB);
242  double dSN = SN*sqrt(dS*dS/(S*S) + dB*dB/(B*B));
243 
244  g[2]->SetPoint(i, emc_rmv[i], SN);
245  g[2]->SetPointError(i, 0, dSN);
246 
247 
248  double Z = S*fS/sqrt(S*fS+B*fB);
249  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));
250 
251  g[3]->SetPoint(i, emc_rmv[i], Z);
252  g[3]->SetPointError(i, 0, dZ);
253 
254 
255  double effS = S/S0;
256  double deffS = effS*sqrt( 1./S + 1./S0 );
257 
258  g[0]->SetPoint(i, emc_rmv[i], effS*100.);
259  g[0]->SetPointError(i, 0, deffS*100.);
260 
261 
262  double effB = B/B0;
263  double deffB = effB*sqrt( 1./B + 1./B0 );
264 
265  g[1]->SetPoint(i, emc_rmv[i], effB*100.);
266  g[1]->SetPointError(i, 0, deffB*100.);
267  }
268 
269  TH1F *h[4];
270 
271  gStyle->SetOptFit(0);
272  gStyle->SetOptStat(0);
273 
274  TString plotnam[4]={"effs","effb","sn","sign"};
275 
276  double hmaxy[4] = {20., 1e-5, 0.25, 7.};
277 
278  for (int i=0;i<4;++i)
279  {
280  c1->cd(i+1);
281  gPad->SetTopMargin(0.10);
282  h[i] = createHistoGraph(g[i],tit+lab[i]);
283  h[i]->SetMaximum(hmaxy[i]);
284  h[i]->Draw();
285  g[i]->Draw("P same");
286 
287  f1->SetParameters(TMath::MaxElement(g[i]->GetN(), g[i]->GetY()),70, 10, TMath::MinElement(g[i]->GetN(), g[i]->GetY()));
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) g[i]->Fit("f1","q");
297  else g[i]->Fit("f1","q","",0,95);
298  //else g[i]->Fit("f2");
299  g[i]->SetName(plotnam[i]+"_"+((TString)fname(5,7)));
300  g[i]->Write();
301  }
302 
303  fname.ReplaceAll(".root","");
304  cut.ReplaceAll("&&","_AND_");
305  cut.ReplaceAll("||","_OR_");
306  cut.ReplaceAll("!","_NOT_");
307  cut.ReplaceAll(">","_lg_");
308  cut.ReplaceAll("<","_sm_");
309  cut.ReplaceAll(".","_");
310 
311  c1->SaveAs(Form("fig/%s__%03d__%s.gif",fname.Data(), sigmode, cut.Data()));
312  c1->SaveAs(Form("fig/%s__%03d__%s.C",fname.Data(), sigmode, cut.Data()));
313 
314  fana.Close();
315 }
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  std::map<long int,long int> evcnts =
338  { {1300,781589138}, {1301,778264034}, {1302,791030134}, {1303,781662143}, {1304,794307759}, {1305,778257467}, {1306,775037895}, {1307,775095496}, // Setup A, pbp -> J/psi (-> e+ e-) pi+ pi-
339  {1320,771776262}, {1321,775016383}, {1322,781427165}, {1323,752754673}, {1324,791054235}, {1325,775064262}, {1326,771961964}, {1327,781628321}
340  };
341 */
342 
343  //data/DPMetac1_A_EMC1 ( 451/ 504 files) = 4096546385
344  //data/DPMetac1_A_EMC2 ( 460/ 500 files) = 4252979964
345  //data/DPMetac1_A_EMC3 ( 464/ 506 files) = 4266061680
346  //data/DPMetac1_A_EMC4 ( 458/ 500 files) = 4208139218
347  //data/DPMetac1_A_EMC5 ( 467/ 500 files) = 4301265578
348  //data/DPMetac1_A_EMC6 ( 459/ 497 files) = 4236964259
349  //data/DPMetac1_A_EMC7 ( 456/ 499 files) = 4202381786
350  //data/DPMetac1_A_EMC8 ( 458/ 499 files) = 4233827498
351  //--------------------
352  //data/DPMetac1_B_EMC1 ( 455/ 498 files) = 4198373087
353  //data/DPMetac1_B_EMC2 ( 454/ 499 files) = 4169903068
354  //data/DPMetac1_B_EMC3 ( 459/ 499 files) = 4224726395
355  //data/DPMetac1_B_EMC4 ( 451/ 497 files) = 4211893563
356  //data/DPMetac1_B_EMC5 ( 462/ 498 files) = 4234207999
357  //data/DPMetac1_B_EMC6 ( 455/ 499 files) = 4185698000
358  //data/DPMetac1_B_EMC7 ( 452/ 500 files) = 4150810055
359  //data/DPMetac1_B_EMC8 ( 459/ 500 files) = 4224406655
360 
361  //std::map<long int,long int> evcnts =
362  //{ {1300, 4096546385}, {1301, 4252979964}, {1302, 4266061680}, {1303, 4208139218}, {1304, 4301265578}, {1305, 4236964259}, {1306, 4202381786}, {1307, 4233827498}, // Setup A, pbp -> J/psi (-> e+ e-) pi+ pi-
363  //{1320, 4198373087}, {1321, 4169903068}, {1322, 4224726395}, {1323, 4211893563}, {1324, 4234207999}, {1325, 4185698000}, {1326, 4150810055}, {1327, 4224406655}
364  //};
365 
366 
367  //data/DPMetac1_A_EMC1 ( 630/ 699 files) = 6962622485
368  //data/DPMetac1_A_EMC2 ( 633/ 695 files) = 7023109828
369  //data/DPMetac1_A_EMC3 ( 631/ 696 files) = 6940762184
370  //data/DPMetac1_A_EMC4 ( 636/ 698 files) = 7058415189
371  //data/DPMetac1_A_EMC5 ( 646/ 700 files) = 7167639452
372  //data/DPMetac1_A_EMC6 ( 628/ 695 files) = 6943108546
373  //data/DPMetac1_A_EMC7 ( 628/ 692 files) = 6956329888
374  //data/DPMetac1_A_EMC8 ( 625/ 694 files) = 6908323139
375  //--------------------
376  //data/DPMetac1_B_EMC1 ( 637/ 695 files) = 7113112437
377  //data/DPMetac1_B_EMC2 ( 622/ 698 files) = 6860231634
378  //data/DPMetac1_B_EMC3 ( 626/ 696 files) = 6899184216
379  //data/DPMetac1_B_EMC4 ( 628/ 697 files) = 7046479715
380  //data/DPMetac1_B_EMC5 ( 636/ 697 files) = 7020440435
381  //data/DPMetac1_B_EMC6 ( 624/ 697 files) = 6892077970
382  //data/DPMetac1_B_EMC7 ( 625/ 698 files) = 6921142848
383  //data/DPMetac1_B_EMC8 ( 627/ 695 files) = 6914980193
384 
385  //std::map<long int,long int> evcnts =
386  //{ {1300, 6962622485}, {1301, 7023109828}, {1302, 6940762184}, {1303, 7058415189}, {1304, 7167639452}, {1305, 6943108546}, {1306, 6956329888}, {1307, 6908323139}, // Setup A, pbp -> J/psi (-> e+ e-) pi+ pi-
387  //{1320, 7113112437}, {1321, 6860231634}, {1322, 6899184216}, {1323, 7046479715}, {1324, 7020440435}, {1325, 6892077970}, {1326, 6921142848}, {1327, 6914980193}
388  //};
analyse_etac1
void analyse_etac1(TString fname="ntp5_etac1_A.root", TString cut="chi24c<50", int sigmode=300, TString tit="")
Definition: analyse_etac1.C:91
printf
printf("RealTime=%f seconds, CpuTime=%f seconds\n", rtime, ctime)
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
cntEvt
long int cntEvt(TTree *t1, TString cut, TString var="ev")
Definition: analyse_etac1.C:3
xmax
Double_t xmax
Definition: convolutionAnalysis.C:11
col
int col
Definition: anaLmdDigi.C:67
g
TFile * g
Definition: GammaSpectraAnalysis.C:22
createHistoGraph
TH1F * createHistoGraph(TGraph *g, TString tit="", double xmin=0, double xmax=0)
Definition: analyse_etac1.C:60
t1
Int_t t1
Definition: hist-t7.C:106
cut
double cut[MAX]
Definition: autocutx.C:36
confgraph
void confgraph(TGraph *g, TString tit, int col=1, int marker=20)
Definition: analyse_etac1.C:47
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 Wed Apr 3 2019 08:43:48 for FairRoot/PandaRoot by   doxygen 1.8.5