FairRoot/PandaRoot
analyse_J_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_J_slc(TString fname="ntp1_slc_Jee_A.root", TString cut="xd0d0pide>0.5&&xd0d1pide>0.5&&chi24c<50", int sigmode=100, TString tit="")
92 {
93  //data/DPMJee_A_SLC1 ( 366/ 500 files) = 29996980207
94  //data/DPMJee_A_SLC2 ( 491/ 616 files) = 27740885200
95  //data/DPMJee_A_SLC3 ( 460/ 604 files) = 26471309572
96  //data/DPMJee_A_SLC4 ( 491/ 617 files) = 27745096509
97  //data/DPMJee_A_SLC5 ( 485/ 624 files) = 27497769591
98  //data/DPMJee_A_SLC6 ( 513/ 624 files) = 28645444363
99  //data/DPMJee_A_SLC7 ( 564/ 659 files) = 25532064241
100  //data/DPMJee_A_SLC8 ( 594/ 662 files) = 24338865274
101  //data/DPMJee_A_SLC9 ( 502/ 503 files) = 20568203496
102  //--------------------
103  //data/DPMJee_B_SLC1 ( 297/ 500 files) = 24341219982
104  //data/DPMJee_B_SLC2 ( 659/ 687 files) = 27009017984
105  //data/DPMJee_B_SLC3 ( 656/ 682 files) = 26881191552
106  //data/DPMJee_B_SLC4 ( 613/ 673 files) = 25119567078
107  //data/DPMJee_B_SLC5 ( 612/ 665 files) = 25082906031
108  //data/DPMJee_B_SLC6 ( 614/ 616 files) = 25163739980
109  //data/DPMJee_B_SLC7 ( 608/ 610 files) = 24914199407
110  //data/DPMJee_B_SLC8 ( 511/ 512 files) = 20942409620
111  //data/DPMJee_B_SLC9 ( 507/ 507 files) = 20772725573
112  //--------------------
113  //--------------------
114  //data/DPMJmm_A_SLC1 ( 272/ 273 files) = 20255703739
115  //data/DPMJmm_A_SLC2 ( 200/ 200 files) = 7446866207
116  //data/DPMJmm_A_SLC3 ( 214/ 214 files) = 7970194795
117  //data/DPMJmm_A_SLC4 ( 202/ 202 files) = 7520702389
118  //data/DPMJmm_A_SLC5 ( 214/ 214 files) = 7965350396
119  //data/DPMJmm_A_SLC6 ( 250/ 250 files) = 9310510861
120  //data/DPMJmm_A_SLC7 ( 244/ 244 files) = 9084144004
121  //data/DPMJmm_A_SLC8 ( 247/ 247 files) = 9197632458
122  //data/DPMJmm_A_SLC9 ( 233/ 233 files) = 8676693188
123  //--------------------
124  //data/DPMJmm_B_SLC1 ( 101/ 101 files) = 7519725318
125  //data/DPMJmm_B_SLC2 ( 231/ 231 files) = 8595694989
126  //data/DPMJmm_B_SLC3 ( 235/ 235 files) = 8751031806
127  //data/DPMJmm_B_SLC4 ( 253/ 253 files) = 9421220844
128  //data/DPMJmm_B_SLC5 ( 242/ 242 files) = 9010126829
129  //data/DPMJmm_B_SLC6 ( 273/ 273 files) = 10164780274
130  //data/DPMJmm_B_SLC7 ( 239/ 239 files) = 8896609302
131  //data/DPMJmm_B_SLC8 ( 248/ 248 files) = 9235338848
132  //data/DPMJmm_B_SLC9 ( 239/ 239 files) = 8900632658
133 
134  std::map<long int,long int> evcnts =
135  { {1100, 29996980207}, {1101, 27740885200}, {1102, 26471309572}, {1103, 27745096509}, {1104, 27497769591}, {1105, 28645444363}, {1106, 25532064241}, {1107, 24338865274}, {1108, 20568203496}, // Setup A, pbp -> J/psi (-> e+ e-) pi+ pi-
136  {1120, 24341219982}, {1121, 27009017984}, {1122, 26881191552 }, {1123, 25119567078}, {1124, 25082906031}, {1125, 25163739980}, {1126, 24914199407}, {1127, 20942409620}, {1128, 20772725573}, // Setup B, pbp -> J/psi (-> e+ e-) pi+ pi- //FIXME modes 1120 ...
137 
138  {1200, 20255703739}, {1201, 7446866207}, {1202, 7970194795}, {1203, 7520702389}, {1204, 7965350396}, {1205, 9310510861}, {1206, 8972578036}, {1207, 9197632458}, {1208, 8676693188}, // Setup A, pbp -> J/psi (-> mu+ mu-) pi+ pi-
139  {1220, 7519725318}, {1221, 8595694989}, {1222, 8751031806}, {1223, 9421220844}, {1224, 8935730839}, {1225, 9084618361}, {1226, 8896609302}, {1227, 9235338848}, {1228, 8900632658} // Setup B, pbp -> J/psi (-> mu+ mu-) pi+ pi- //FIXME modes 1220...
140  };
141 
142 
143 
144  //for ( auto x:evcnts) cout <<"mode "<<x.first<<" : "<<x.second<<endl;
145  //cout <<endl;
146 
147  int bkgmode = sigmode+1000;
148  //if (sigmode%100>19) bkgmode+=1;
149 
150  TFile *f = new TFile(fname);
151  TTree *t = (TTree*)f->Get("ntp1");
152 
153  TFile fana("anaJ_slc.root","UPDATE");
154 
155  int Nred = 9;
156 
157  TGraphErrors *g[4];
158  g[0] = new TGraphErrors(Nred);
159  g[1] = new TGraphErrors(Nred);
160  g[2] = new TGraphErrors(Nred);
161  g[3] = new TGraphErrors(Nred);
162 
163  confgraph(g[0], "signal to noise");
164  confgraph(g[1], "significance");
165  confgraph(g[2], "signal efficiency");
166  confgraph(g[3], "background efficiency");
167 
168  TF1 *f1[4];
169 
170  f1[0]=new TF1("f0","0.5*[0]*(1.0-TMath::Erf((x-[1])/[2]))+[3]",0,100.);
171  f1[0]->SetParameters(10,70, 10, 5);
172 
173  f1[1]=new TF1("f1","0.5*[0]*(1.0-TMath::Erf((x-[1])/[2]))+[3]",0,100.);
174  //1 p0 9.25779e-05 1.53699e-06 1.20127e-09 1.87008e+02
175  //2 p1 5.00805e+01 2.94904e+00 6.21765e-06 -1.66797e-02
176  //3 p2 -6.50046e+01 2.02989e+00 2.97643e-04** at limit **
177  //4 p3 6.43684e-05 1.96720e-06 6.21319e-10 1.95246e+02
178  f1[1]->SetParameters(1e-4, 50, -65, 6e-5);
179  f1[1]->SetParLimits(1,20,200);
180  f1[1]->SetParLimits(2,20,200);
181 
182  //f1[2]=new TF1("f2","0.5*[0]*(1.0-TMath::Erf((x-[1])/[2]))+[3]",0,100.);
183  //f1[2]->SetParameters(1e-7, 120, -10, 1e-8);
184  //f1[2]->SetParLimits(1,20,200);
185 
186  f1[2]=new TF1("f2","expo(0)+pol1(2)",0,100.);
187  f1[2]->SetParameters(-0.5,-0.01,1,1);
188 
189  f1[3]=new TF1("f3","expo(0)+[2]",0,100.);
190  f1[3]->SetParameters(1.2,-0.02,1.6); //
191 
192 
193  //f1->SetParLimits(2,15,100);
194  //TF1 *f1=new TF1("f1","pol2(0)",0,100.);
195  //f1->SetParameters(1,1,1);
196 
197  // -------------------------------------------------------------------
198  // 8x2 phi slices, symmetric
199  // -------------------------------------------------------------------
200 
201  // missing fraction in [%] of EMC = [ 1.0 - (tht_max - 22°)/118° ]* 100
202  double emc_rmv[9] = { 0., 12.5, 25., 37.5, 50., 62.5, 75., 87.5, 100. };
203 
204  //TString lab[4] = { ";supermodules missing;S/B", ";supermodules missing;significance [#sigma]", ";supermodules missing;signal efficiency [%]", ";supermodules missing;background efficiency [%]"};
205  TString lab[4] = { ";EMC missing [%];signal efficiency (#phi) [%]", ";EMC missing (#phi) [%];background efficiency [%]", ";EMC missing (#phi) [%];S/B", ";EMC missing (#phi) [%];significance [#sigma]"};
206 
207 
208  TCanvas *c1 = new TCanvas("c1","c1",1000,800);
209  c1->Divide(2,2,0.0001,0.0001);
210 
211  double sig_S = 50;
212  double sig_B = 46e6;
213 
214  double Lint = 1170*2;
215  double fBR = 0.05*0.06;
216 
217  double S_dat = sig_S * fBR * Lint;
218  double B_dat = sig_B * Lint;
219 
220  cout <<"S:B = "<<S_dat/B_dat<<" mode_S="<<sigmode<<" mode_B="<<bkgmode<<endl;
221 
222  printf (" i | S0 | B0 | S | B | fS | fB | S* | B* |\n");
223  printf ("----+----------+------------+----------+---------+----------+-----------+---------+----------+\n");
224  for (int i=0;i<Nred;++i)
225  {
226  long int S0 = 1e5;
227  long int B0 = evcnts[sigmode+i+1000];
228 
229  //cout <<"S0 = "<<S0<<" B0 = "<<B0<<endl;
230 
231  double fS = S_dat/S0;
232  double fB = B_dat/B0;
233 
234  TString sigcut = Form("%s && mode==%d", cut.Data(), sigmode+i);
235  TString bkgcut = Form("%s && mode==%d", cut.Data(), bkgmode+i);
236 
237 
238  double S = (double) cntEvt(t, sigcut);
239  double B = (double) cntEvt(t, bkgcut);
240  if (B==0) B=1;
241  double dS = sqrt(S);
242  double dB = sqrt(B);
243 
244  printf("%2d | %6ld | %6.3f G | %6d | %5d | %6.4f | %7.4f | %5.0f | %6.0f |\n", i, S0, (double)B0/1e9, (int)S, (int)B, fS, fB, S*fS, B*fB);
245  //cout <<i<<": S="<<S<<" B="<<B<<" --> S*"<<fS<<"="<<S*fS<<" B*"<<fB<<"="<<B*fB<<endl;
246 
247  double SN = S*fS/(B*fB);
248  double dSN = SN*sqrt(dS*dS/(S*S) + dB*dB/(B*B));
249 
250  g[2]->SetPoint(i, emc_rmv[i], SN);
251  g[2]->SetPointError(i, 0, dSN);
252 
253 
254  double Z = S*fS/sqrt(S*fS+B*fB);
255  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));
256 
257  g[3]->SetPoint(i, emc_rmv[i], Z);
258  g[3]->SetPointError(i, 0, dZ);
259 
260 
261  double effS = S/S0;
262  double deffS = effS*sqrt( 1./S + 1./S0 );
263 
264  g[0]->SetPoint(i, emc_rmv[i], effS*100.);
265  g[0]->SetPointError(i, 0, deffS*100.);
266 
267 
268  double effB = B/B0;
269  double deffB = effB*sqrt( 1./B + 1./B0 );
270 
271  g[1]->SetPoint(i, emc_rmv[i], effB*100.);
272  g[1]->SetPointError(i, 0, deffB*100.);
273  }
274 
275  TH1F *h[4];
276 
277  gStyle->SetOptFit(0);
278  gStyle->SetOptStat(0);
279 
280  TString plotnam[4]={"effs","effb","sn","sign"};
281 
282  for (int i=0;i<4;++i)
283  {
284  c1->cd(i+1);
285  h[i] = createHistoGraph(g[i],tit+lab[i]);
286  h[i]->Draw();
287  g[i]->Draw("P same");
288 
289  //f1[i]->SetParameters(TMath::MaxElement(g[i]->GetN(), g[i]->GetY()),70, 10, TMath::MinElement(g[i]->GetN(), g[i]->GetY()));
290  //if (i==1)
291  //{
292  //f1[i]->SetParameters(TMath::MaxElement(g[i]->GetN(), g[i]->GetY()),120, -10, TMath::MinElement(g[i]->GetN(), g[i]->GetY()));
293  //f1[i]->SetParLimits(1,20,200);
294  //}
295 
296  g[i]->Fit(Form("f%d",i),"q");
297 
298  g[i]->SetName(plotnam[i]+"_"+((TString)fname(9,5)));
299  g[i]->Write();
300  }
301 
302  fname.ReplaceAll(".root","");
303  cut.ReplaceAll("&&","_AND_");
304  cut.ReplaceAll("||","_OR_");
305  cut.ReplaceAll("!","_NOT_");
306  cut.ReplaceAll(">","_lg_");
307  cut.ReplaceAll("<","_sm_");
308  cut.ReplaceAll(".","_");
309 
310  c1->SaveAs(Form("fig/%s__%03d__%s.gif",fname.Data(), sigmode, cut.Data()));
311  c1->SaveAs(Form("fig/%s__%03d__%s.C",fname.Data(), sigmode, cut.Data()));
312 
313  fana.Close();
314 }
315 
316 
317  //data/DPMJee_A_SLC1 ( 366/ 500 files) = 29996980207
318  //data/DPMJee_A_SLC2 ( 186/ 500 files) = 15245616034
319  //data/DPMJee_A_SLC3 ( 186/ 500 files) = 15245616034
320  //data/DPMJee_A_SLC4 ( 186/ 500 files) = 15245616034
321  //data/DPMJee_A_SLC5 ( 186/ 500 files) = 15245616034
322  //data/DPMJee_A_SLC6 ( 186/ 500 files) = 15245616034
323  //data/DPMJee_A_SLC7 ( 59/ 500 files) = 4836581051
324  //data/DPMJee_A_SLC8 ( 185/ 500 files) = 7580713308
325  //data/DPMJee_A_SLC9 ( 211/ 211 files) = 8644837141
326  //--------------------
327  //data/DPMJee_B_SLC1 ( 297/ 500 files) = 24341219982
328  //data/DPMJee_B_SLC2 ( 216/ 500 files) = 8855594650
329  //data/DPMJee_B_SLC3 ( 208/ 500 files) = 8524796869
330  //data/DPMJee_B_SLC4 ( 182/ 500 files) = 7456922637
331  //data/DPMJee_B_SLC5 ( 224/ 500 files) = 9182421996
332  //data/DPMJee_B_SLC6 ( 205/ 205 files) = 8400473125
333  //data/DPMJee_B_SLC7 ( 207/ 208 files) = 8484522825
334  //data/DPMJee_B_SLC8 ( 196/ 196 files) = 8034610317
335  //data/DPMJee_B_SLC9 ( 190/ 190 files) = 7785364661
336 
337  //--------------------
338  //--------------------
339  //data/DPMJmm_A_SLC1 ( 272/ 273 files) = 20255703739
340  //data/DPMJmm_A_SLC2 ( 200/ 200 files) = 7446866207
341  //data/DPMJmm_A_SLC3 ( 214/ 214 files) = 7970194795
342  //data/DPMJmm_A_SLC4 ( 202/ 202 files) = 7520702389
343  //data/DPMJmm_A_SLC5 ( 214/ 214 files) = 7965350396
344  //data/DPMJmm_A_SLC6 ( 250/ 250 files) = 9310510861
345  //data/DPMJmm_A_SLC7 ( 241/ 241 files) = 8972578036
346  //data/DPMJmm_A_SLC8 ( 247/ 247 files) = 9197632458
347  //data/DPMJmm_A_SLC9 ( 233/ 233 files) = 8676693188
348  //--------------------
349  //data/DPMJmm_B_SLC1 ( 101/ 101 files) = 7519725318
350  //data/DPMJmm_B_SLC2 ( 231/ 231 files) = 8595694989
351  //data/DPMJmm_B_SLC3 ( 235/ 235 files) = 8751031806
352  //data/DPMJmm_B_SLC4 ( 253/ 253 files) = 9421220844
353  //data/DPMJmm_B_SLC5 ( 240/ 240 files) = 8935730839
354  //data/DPMJmm_B_SLC6 ( 244/ 244 files) = 9084618361
355  //data/DPMJmm_B_SLC7 ( 239/ 239 files) = 8896609302
356  //data/DPMJmm_B_SLC8 ( 248/ 248 files) = 9235338848
357  //data/DPMJmm_B_SLC9 ( 239/ 239 files) = 8900632658
358 
359 
360  //std::map<long int,long int> evcnts =
361  //{ {1100, 29996980207}, {1101, 15245616034}, {1102, 15245616034}, {1103, 15245616034}, {1104, 15245616034}, {1105, 15245616034}, {1106, 4836581051}, {1107, 7580713308}, {1108, 8644837141}, // Setup A, pbp -> J/psi (-> e+ e-) pi+ pi-
362  //{1120, 24341219982}, {1121, 8855594650}, {1122, 8524796869 }, {1123, 7456922637}, {1124, 9182421996}, {1125, 8400473125}, {1126, 8484522825}, {1127, 8034610317}, {1128, 7785364661}, // Setup B, pbp -> J/psi (-> e+ e-) pi+ pi- //FIXME modes 1120 ...
363 
364  //{1200, 1}, {1201, 1}, {1202, 1}, {1203, 1}, {1204, 1}, {1205, 1}, {1206, 1}, {1207, 1}, {1208, 1}, // Setup A, pbp -> J/psi (-> mu+ mu-) pi+ pi-
365  //{1220, 1}, {1221, 1}, {1222, 1}, {1223, 1}, {1224, 1}, {1225, 1}, {1226, 1}, {1227, 1}, {1228, 1} // Setup B, pbp -> J/psi (-> mu+ mu-) pi+ pi- //FIXME modes 1220...
366  //};
367 
368  //std::map<long int,long int> evcnts =
369  //{ {1100, 29996980207}, {1101, 15245616034}, {1102, 15245616034}, {1103, 15245616034}, {1104, 15245616034}, {1105, 15245616034}, {1106, 4836581051}, {1107, 7580713308}, {1108, 8644837141}, // Setup A, pbp -> J/psi (-> e+ e-) pi+ pi-
370  //{1120, 24341219982}, {1121, 8855594650}, {1122, 8524796869 }, {1123, 7456922637}, {1124, 9182421996}, {1125, 8400473125}, {1126, 8484522825}, {1127, 8034610317}, {1128, 7785364661}, // Setup B, pbp -> J/psi (-> e+ e-) pi+ pi- //FIXME modes 1120 ...
371 
372  //{1200, 20255703739}, {1201, 7446866207}, {1202, 7970194795}, {1203, 7520702389}, {1204, 7965350396}, {1205, 9310510861}, {1206, 8972578036}, {1207, 9197632458}, {1208, 8676693188}, // Setup A, pbp -> J/psi (-> mu+ mu-) pi+ pi-
373  //{1220, 7519725318}, {1221, 8595694989}, {1222, 8751031806}, {1223, 9421220844}, {1224, 8935730839}, {1225, 9084618361}, {1226, 8896609302}, {1227, 9235338848}, {1228, 8900632658} // Setup B, pbp -> J/psi (-> mu+ mu-) pi+ pi- //FIXME modes 1220...
374  //};
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
analyse_J_slc
void analyse_J_slc(TString fname="ntp1_slc_Jee_A.root", TString cut="xd0d0pide>0.5&&xd0d1pide>0.5&&chi24c<50", int sigmode=100, TString tit="")
Definition: analyse_J_slc.C:91
Generated on Wed Apr 3 2019 08:43:48 for FairRoot/PandaRoot by   doxygen 1.8.5