FairRoot/PandaRoot
Functions
analyse_etac1.C File Reference

Go to the source code of this file.

Functions

long int cntEvt (TTree *t1, TString cut, TString var="ev")
 
void confgraph (TGraph *g, TString tit, int col=1, int marker=20)
 
TH1F * createHistoGraph (TGraph *g, TString tit="", double xmin=0, double xmax=0)
 
void analyse_etac1 (TString fname="ntp5_etac1_A.root", TString cut="chi24c<50", int sigmode=300, TString tit="")
 

Function Documentation

void analyse_etac1 ( TString  fname = "ntp5_etac1_A.root",
TString  cut = "chi24c<50",
int  sigmode = 300,
TString  tit = "" 
)

Definition at line 91 of file analyse_etac1.C.

References c1, cntEvt(), confgraph(), createHistoGraph(), cut, f, f1, g, h, i, printf(), sqrt(), t, TString, and Z.

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 }
printf("RealTime=%f seconds, CpuTime=%f seconds\n", rtime, ctime)
Int_t i
Definition: run_full.C:25
TF1 * f1
Definition: reco_analys2.C:50
friend F32vec4 sqrt(const F32vec4 &a)
Definition: P4_F32vec4.h:29
long int cntEvt(TTree *t1, TString cut, TString var="ev")
Definition: analyse_etac1.C:3
TFile * g
TH1F * createHistoGraph(TGraph *g, TString tit="", double xmin=0, double xmax=0)
Definition: analyse_etac1.C:60
double cut[MAX]
Definition: autocutx.C:36
void confgraph(TGraph *g, TString tit, int col=1, int marker=20)
Definition: analyse_etac1.C:47
TFile * f
Definition: bump_analys.C:12
c1
Definition: plot_dirc.C:35
double Z
Definition: anaLmdDigi.C:68
TTree * t
Definition: bump_analys.C:13
long int cntEvt ( TTree *  t1,
TString  cut,
TString  var = "ev" 
)

Definition at line 3 of file analyse_etac1.C.

References cnt, i, and TString.

Referenced by analyse_etac1(), analyse_etac1_slc(), analyse_J(), analyse_J_slc(), and analyse_phi_base().

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 }
Int_t i
Definition: run_full.C:25
Int_t t1
Definition: hist-t7.C:106
double cut[MAX]
Definition: autocutx.C:36
Int_t cnt
Definition: hist-t7.C:106
void confgraph ( TGraph *  g,
TString  tit,
int  col = 1,
int  marker = 20 
)

Definition at line 47 of file analyse_etac1.C.

References col.

Referenced by analyse_etac1(), analyse_etac1_slc(), analyse_J(), analyse_J_slc(), analyse_phi(), analyse_phi2(), analyse_phi_base(), combinePlotsEtac(), combinePlotsEtac_slc(), combinePlotsJ(), combinePlotsJ_slc(), and combinePlotsPhi().

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 }
int col
Definition: anaLmdDigi.C:67
TFile * g
TH1F* createHistoGraph ( TGraph *  g,
TString  tit = "",
double  xmin = 0,
double  xmax = 0 
)

Definition at line 60 of file analyse_etac1.C.

References cnt, h, xmax, and xmin.

Referenced by analyse_etac1(), analyse_etac1_slc(), analyse_J(), analyse_J_slc(), analyse_phi(), analyse_phi2(), analyse_phi_base(), combinePlotsEtac(), combinePlotsEtac_slc(), combinePlotsJ(), combinePlotsJ_slc(), and combinePlotsPhi().

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 }
Double_t xmax
TFile * g
Int_t cnt
Definition: hist-t7.C:106
Double_t xmin