FairRoot/PandaRoot
analyse_J.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(TString fname="ntp1_Jee_A.root", TString cut="xd0d0pide>0.5&&xd0d1pide>0.5&&chi24c<50", int sigmode=100, TString tit="")
92 {
93  //data/DPMJee_A_EMC1 ( 421/ 423 files) = 24342651889
94  //data/DPMJee_A_EMC2 ( 412/ 413 files) = 23723523588
95  //data/DPMJee_A_EMC3 ( 392/ 394 files) = 21921465124
96  //data/DPMJee_A_EMC4 ( 382/ 386 files) = 21192161015
97  //data/DPMJee_A_EMC5 ( 379/ 380 files) = 20899111775
98  //data/DPMJee_A_EMC6 ( 382/ 383 files) = 21264224163
99  //data/DPMJee_A_EMC7 ( 391/ 394 files) = 21918977223
100  //data/DPMJee_A_EMC8 ( 371/ 375 files) = 20327238649
101  //--------------------
102  //data/DPMJee_B_EMC1 ( 383/ 386 files) = 21266925718
103  //data/DPMJee_B_EMC2 ( 388/ 392 files) = 21721036293
104  //data/DPMJee_B_EMC3 ( 381/ 383 files) = 21144322234
105  //data/DPMJee_B_EMC4 ( 373/ 375 files) = 20364376890
106  //data/DPMJee_B_EMC5 ( 379/ 379 files) = 20899648526
107  //data/DPMJee_B_EMC6 ( 387/ 389 files) = 21667979271
108  //data/DPMJee_B_EMC7 ( 393/ 397 files) = 22088200372
109  //data/DPMJee_B_EMC8 ( 383/ 386 files) = 21308984620
110  //--------------------
111  //--------------------
112  //data/DPMJmm_A_EMC1 ( 403/ 406 files) = 20850842535
113  //data/DPMJmm_A_EMC2 ( 398/ 399 files) = 20440294344
114  //data/DPMJmm_A_EMC3 ( 392/ 394 files) = 20333763226
115  //data/DPMJmm_A_EMC4 ( 361/ 362 files) = 18321499510
116  //data/DPMJmm_A_EMC5 ( 343/ 343 files) = 16868152786
117  //data/DPMJmm_A_EMC6 ( 336/ 338 files) = 16530666389
118  //data/DPMJmm_A_EMC7 ( 329/ 329 files) = 15753960106
119  //data/DPMJmm_A_EMC8 ( 316/ 317 files) = 14926218316
120  //--------------------
121  //data/DPMJmm_B_EMC1 ( 332/ 333 files) = 15939265122
122  //data/DPMJmm_B_EMC2 ( 330/ 333 files) = 15935942596
123  //data/DPMJmm_B_EMC3 ( 327/ 328 files) = 15490042463
124  //data/DPMJmm_B_EMC4 ( 237/ 238 files) = 8823408092
125  //data/DPMJmm_B_EMC5 ( 236/ 237 files) = 8789625908
126  //data/DPMJmm_B_EMC6 ( 236/ 236 files) = 8787703243
127  //data/DPMJmm_B_EMC7 ( 235/ 235 files) = 8748201639
128  //data/DPMJmm_B_EMC8 ( 232/ 232 files) = 8639730623
129 
130  std::map<long int,long int> evcnts =
131  { {1100,24342651889}, {1101,23723523588}, {1102,21921465124}, {1103,21192161015}, {1104,20899111775}, {1105,21264224163}, {1106,21918977223}, {1107,20327238649}, // Setup A, pbp -> J/psi (-> e+ e-) pi+ pi-
132  {1120,21266925718}, {1121,21721036293}, {1122,21144322234}, {1123,20364376890}, {1124,20899648526}, {1125,21667979271}, {1126,22088200372}, {1127,21308984620}, // Setup B, pbp -> J/psi (-> e+ e-) pi+ pi- //FIXME modes 1120 ...
133 
134  {1200,8821819045}, {1201,9010185209}, {1202,7969422935}, {1203,8081829652}, {1204,7707163073}, {1205,7482700194}, {1206,7781565342}, {1207,7967077007}, // Setup A, pbp -> J/psi (-> mu+ mu-) pi+ pi-
135  {1220,8229011119}, {1221,7705315505}, {1222,8043282415}, {1223,7780855820}, {1224,8230265044}, {1225,8005672265}, {1226,7407676699}, {1227,7932457375} // Setup B, pbp -> J/psi (-> mu+ mu-) pi+ pi- //FIXME modes 1220...
136  };
137 
138 
139  //for ( auto x:evcnts) cout <<"mode "<<x.first<<" : "<<x.second<<endl;
140  //cout <<endl;
141 
142  int bkgmode = sigmode+1000;
143  //if (sigmode%100>19) bkgmode+=1;
144 
145  TFile *f = new TFile(fname);
146  TTree *t = (TTree*)f->Get("ntp1");
147 
148  TFile fana("anaJ.root","UPDATE");
149 
150  TGraphErrors *g[4];
151  g[0] = new TGraphErrors(8);
152  g[1] = new TGraphErrors(8);
153  g[2] = new TGraphErrors(8);
154  g[3] = new TGraphErrors(8);
155 
156  confgraph(g[0], "signal to noise");
157  confgraph(g[1], "significance");
158  confgraph(g[2], "signal efficiency");
159  confgraph(g[3], "background efficiency");
160 
161  TF1 *f1=new TF1("f1","0.5*[0]*(1.0-TMath::Erf((x-[1])/[2]))+[3]",0,100.);
162  f1->SetParLimits(2,15,100);
163 
164  // -------------------------------------------------------------------
165  // Supermodul | 1 | 2 | 3 | 4 | 5 | 6 | 7 |
166  // num Alveolen in theta | 1 | 3 | 3 | 3 | 3 | 3 | 2 |
167  // coverage from 22 to |140.0|133.4|113.8| 94.1| 74.4| 54.8| 35.1|
168  // -------------------------------------------------------------------
169 
170  // missing fraction in [%] of EMC = [ 1.0 - (tht_max - 22°)/118° ]* 100
171  double emc_rmv[8] = { 0., 6., 23., 39., 56., 72., 89., 100.};
172 
173  //TString lab[4] = { ";supermodules missing;S/B", ";supermodules missing;significance [#sigma]", ";supermodules missing;signal efficiency [%]", ";supermodules missing;background efficiency [%]"};
174  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]"};
175 
176 
177  TCanvas *c1 = new TCanvas("c1","c1",1000,800);
178  c1->Divide(2,2,0.0001,0.0001);
179 
180  double sig_S = 50;
181  double sig_B = 46e6;
182 
183  double Lint = 1170*2;
184  double fBR = 0.05*0.06;
185 
186  double S_dat = sig_S * fBR * Lint;
187  double B_dat = sig_B * Lint;
188 
189  cout <<"S:B = "<<S_dat/B_dat<<" mode_S="<<sigmode<<" mode_B="<<bkgmode<<endl;
190 
191  printf (" i | S0 | B0 | S | B | fS | fB | S* | B* |\n");
192  printf ("----+----------+------------+----------+---------+----------+------------+---------+---------+\n");
193  for (int i=0;i<8;++i)
194  {
195  long int S0 = 1e5;
196  long int B0 = evcnts[sigmode+i+1000];
197 
198  //cout <<"S0 = "<<S0<<" B0 = "<<B0<<endl;
199 
200  double fS = S_dat/S0;
201  double fB = B_dat/B0;
202 
203  TString sigcut = Form("%s && mode==%d", cut.Data(), sigmode+i);
204  TString bkgcut = Form("%s && mode==%d", cut.Data(), bkgmode+i);
205 
206 
207  double S = (double) cntEvt(t, sigcut);
208  double B = (double) cntEvt(t, bkgcut);
209  double dS = sqrt(S);
210  double dB = sqrt(B);
211 
212  printf("%2d | %6ld | %6.3f G | %6d | %5d | %6.4f | %7.4f | %5.0f | %5.0f |\n", i, S0, (double)B0/1e9, (int)S, (int)B, fS, fB, S*fS, B*fB);
213  //cout <<i<<": S="<<S<<" B="<<B<<" --> S*"<<fS<<"="<<S*fS<<" B*"<<fB<<"="<<B*fB<<endl;
214 
215  double SN = S*fS/(B*fB);
216  double dSN = SN*sqrt(dS*dS/(S*S) + dB*dB/(B*B));
217 
218  g[2]->SetPoint(i, emc_rmv[i], SN);
219  g[2]->SetPointError(i, 0, dSN);
220 
221 
222  double Z = S*fS/sqrt(S*fS+B*fB);
223  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));
224 
225  g[3]->SetPoint(i, emc_rmv[i], Z);
226  g[3]->SetPointError(i, 0, dZ);
227 
228 
229  double effS = S/S0;
230  double deffS = effS*sqrt( 1./S + 1./S0 );
231 
232  g[0]->SetPoint(i, emc_rmv[i], effS*100.);
233  g[0]->SetPointError(i, 0, deffS*100.);
234 
235 
236  double effB = B/B0;
237  double deffB = effB*sqrt( 1./B + 1./B0 );
238 
239  g[1]->SetPoint(i, emc_rmv[i], effB*100.);
240  g[1]->SetPointError(i, 0, deffB*100.);
241  }
242 
243  TH1F *h[4];
244 
245  gStyle->SetOptFit(0);
246  gStyle->SetOptStat(0);
247 
248  TString plotnam[4]={"effs","effb","sn","sign"};
249 
250  for (int i=0;i<4;++i)
251  {
252  c1->cd(i+1);
253  h[i] = createHistoGraph(g[i],tit+lab[i]);
254  h[i]->Draw();
255  g[i]->Draw("P same");
256 
257  f1->SetParameters(TMath::MaxElement(g[i]->GetN(), g[i]->GetY()),70, 10, TMath::MinElement(g[i]->GetN(), g[i]->GetY()));
258  if (i==1)
259  {
260  f1->SetParameters(TMath::MaxElement(g[i]->GetN(), g[i]->GetY()),120, -10, TMath::MinElement(g[i]->GetN(), g[i]->GetY()));
261  f1->SetParLimits(1,110,200);
262  }
263 
264 
265  //if (i<3)
266  g[i]->Fit("f1","q");
267  //else g[i]->Fit("f2");
268  g[i]->SetName(plotnam[i]+"_"+((TString)fname(5,5)));
269  g[i]->Write();
270  }
271 
272  fname.ReplaceAll(".root","");
273  cut.ReplaceAll("&&","_AND_");
274  cut.ReplaceAll("||","_OR_");
275  cut.ReplaceAll("!","_NOT_");
276  cut.ReplaceAll(">","_lg_");
277  cut.ReplaceAll("<","_sm_");
278  cut.ReplaceAll(".","_");
279 
280  c1->SaveAs(Form("fig/%s__%03d__%s.gif",fname.Data(), sigmode, cut.Data()));
281  c1->SaveAs(Form("fig/%s__%03d__%s.C",fname.Data(), sigmode, cut.Data()));
282 
283  fana.Close();
284 }
285 
286 
287 /*
288 Found 421 files with pattern data/DPMJee_A_EMC1 containing event info Sum of events = 24342651889
289 Found 412 files with pattern data/DPMJee_A_EMC2 containing event info Sum of events = 23723523588
290 Found 392 files with pattern data/DPMJee_A_EMC3 containing event info Sum of events = 21921465124
291 Found 382 files with pattern data/DPMJee_A_EMC4 containing event info Sum of events = 21192161015
292 Found 379 files with pattern data/DPMJee_A_EMC5 containing event info Sum of events = 20899111775
293 Found 382 files with pattern data/DPMJee_A_EMC6 containing event info Sum of events = 21264224163
294 Found 391 files with pattern data/DPMJee_A_EMC7 containing event info Sum of events = 21918977223
295 Found 371 files with pattern data/DPMJee_A_EMC8 containing event info Sum of events = 20327238649
296 --------------------
297 Found 383 files with pattern data/DPMJee_B_EMC1 containing event info Sum of events = 21266925718
298 Found 388 files with pattern data/DPMJee_B_EMC2 containing event info Sum of events = 21721036293
299 Found 381 files with pattern data/DPMJee_B_EMC3 containing event info Sum of events = 21144322234
300 Found 373 files with pattern data/DPMJee_B_EMC4 containing event info Sum of events = 20364376890
301 Found 379 files with pattern data/DPMJee_B_EMC5 containing event info Sum of events = 20899648526
302 Found 387 files with pattern data/DPMJee_B_EMC6 containing event info Sum of events = 21667979271
303 Found 393 files with pattern data/DPMJee_B_EMC7 containing event info Sum of events = 22088200372
304 Found 383 files with pattern data/DPMJee_B_EMC8 containing event info Sum of events = 21308984620
305 
306 
307 //Found 389 files with pattern data/DPMJee_A_EMC1 containing event info Sum of events = 21720072869
308 //Found 386 files with pattern data/DPMJee_A_EMC2 containing event info Sum of events = 21593645357
309 //Found 335 files with pattern data/DPMJee_A_EMC3 containing event info Sum of events = 17249102681
310 //Found 332 files with pattern data/DPMJee_A_EMC4 containing event info Sum of events = 17094178566
311 //Found 320 files with pattern data/DPMJee_A_EMC5 containing event info Sum of events = 16064391492
312 //Found 328 files with pattern data/DPMJee_A_EMC6 containing event info Sum of events = 16840054994
313 //Found 339 files with pattern data/DPMJee_A_EMC7 containing event info Sum of events = 17659373515
314 //Found 320 files with pattern data/DPMJee_A_EMC8 containing event info Sum of events = 16148632260
315 //--------------------
316 //Found 332 files with pattern data/DPMJee_B_EMC1 containing event info Sum of events = 17087475331
317 //Found 331 files with pattern data/DPMJee_B_EMC2 containing event info Sum of events = 17047402087
318 //Found 323 files with pattern data/DPMJee_B_EMC3 containing event info Sum of events = 16391336239
319 //Found 330 files with pattern data/DPMJee_B_EMC4 containing event info Sum of events = 16841138917
320 //Found 323 files with pattern data/DPMJee_B_EMC5 containing event info Sum of events = 16310034426
321 //Found 327 files with pattern data/DPMJee_B_EMC6 containing event info Sum of events = 16750649874
322 //Found 346 files with pattern data/DPMJee_B_EMC7 containing event info Sum of events = 18236171676
323 //Found 327 files with pattern data/DPMJee_B_EMC8 containing event info Sum of events = 16721268259
324 
342 --------------------
343 --------------------
344 Found 238 files with pattern data/DPMJmm_A_EMC1 Sum of events = 8821819045
345 Found 242 files with pattern data/DPMJmm_A_EMC2 Sum of events = 9010185209
346 Found 214 files with pattern data/DPMJmm_A_EMC3 Sum of events = 7969422935
347 Found 218 files with pattern data/DPMJmm_A_EMC4 Sum of events = 8081829652
348 Found 208 files with pattern data/DPMJmm_A_EMC5 Sum of events = 7707163073
349 Found 203 files with pattern data/DPMJmm_A_EMC6 Sum of events = 7482700194
350 Found 209 files with pattern data/DPMJmm_A_EMC7 Sum of events = 7781565342
351 Found 214 files with pattern data/DPMJmm_A_EMC8 Sum of events = 7967077007
352 --------------------
353 Found 221 files with pattern data/DPMJmm_B_EMC1 Sum of events = 8229011119
354 Found 208 files with pattern data/DPMJmm_B_EMC2 Sum of events = 7705315505
355 Found 217 files with pattern data/DPMJmm_B_EMC3 Sum of events = 8043282415
356 Found 210 files with pattern data/DPMJmm_B_EMC4 Sum of events = 7780855820
357 Found 222 files with pattern data/DPMJmm_B_EMC5 Sum of events = 8230265044
358 Found 215 files with pattern data/DPMJmm_B_EMC6 Sum of events = 8005672265
359 Found 199 files with pattern data/DPMJmm_B_EMC7 Sum of events = 7407676699
360 Found 213 files with pattern data/DPMJmm_B_EMC8 Sum of events = 7932457375
361 */
printf
printf("RealTime=%f seconds, CpuTime=%f seconds\n", rtime, ctime)
t1
Int_t t1
Definition: hist-t7.C:106
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
analyse_J
void analyse_J(TString fname="ntp1_Jee_A.root", TString cut="xd0d0pide>0.5&&xd0d1pide>0.5&&chi24c<50", int sigmode=100, TString tit="")
Definition: analyse_J.C:91
cut
double cut[MAX]
Definition: autocutx.C:36
f
TFile * f
Definition: bump_analys.C:12
confgraph
void confgraph(TGraph *g, TString tit, int col=1, int marker=20)
Definition: analyse_J.C:47
c1
c1
Definition: plot_dirc.C:35
createHistoGraph
TH1F * createHistoGraph(TGraph *g, TString tit="", double xmin=0, double xmax=0)
Definition: analyse_J.C:60
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
cntEvt
long int cntEvt(TTree *t1, TString cut, TString var="ev")
Definition: analyse_J.C:3
Generated on Sun Aug 18 2019 03:01:05 for FairRoot/PandaRoot by   doxygen 1.8.8