StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
addFills.C
1 TFile *outH;
2 const int mxR=3;
3 const int mxC=16;// was 16
4 char *cutL[mxC]={"All", "Pt1", "Pt2", "Pt3", "Pt4", "Pt5", "Pt6", "PtL", "PtM", "PtH", "EtaBB", "EtaBc", "EtaFc", "EtaFF", "Qpos", "Qneg"};
5 
6 const int mxAmp=9;
7 char *ampL[mxAmp]={"a1:P", "-b1:Q", "c0:PQ", "c2:PQ", "a0", "a2", "b0", "b2", "c1"};
8 
9 const int mxGr=200;
10 TGraphErrors * grL[mxGr];
11 int nGr=0;
12 
13 // logic of Lum & pol
14 float minPol=0.07; // was 0.07
15 float minPol2=0.01; // was 0.01
16 
17 
18 addFills( TString outPath="defaultB-H/") {
19  TString inpDir="/star/data04/sim/balewski/LcpRun2/maxEta1.0/";
20 
21  char *fillL=" F2258 F2266 F2161 ";
22  // fillL=" F2053 F2075 F2076 F2083 F2095 F2102 F2105 F2110 F2116 F2127 F2132 F2134 F2135 F2136 F2147 F2153 F2161 ";
23 
24  fillL=" F2053 F2075 F2076 F2083 F2095 F2102 F2105 F2110 F2116 F2127 F2132 F2134 F2135 F2136 F2147 F2153 F2161 F2162 F2175 F2178 F2181 F2185 F2187 F2192 F2193 F2196 F2201 F2207 F2208 F2212 F2216 F2222 F2235 F2246 F2251 F2257 F2258 F2266 F2269 F2275 F2277 F2281 F2289 F2290 F2301 F2303 F2304";
25 
26  gStyle->SetPalette(1,0);
27  gStyle->SetOptStat(1111);
28  gROOT->LoadMacro("getPol.C");
29 
30  createHistTgr("final/"+outPath+"asyVer1.hist.root");
31 
32  char *fill=strtok(fillL," "); // init 'strtok'
33  int nFill=0;
34  do {
35  //.................. access R1,2,3 histos
36  TString fname=inpDir+outPath+"r"+fill+".hist.root";
37  TFile * inpH=new TFile(fname);
38  assert(inpH->IsOpen());
39  if(!inpH->IsOpen()) {
40  printf("#fail %s does not open, skip\n",fname.Data());
41  continue;
42  }
43  //inpH->ls();
44  nFill++;
45  //printf("\nprocess fill %02d '%s' \n",nFill,fill);
46  int xFill=atoi(fill+1);
47  assert(xFill);
48  sortCoef(inpH,xFill);
49 
50  }while(fill=strtok(0," "));
51 
52  //nicePlot(); return;
53  // plotTG("Pt1",3); return;
54 
55  int icut;
56  for(icut=0;icut<mxC;icut++) {
57  char *cut=cutL[icut];
58  plotTG(cut,3,"final/"+outPath);
59  }
60 
61  // save histo+TGraph's
62  outH->cd();
63  int i;
64  for(i=0;i<nGr;i++) {// add TGrpahs only
65  if(outH->Get(grL[i]->GetName())) continue;
66  grL[i]->Write();
67  }
68 
69  // outH->ls();
70  outH->Write();
71  return;
72 }
73 
74 
75 //==========================================
76 //==========================================
77 TGraphErrors * fetchTG(TString name) {
78  //printf("fetch-->'%s'\n",name.Data());
79  TGraphErrors *gr=0;
80  int i;
81  for(i=0;i<nGr && gr==0;i++) {
82  if(!name.Contains(grL[i]->GetName())) continue;
83  gr=grL[i];
84  }
85  assert(gr);
86  return gr;
87 }
88 
89 
90 //==========================================
91 //==========================================
92 void sortCoef(TFile *inp, int xFill){
93  assert(inp->IsOpen());
94 
95  float P,Q,eP,eQ;
96  getPol(xFill,P,Q,eP,eQ);
97  float PQ, ePQ;
98  PQ=P*Q;
99  ePQ=PQ*sqrt(eP*eP/P/P +eQ*eQ/Q/Q);
100  //printf("%d P*Q=%f /- %f\n",xFill,PQ,ePQ);
101 
102  int cLum=true, cP=true, cQ=true,cPQ=true;
103  if(xFill<2189) cLum=false;
104  if(P<minPol) cP=false;
105  if(Q<minPol) cQ=false;
106  if(PQ<minPol2) cPQ=false;
107  printf(" logic cL=%d cP=%d cQ=%d cPQ=%d\n",cLum,cP,cQ,cPQ);
108 
109 
110  int icut;
111  for(icut=0;icut<mxC;icut++) {
112  // for(icut=0;icut<2;icut++) {
113  char *cut=cutL[icut];
114  // fetch histos with fitted ratios
115  TH1F * h[mxR];
116  int k;
117  for(k=0;k<mxR;k++) {
118  char name[100];
119  sprintf(name,"r%d*%s",k+1,cut);
120  //printf("nn2=%s=\n",name);
121  h[k]=(TH1F *)inp->Get(name);
122  //assert(h[k]);
123  }
124 
125  // extract fit params and add to Tgraph's
126  int ir;
127  for(ir=0;ir<mxR;ir++) {
128  if( h[ir]==0) continue;
129  TF1 *ff=h[ir]->GetFunction("fCos012");
130  assert(ff);
131 
132  const int mxPar=3; // out
133  int ip;
134  for(ip=0;ip<mxPar;ip++) {
135  TString grName;
136  float pol=1, ePol=0;
137 
138  // manipulate with fite coef depending on observable
139 
140  if(ir==0 && ip==0 ) { grName="a0";
141  } else if (ir==0 && ip==1 && cP ) { grName="a1:P";
142  pol=P;
143  ePol=eP;
144  } else if (ir==0 && ip==2 ) { grName="a2";
145  } else if (ir==1 && ip==0 ) { grName="b0";
146  } else if (ir==1 && ip==1 && cQ ) { grName="-b1:Q";
147  pol=-Q;
148  ePol=eQ;
149  } else if (ir==1 && ip==2 ) { grName="b2";
150  } else if (ir==2 && ip==0 && cPQ && cLum) { grName="c0:PQ";
151  pol=PQ;
152  ePol=ePQ;
153  } else if (ir==2 && ip==1 ) { grName="c1";
154  } else if (ir==2 && ip==2 & cPQ ) { grName="c2:PQ";
155  pol=PQ;
156  ePol=ePQ;
157  } else {
158  continue;
159  }
160 
161  // add result to the proper tgraph
162  grName+="*";
163  grName+=cut;
164  TGraphErrors * gr=fetchTG(grName);
165  // printf("ir=%d ip=%d name=%s=%p\n",ir,ip,grName.Data(),gr);
166  double val=ff->GetParameter(ip);
167  double err=ff->GetParError(ip);
168  storeVal(gr,xFill, val,err,pol,ePol);
169 
170  }
171 
172  // accumulate chi2
173  double chi2=ff->GetChisquare();
174  double ndf=ff->GetNDF();
175  //printf("ir%d %f %f \n",ir,chi2,ndf);
176  char name[100];
177  sprintf(name,"chi2R%d*%s",ir+1,cut);
178  TH1F *ho=(TH1F*)fetchTG(name);
179  assert(ho);
180  ho->Fill(chi2/ndf);
181  }// end of loop over ir
182  }// end of loop over cuts
183 
184 }
185 
186 
187 //==========================================
188 //==========================================
189 void createHistTgr(TString fname) {
190  outH=new TFile(fname,"RECREATE");
191  assert(outH->IsOpen());
192  printf("save outH -->%s\n", fname.Data());
193 
194  // create TGraphis with fit amplitudes
195 
196  char *ampTit[mxAmp]={"An, ~ cos(#phi), Blue","An, ~ cos(#phi), Yellow","A#Sigma, no #phi ","A#Delta, ~ cos(2#phi)",
197  "a0 -instrumental", "a2 -instrumental",
198  "b0 -instrumental", "b2 -instrumental",
199  "c1 -instrumental"};
200 
201  int ampCol[mxAmp]={kBlue,kYellow,kGreen,kMagenta,kBlack,
202  kBlack,kBlack,kBlack,kBlack};
203  int ic,iam=0;
204  for (ic=0;ic<mxC;ic++) {
205  for (iam=0;iam<mxAmp;iam++) {
206  char name[100];
207  sprintf(name,"%s*%s",ampL[iam],cutL[ic]);
208  //printf("ic=%d iam=%d name=%s=\n",ic,iam,name);
209  TGraphErrors *gr =new TGraphErrors;
210  gr->SetName(name);
211  gr->SetTitle(ampTit[iam]);
212  gr->SetMarkerColor(ampCol[iam]);
213  gr->SetMarkerSize(0.8);
214  gr->SetMarkerStyle(21);
215  assert(nGr<mxGr);
216  grL[nGr++]=gr;
217  }
218  int ir=0;
219  for(ir=0;ir<mxR;ir++) {
220  char name[100], tit[200];
221  sprintf(name,"chi2R%d*%s",ir+1,cutL[ic]);
222  sprintf(tit,"chi2 for R%d(#Phi), LCP cut=%s",ir+1,cutL[ic]);
223  //printf("ic=%d ir=%d name=%s=\n",ic,ir,name);
224  TH1F *h=new TH1F(name, tit,20,0,4);
225  assert(nGr<mxGr);
226  grL[nGr++]=(TGraphErrors*)h;
227  }
228  }
229 
230 }
231 
232 //==========================================
233 //==========================================
234 void plotTG(char *cut, int flag=0, TString plotDir="./") {
235  printf("plot cut -->%s\n",cut);
236 
237  c=new TCanvas(cut,cut,900,700);
238  c->Divide(4,3);
239  int ampPan[mxAmp]={2,6,9,11,1,3,5,7,10};
240  int chi2Pan[mxR]={4,8,12};
241 
242  int iam=0;
243  for (iam=0;iam<mxAmp;iam++) {
244  char name[100];
245  sprintf(name,"%s*%s",ampL[iam],cut);
246  TGraphErrors * gr=fetchTG(name);
247  int n=gr->GetN();
248  printf("\niam=%d name='%s', N=%d\n",iam,name,n);
249  // gr->Print();
250  if(n<=0) continue;
251  c->cd(ampPan[iam]);
252  gr->Draw("AP");
253  gr->Fit("pol0");
254 
255  // draw +/- sigma of the fit
256  TF1 *ff=gr->GetFunction("pol0");
257  assert(ff);
258  float val=ff->GetParameter(0);
259  float err=ff->GetParError(0);
260 
261  float chi2ndf=0;
262  if (ff->GetNDF()>0)
263  chi2ndf=ff->GetChisquare()/ff->GetNDF();
264 
265  drawFitError(gr); // draw +/- sigma of the fit
266  c->Update();
267  modifyLabels();
268 
269  float xSig=fabs(val)/err;
270  char cKey='N';
271  if(xSig>2) cKey='Y';
272 
273  printf("#ampl= %s , cut= %s , nFill= %d , <y>= %f , sig<y>= %f , xSig= %.1f , nonZero= %c\n",
274  ampL[iam], cut,gr->GetN(), val,err, xSig,cKey);
275  printf("## , %s , %s , %f , %f , %f\n",
276  ampL[iam], cut, val,err, chi2ndf);
277  }
278 
279  int ir=0;
280  for(ir=0;ir<mxR;ir++) {
281  char name[100];
282  sprintf(name,"chi2R%d*%s",ir+1,cut);
283  TH1F *ho=(TH1F*)fetchTG(name);
284  assert(ho);
285  c->cd(chi2Pan[ir]);
286  ho->Draw();
287  printf("%s nE=%d\n",ho->GetName(),ho->GetEntries());
288  c->Update();
289  modifyLabels();
290  }
291 
292  if(flag<=0) return;
293  TString outFig=plotDir+"asy";
294  outFig+=cut;
295  if(flag&1) c->Print(outFig+".ps");
296  if(flag&2) c->Print(outFig+".gif");
297 
298 
299 }
300 
301 //--------------------------------------------------------
302 //--------------------------------------------------------
303 
304 void storeVal(TGraphErrors *gr, float x, float y,float ey,float fac, float efac=0){
305  assert(gr);
306 
307  int n=gr->GetN();
308  double val=y/fac;
309 
310  double e1=ey/y;
311  double e2=efac/fac;
312  double err=fabs(val)*sqrt(e1*e1+e2*e2);
313  // printf("%f %f %f %f\n",e1,e2,ey/fabs(fac),err);
314  gr->SetPoint(n,x,val);
315  gr->SetPointError(n,.0,err);
316 
317  // printf("add: %s %f=fac y=%f ey=%f val=%f err=%f\n",gr->GetName(),fac,y,ey,val,err);
318 }
319 
320 
321 //==========================================
322 //==========================================
323 void nicePlot(){
324 
325  c=new TCanvas();
326 
327  TGraphErrors * gr=fetchTG("a1:P*All");
328  //TGraphErrors * gr=fetchTG("c2:PQ*All");
329  int n=gr->GetN();
330  printf("\n name='%s', N=%d\n",gr->GetName(),n);
331  assert(n>0);
332  // gr->Print();
333  gr->Draw("AP");
334  gr->Fit("pol0");
335 
336  drawFitError(gr); // draw +/- sigma of the fit
337 
338  return;
339  c->Update();
340  TPaveStats *st1 =( TPaveStats *)gPad->GetPrimitive("stats");
341  st1->SetX1NDC(0.4);
342  st1->SetX2NDC(1.);
343  st1->SetY1NDC(.7);
344  st1 =( TPaveStats *)gPad->GetPrimitive("title");
345  st1->SetX2NDC(.4);
346  st1->SetY1NDC(.8);
347  c->Update();
348 
349 }
350 
351 
352 //==========================================
353 //==========================================
354 void drawFitError( TGraphErrors * gr) {
355  int col=kRed;
356  TF1 *ff=gr->GetFunction("pol0");
357  assert(ff);
358 
359  ff->SetLineColor(col);
360  ff->SetParName(0,"avr");
361 
362  float val=ff->GetParameter(0);
363  float err=ff->GetParError(0);
364  float upY=val+err, dwY=val-err;
365 
366  TAxis *ax=gr->GetXaxis();
367  float x1=ax->GetXmin();
368  float x2=ax->GetXmax();
369  ax->SetTitle("RHIC FILL \#");
370  //printf("x1=%f x2=%f\n",x1,x2);
371 
372  TLine *ln0=new TLine(x1,0.,x2,0.);
373  x1+=2;
374  x2-=2;
375  TLine *upL=new TLine(x1,upY,x2,upY);
376  TLine *dwL=new TLine(x1,dwY,x2,dwY);
377  upL->SetLineStyle(3);
378  dwL->SetLineStyle(3);
379 
380  upL->SetLineColor(col);
381  dwL->SetLineColor(col);
382 
383  upL->Draw();
384  dwL->Draw();
385  ln0->Draw();
386 
387  float xSig=fabs(val)/err;
388  char cKey='N';
389  if(xSig>2) cKey='Y';
390  printf(" nData= %d , <y>= %f , sig<y>= %f , xSig= %.1f , nonZero= %c\n",
391  gr->GetN(), val,err, xSig,cKey);
392 
393 }
394 
395 
396 
397 //==========================================
398 //==========================================
399 void modifyLabels(){
400  TPaveStats *st1 =( TPaveStats *)gPad->GetPrimitive("stats");
401  st1->SetX1NDC(0.3);
402  st1->SetX2NDC(1.);
403  st1->SetY1NDC(.8);
404  return;
405  st1 =( TPaveStats *)gPad->GetPrimitive("title");
406  st1->SetX1NDC(.6);
407  st1->SetX2NDC(1.0);
408  st1->SetY1NDC(.1);
409  st1->SetY1NDC(.3);
410 
411 
412 }
413 
414 
415 
Generated by   doxygen 1.8.5