StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
rdN2AL2012.C
1 /* converts raw inputs yields in to physics asymetries vs. polBeam-eta
2  version 2.2
3 
4 INPUT:
5  - BGmodel(year) as root hists
6  - Wyield+lumis(year) as root histo (STAR data or toy MC)
7  - polarization(year) as input params
8 
9 Computation is done in 2 steps:
10 1) in loop over starPhys-eta bins
11 AL,ALL,EpsNull are computed for both beams * background and stored in tempHistos bins [1-8]
12 
13 2) SSA, DSA from tempHistos are combined according to polBeam-eta index [10-20] and store in final histos
14 
15 OUTPUT : final histos are saved vs. polBeam eta float.
16 
17 Now WB & WE are used with the following hardcoded rules:
18 bins 1-4 : only WB
19 bin 5 : WB+WE
20 bin 6,7 : only WE
21 bin 8 : only WB
22 */
23 
24 /* mapping of spin4-index to helicity at STAR */
25 enum { ka=10, /* STAR pol B+ Y + */
26  kb=9, /* STAR pol B+ Y - */
27  kc=6, /* STAR pol B- Y + */
28  kd=5, /* STAR pol B- Y - */
29  mxSS=4,
30  mxEta=8, /* # of used STAR eta bins */
31  mxQ=2 /*0=W+, 1=W- */
32 };
33 
34 int kA[mxSS]={ka,kb,kc,kd};
35 int minYieldPerSpin=1; // cut-off to drop whole eta bin due to not reliable error propagation
36 
37 // working histos to store AL, ALL from single beams vs. starEtaBin
38 TH1F * hAL1[mxQ], *hAL2[mxQ], *hALL[mxQ], *hEpsNull[mxQ], *hNused[mxQ], *hBeta[mxQ];
39 // final histos to store AL from both beams vs. polBeamEtaBin
40 TH1F * hALfin[mxQ],* hLumi;
41 double SumYieldPerSpin[mxQ][mxEta][mxSS]; // only for nice printing at the end, not saved
42 
43 TFile *fdWAL=0;
44 TFile *fdBG=0;
45 TFile *fdData=0;
46 int addEndcap=1; // include WE events in AL computation
47 int skipWminus=0; // used '1' for testing
48 
49 //=================================
50 void rdN2AL2012(TString inpCore="run12toySetB", TString bgYear="2012", double pol1=0.5,double pol2=0.99 ) {
51  // inpCore="8.15.12/run9long";pol1=0.376; pol2=0.395 ;
52  //inpCore="../8.24.12/run12long";pol1=0.55; pol2=0.57 ;
53  inpCore="../results_August28/8.29.12/run12long";pol1=0.549; pol2=0.571;
54  TString bgName="8.15.12/bkgdBetaModel_"+bgYear+".hist.root";
55 
56  // fdWAL=new TFile("WALModel.hist.root"); assert(fdWAL->IsOpen()); // activate to plot on final page
57 
58  gROOT->LoadMacro("setWALmodel.C");
59 
60  fdBG=new TFile(bgName); assert(fdBG->IsOpen());
61  //fdBG->ls();
62 
63  TString dataName=inpCore+".wana.hist.root";
64  fdData=new TFile(dataName); assert(fdData->IsOpen());
65  printf("Opened %s\n",fdData->GetName());
66  //fdData->ls();
67 
68  memset(SumYieldPerSpin,0,sizeof(SumYieldPerSpin)); // just in case
69  //......... access & store lumis (works: )
70  TH1F *hAbsLum=(TH1F *)fdData->Get("Eta8/AEta8spinY1"); assert(hAbsLum); //hAbsLum->Draw();
71  hLumi=newSpin4Histo("relLumi","renormalized relative lums");// just monitor
72  double sum=0;
73  for(int is=0; is<mxSS; is++){
74  int bin=kA[is]+1;
75  double y= hAbsLum->GetBinContent(bin);
76  sum+=y;
77  hLumi->SetBinContent(bin,y);
78  //printf("y=%.0f\n",y);
79  }
80  hLumi->Scale(mxSS/sum);
81  for(int is=0; is<mxSS; is++)
82  hLumi->SetBinError(kA[is]+1,1./2/sqrt(sum));
83 
84  double rLumi[mxSS],sLumi[mxSS];
85  for(int is=0; is<mxSS; is++) {
86  rLumi[is]= hLumi->GetBinContent(kA[is]+1);
87  sLumi[is]= hLumi->GetBinError(kA[is]+1);
88  }
89  //hLumi->Fit("pol0");
90  printf("renorm lumis: %.4f %.4f %.4f %.4f QCD eve count=%.1e\n", rLumi[0],rLumi[1],rLumi[2],rLumi[3],sum);
91  printf(" have errors: %.4f %.4f %.4f %.4f 1/2/sqrt(N)=%.1e\n", sLumi[0],sLumi[1],sLumi[2],sLumi[3],1./2/sqrt(sum));
92 
93 
94  TH1F *h; TGraphErrors *gr;
95 
96  for( int iq=0; iq<mxQ; iq++) {
97  if(skipWminus && iq==1) continue; // tmp work only with + charge
98  TString Q="P"; if (iq) Q="N";
99 
100  hNused[iq]=newEtaBinHisto("Nused_"+Q, " Lum-norm. events, Q="+Q,"EVENTS");
101 
102  hAL1[iq]=h=newEtaBinHisto("AL1_"+Q, "beam 1 AL (starEta) , Q="+Q, "A_L");
103  h->SetMarkerStyle(22+iq*4); // tringleUp 22:26
104  h->SetMarkerSize(1.3);
105 
106  hAL2[iq]=h=newEtaBinHisto("AL2_"+Q, " beam 2 AL (starEta) , Q="+Q, "A_L");
107  h->SetMarkerStyle(21+iq*4); // square 21:25
108 
109  hALL[iq]=h=newEtaBinHisto("ALL_"+Q, " beam 1*2 ALL, Q="+Q, "A_LL");
110  h->SetMarkerStyle(23+iq*9); // tringleDown 23:32
111  h->SetMarkerSize(1.3);
112 
113  hEpsNull[iq]=h=newEtaBinHisto("EpsNull_"+Q, "EpsNull , Q="+Q, "Eps_null");
114  h->SetMarkerStyle(21+iq*4); // square 21:25
115 
116  hALfin[iq]=h=newEtaBinHisto("finAL_"+Q, "AL (polBeamEta) , Q="+Q, "A_L");
117  h->SetMarkerStyle(21+iq*4); // square 21:25
118 
119  //....... access background
120  // TH1F *hAlpha=(TH1F *)fdBG->Get("alpha"+Q+"_"+bgYear); assert(hAlpha);
121  hBeta[iq]=(TH1F *)fdBG->Get("beta"+Q+"_"+bgYear); assert(hBeta);
122 
123 
124  printf("\n\n ========== PHASE_I =========== Q=%s\n",Q.Data());
125  // compute AL, ALL for each beam separately
126 
127  //.... ADD loop over starPhysEta bins here ......
128  for( int starPhysEtaBin=1; starPhysEtaBin<=mxEta; starPhysEtaBin++) {
129  printf("\n------------ processing Q=%s ------- starEtaBin=%d ------\n",Q.Data(),starPhysEtaBin);
130  // if(starPhysEtaBin!=8) continue;
131  //..... access spin sorted yields .....
132  TH1F *hData=(TH1F *)fdData->Get(Form("Eta%d/AEta%dspinY2_",starPhysEtaBin,starPhysEtaBin)+Q); assert(hData);
133  TH1F *hDataE=(TH1F *)fdData->Get(Form("Eta%d/AEta%dspinEY2_",starPhysEtaBin,starPhysEtaBin)+Q); assert(hDataE);
134 
135  //hData->Draw();
136  if(addEndcap && ( starPhysEtaBin==5 )) hData->Add(hDataE);
137  if( starPhysEtaBin==7 || starPhysEtaBin==6) {
138  if (addEndcap ) hData=hDataE;
139  else hData->Reset();
140  }
141  // eta histo : AEta7spinELepEta
142  // Endcap yiled : AEta7spinEY2_N
143 
144  double beta =hBeta[iq] ->GetBinContent(starPhysEtaBin);
145  double betaErr=hBeta[iq] ->GetBinError(starPhysEtaBin);
146 
147  //tmp-start
148  if(beta<0.1) {printf("changed beta WARN\n"); beta=0.9;}
149  double alpha1=0, alpha2=0;
150  double alphaErr1=0, alphaErr2=0;
151  //tmp-end
152 
153  double M[mxSS]; //yields corrected for lumi
154  double VM[mxSS]; // variances corrected for lumi
155  //...... access Walgo yields
156  printf("raw Wyields iQ=%d starEtaBin%d: ",iq,starPhysEtaBin);
157  double sum=0;
158  int skipBin=0;
159  for(int is=0; is<mxSS; is++) {
160  double val= hData->GetBinContent(kA[is]+1); // raw yield , events
161  if (val<minYieldPerSpin) skipBin++;
162  M[is]=val/ rLumi[is];
163  VM[is]=val/ rLumi[is]/ rLumi[is];
164  sum+=val;
165  SumYieldPerSpin[iq][starPhysEtaBin-1][is]=M[is];
166  //printf(" M[%d]=%.1f V=%.1f; ",is,M[is], VM[is]);
167  printf(" %.0f ",val);
168  }
169 
170  printf(" sumN=%.1f skipBin=%d\n",sum,skipBin);
171  if(skipBin) { printf(" SKIP this starEtaBin due to low stats\n"); continue;}
172  printf("alpha1: %.4f +/- %.4f, ",alpha1, alphaErr1);
173  printf("alpha2: %.4f +/- %.4f, ",alpha2, alphaErr2);
174  printf("beta: %.4f +/- %.4f\n",beta,betaErr);
175 
176  hNused[iq]->SetBinContent(starPhysEtaBin,sum);
177 
178  double A,sA; // physics asymmetry + error
179  double eps, sEps; //raw asymmetry + error
180 
181  if(starPhysEtaBin!=8) {// do sth else for full-barrel
182  printf("..................... AL for beam 1 \n");
183  computeAL(pol1, alpha1,alphaErr1, beta,betaErr,
184  M[0]+ M[1], M[2]+ M[3],
185  VM[0]+VM[1], VM[2]+VM[3], A, sA);
186  hAL1[iq]->SetBinContent(starPhysEtaBin,A);
187  hAL1[iq]->SetBinError (starPhysEtaBin,sA);
188 
189  printf("..................... AL for beam 2 \n");
190  computeAL(pol2, alpha2,alphaErr2, beta,betaErr,
191  M[0]+ M[2], M[1]+ M[3],
192  VM[0]+VM[2], VM[1]+VM[3], A, sA);
193  hAL2[iq]->SetBinContent(starPhysEtaBin,A);
194  hAL2[iq]->SetBinError (starPhysEtaBin,sA);
195  } else { // special case for full barrel
196  printf("..................... AL for beam 1+2 over Barrel \n");
197  doEps_II( M[0], M[3], M[1]+ M[2],
198  VM[0], VM[3], VM[1]+VM[2], eps, sEps);
199  printf("epsALsym= %.3f $\\pm$ %.3f nSig=%.2f \n", eps, sEps, eps/ sEps);
200  // use alpha from beam 1 only, should be identical for this case
201  doPolBckgCorr(eps, sEps,(pol1+pol2)/2., alpha1, alphaErr1, beta, betaErr, A, sA);
202  printf("ALsym=%.3f +/- %.3f nSig=%.2f\n", A, sA, A/sA);
203  hAL1[iq]->SetBinContent(starPhysEtaBin,A);
204  hAL1[iq]->SetBinError (starPhysEtaBin,sA);
205  }
206 
207  printf("..................... ALL for beam 1*2 \n");
208  // (sic!) - the same formula is used for AL() and ALL()
209  computeAL(pol1*pol2, 0. , 0., beta,betaErr,
210  M[0]+ M[3], M[1]+ M[2],
211  VM[0]+VM[3], VM[1]+VM[2], A, sA);
212 
213  hALL[iq]->SetBinContent(starPhysEtaBin,A);
214  hALL[iq]->SetBinError (starPhysEtaBin,sA);
215 
216 
217  printf("..................... EpsNULL \n");
218  doEps_II( M[1], M[2], M[0]+ M[3],
219  VM[1], VM[2], VM[0]+VM[3], eps, sEps);
220  printf("epsNull_not= %.3f $\\pm$ %.3f nSig=%.2f \n", eps, sEps, eps/ sEps);
221  hEpsNull[iq]->SetBinContent(starPhysEtaBin,eps);
222  hEpsNull[iq]->SetBinError (starPhysEtaBin,sEps);
223 
224  } // ----- end of starEtaBins loop
225 
226  // now SSA, DSA are computed per beam , merging of results is done below
227  // combine AL, ALL from 2 beams in to single value
228 
229  printf("\n\n ========== PHASE_II =========== Q=%s\n",Q.Data());
230 
231  // work on polBeam-eta slices
232  int k1,k2;
233  for (int k=10;k <=20; k++) { // index polBeamEtaBin
234  printf("\n........ processing iQ=%d ..... polBeamEtaBin=%d ....\n",iq,k);
235  k1=k-12; // for beam 1
236  k2=17-k; // for beam 2
237  if(k==20) { k1=8; k2=-999; } // special case: barrel Ws
238 
239  double N1=0, N2=0;
240  if(k1>0) N1 = hNused[iq]->GetBinContent(k1);
241  if(k2>0) N2 = hNused[iq]->GetBinContent(k2);
242  hNused[iq]->SetBinContent(k,N1+N2);
243  if(N1+N2<=0) continue; // skip if neither beam have any contribution
244 
245  printf("..................... AL from both beams \n");
246  doOneBinAverage(k1,k2, hAL1[iq], hAL2[iq],k, hALfin[iq]);
247  if(k>=15) {
248  printf("..................... ALL from both beams \n");
249  doOneBinAverage(k1,k2, hALL[iq], hALL[iq],k, hALL[iq]);
250  printf("..................... EpsNull from both beams \n");
251  doOneBinAverage(k1,k2, hEpsNull[iq], hEpsNull[iq],k, hEpsNull[iq]);
252  }
253  } // end of loop over polBeamEtaBins
254 
255  // hALfin[iq]->Fit("pol0","","R",10.5,18.5);
256  }// ---- end of 2nd + - charge loop
257 
258  // return;
259  TString outName=inpCore+".wasy.hist.root";
260  TFile *fdOut=new TFile(outName,"recreate"); assert(fdOut->IsOpen());
261  //-------- final plots and printouts
262  for( int iq=0; iq<mxQ; iq++) {
263  if(skipWminus && iq==1) continue; // tmp work only with + charge
264 
265  makeNiceSummaryPlotPerBeam(inpCore, iq);
266  makeNiceSummaryPlotSumBeam(inpCore, iq);
267  printNiceTable(inpCore, iq, pol1, pol2);
268  hNused[iq]->Write();
269  hEpsNull[iq]->Write();
270  hALL[iq]->Write();
271  hBeta[iq]->Write();
272  hALfin[iq]->Write();
273  saveEtaSpectra(iq);
274 
275  }
276  hLumi->Write();
277  hAbsLum->Write();
278 
279 }
280 
281 
282 
283 
284 //-----------------------------------
285 void doOneBinAverage( int k1, int k2, TH1F *h1, TH1F *h2,int k, TH1F *hfin) {
286  double A1=0, A2=0, sA1=0, sA2=0, N1=0, N2=0;
287  double A1=0, A2=0, sA1=0, sA2=0;
288  if(k1>0) {
289  A1 = h1->GetBinContent(k1);
290  sA1= h1->GetBinError(k1);
291  }
292  if(k2>0) {
293  A2 = h2->GetBinContent(k2);
294  sA2= h2->GetBinError(k2);
295  }
296  double A=999,sA=0.01; // asymmetry + error
297  averageTwoValues(A1,sA1, A2,sA2, A,sA); // if one is missing it gets 0-weight = OK
298  printf(" k=%d k1=%d k2=%d -> A=%.3f +/-%.3f\n",k,k1,k2,A,sA);
299  hfin->SetBinContent(k,A);
300  hfin->SetBinError(k,sA);
301 }
302 
303 
304 //-----------------------------------
305 void averageTwoValues(double x, double sx,double y, double sy , double &A, double &sA) {
306  printf("X=%.3f sX=%.3f; Y=%.3f sY=%.3f\n",x,sx,y,sy);
307  if( sx<=0.) { A=y; sA=sy; return; }
308  if( sy<=0.) { A=x; sA=sx; return; }
309  double wx=1./sx/sx;
310  double wy=1./sy/sy;
311  double wsum=wx+wy;
312  wx/=wsum; wy/=wsum;
313  A=wx*x+wy*y;
314  sA=sqrt(wx*wx*sx*sx +wy*wy*sy*sy);
315  printf("wx=%f wy=%f --> A=%f\n",wx,wy,A);
316 }
317 
318 
319 
320 //-----------------------------------
321 void computeAL(double pol, double alpha,double alphaErr, double beta,double betaErr,
322  double Ma, double Mb,double VMa, double VMb, double &A, double &sA) {
323  printf("computeAL: pol=%.3f Ma+Mb=%.1f\n",pol, Ma+Mb);
324  printf("alpha: %.4f +/- %.4f , ",alpha, alphaErr);
325  printf("beta: %.4f +/- %.4f\n",beta,betaErr);
326  printf("Ma=%.1f Va=%.1f; ",Ma, VMa); printf(" Mb=%.1f Vb=%.1f\n",Mb, VMb);
327 
328  double eps, sEps;
329  doEps_I( Ma, Mb, VMa, VMb , eps, sEps );
330  printf("eps= %.3f $\\pm$ %.3f nSig=%.2f sig(eps)*sqrt(Ma+Mb)=%.3f\n", eps, sEps, eps/ sEps , sEps*sqrt(Ma+Mb));
331 
332  doPolBckgCorr(eps, sEps, pol, alpha, alphaErr, beta, betaErr, A, sA);
333  printf("A=%.3f +/- %.3f nSig=%.2f\n", A, sA, A/sA);
334 
335 }
336 
337 //-----------------------------
338 void doEps_I( double a, double b, double va, double vb,
339  double &eps, double &sEps) {
340  printf("doEps_I: Ma=%f Mb=%f\n",a,b);
341  double sum=a+b;
342  eps= (a-b)/sum;
343  double xx=b*b*va+ a*a*vb;
344  sEps= sqrt( 4 * xx/sum/sum/sum/sum);
345 }
346 
347 
348 
349 //------------
350 void doEps_II( double a, double b,double c, double Va, double Vb, double Vc,
351  double &eps, double &sEps) {
352  printf("doEps_II: Ma=%f Mb=%f Mc=%f\n",a,b,c);
353  // printf("doEps_II: Va=%f Vb=%f Vc=%f\n",Va,Vb,Vc);
354  double sum=a+b+c;
355  eps= (a-b)/sum;
356  double vEps=(4*a*c*Vb + c*c*(Va + Vb) + b*b*(4*Va + Vc) +
357  a*a*(4*Vb + Vc) + b*(4*c*Va - 2*a*Vc))/sum/sum/sum/sum;
358  sEps= sqrt(vEps);
359 }
360 
361 
362 
363 //------------------------------
364 void doPolBckgCorr( double eps, double sEps,double P,
365  double a, double sa, double b, double sb,
366  double &AL, double &sAL) {
367  // printf("doPolBckgCorr eps=%f P=%f a=%f b=%f\n",eps,P,a,b);
368  AL=(eps/P-a)/b;
369  double v1=sEps*sEps/P/P;
370  double v2=sa*sa;
371  double v3=AL*AL*sb*sb;
372  sAL=sqrt(v1+v2+v3)/b;
373  //printf("v1=%f v2=%f v3=%f sEps/P=%f sqrt(v1)=%f sAL=%f\n",v1,v2,v3,sEps/P, sqrt(v1),sAL);
374 }
375 
376 //------------------------
377 void splitPadX(float x, TPad **cL, TPad **cR) {
378  (*cL) = new TPad("padL", "apdL",0.0,0.,x,0.95);
379  (*cL)->Draw();
380  (*cR) = new TPad("padL", "apdL",x+0.005,0.,1.0,0.95);
381  (*cR)->Draw();
382 }
383 
384 //------------------------
385 void splitPadY(float y, TPad **cU, TPad **cD) {
386  (*cU) = new TPad("padD", "apdD",0,y+0.005,1.0,1.);
387  (*cU)->Draw();
388  (*cD) = new TPad("padU", "apdU",0.0,0.,1.,y);
389  (*cD)->Draw();
390 
391  /* use case:
392  TPad *cU,*cD; splitPadY(0.4,&cU,&cD); cU->cd(); h->Draw()
393  */
394 }
395 
396 
397 //------------------------
398 TPad *makeTitle(TCanvas *c,TString core) {
399  c->Range(0,0,1,1);
400  TPad *pad0 = new TPad("pad0", "apd0",0.0,0.95,1.,1.);
401  pad0->Draw();
402  pad0->cd();
403 
404  TPaveText *pt = new TPaveText(0,0.,1,1,"br");
405  pt->Draw();
406  TDatime dt;
407  TString txt2=core;
408 
409  txt2+=", ";
410  txt2+=dt.AsString();
411  pt->AddText(txt2);
412  txt2="--";
413  pt->AddText(txt2);
414 
415  c->cd();
416  pad = new TPad("pad1", "apd1",0.0,0.0,1,.95);
417  pad->Draw();
418  return pad;
419 }
420 
421 
422 
423 //======================================
424 //======================================
425 //======================================
426 void makeNiceSummaryPlotPerBeam(TString inpCore, int iq) {
427  TString Q="P"; if (iq) Q="N";
428  //........ nice plot with summary
429  TString titC1=inpCore+"one Q="+Q;
430  can=new TCanvas(titC1, titC1,600,720);
431  gStyle->SetOptStat(0);
432  TPad *c=makeTitle(can,inpCore+" Q="+Q+", perBeam");
433  // c->SetFillColor(kWhite);
434  c->Divide(2,2);
435  ln=new TLine(0,0,10,0); ln->SetLineStyle(1);ln->SetLineColor(kBlue);
436 
437  float yMx=1.05, xMx=8.4;
438  bxEE=new TBox(6.5,-yMx,xMx,yMx);
439  bxEE->SetFillColor(15); bxEE->SetFillStyle(3944);
440 
441 
442  TH1F *hX[4]= {hAL1[iq],hALL[iq],hAL2[iq],hEpsNull[iq]};
443  for(int ih=0;ih<4;ih++) {
444  if(ih==3) yMx/=3.;
445  c->cd(1+ih);
446  hX[ih]->Draw(); hX[ih]->SetAxisRange(0,xMx);
447  hX[ih]->SetMinimum(-yMx); hX[ih]->SetMaximum(yMx);
448  TAxis *ax=hX[ih]->GetXaxis();
449  ax->SetTitle("STAR #eta bins "); ax->SetTitleSize(0.06); ax->SetTitleOffset(0.8);
450  ax->SetLabelSize(0.05);
451  ln->Draw();
452  bxEE->Draw();
453 
454  if(ih!=2) {
455  ar=new TArrow(0.6,-yMx*.9,6.3,-yMx*.9, 0.025); ar->Draw(); ar->SetLineColor(8);
456  tx1=new TLatex(0.8,-yMx*.8,"-0.9"); tx1->Draw(); tx1->SetTextColor(8);
457  tx1=new TLatex(5.5,-yMx*.8,"+1.3"); tx1->Draw();tx1->SetTextColor(8);
458  tx1=new TLatex(8.1,yMx*.8,"Barrel"); tx1->Draw();tx1->SetTextAngle(90);
459  tx1=new TLatex(7.1,yMx*.74,"+Endcap"); tx1->Draw();tx1->SetTextAngle(90);
460  } else {
461  ar=new TArrow(0.6,-yMx*.9,6.3,-yMx*.9, 0.025,"<"); ar->Draw(); ar->SetLineColor(8);
462  tx1=new TLatex(0.8,-yMx*.8,"+0.9"); tx1->Draw();tx1->SetTextColor(8);
463  tx1=new TLatex(5.5,-yMx*.8,"-1.3"); tx1->Draw();tx1->SetTextColor(8);
464  tx1=new TLatex(7.1,yMx*.75,"-Endcap"); tx1->Draw();tx1->SetTextAngle(90);
465  }
466 
467  tx1=new TLatex(3.1,-yMx*.8,"polBeam #eta"); tx1->Draw();tx1->SetTextColor(8);
468 
469  gPad->SetGridy();
470  }
471 }
472 
473 
474 
475 
476 //======================================
477 //======================================
478 //======================================
479 void makeNiceSummaryPlotSumBeam(TString inpCore, int iq) {
480 
481  TString Q="P"; if (iq) Q="N";
482  ln=new TLine(9,0,21,0); ln->SetLineStyle(1);ln->SetLineColor(kBlue);
483  //........ nice plot with summary
484  TString titC1=inpCore+"sum Q="+Q;
485  can=new TCanvas(titC1, titC1,700,700);
486  gStyle->SetOptStat(0);
487  TPad *c=makeTitle(can,inpCore+" Q="+Q+", sumBeam");
488  c->cd();
489  TPad *cL,*cR; splitPadX(0.5,&cL,&cR);
490 
491  //------- AL plot -------
492  TH1F *h=hALfin[iq];
493  TAxis *ax=h->GetXaxis();
494  ax->SetTitle("polBeam #eta bins "); ax->SetTitleSize(0.06); ax->SetTitleOffset(0.73);
495  ax->SetLabelSize(0.05);
496 
497  cL->cd();
498  float yMx=1.05, x1=9.5, x2=20.4;
499  bxEE=new TBox(18.5,-yMx,x2,yMx); bxEE->SetFillColor(15); bxEE->SetFillStyle(3944);
500  bxEF=new TBox(x1,-yMx,10.5,yMx); bxEF->SetFillColor(15); bxEF->SetFillStyle(3944);
501 
502  h->Draw(); h->SetAxisRange(x1,x2);
503  h->SetMinimum(-yMx); h->SetMaximum(yMx);
504 
505  bxEE->Draw(); bxEF->Draw(); ln->Draw();
506  ar=new TArrow(10.6,-yMx*.85,18.3,-yMx*.85, 0.025); ar->Draw(); ar->SetLineColor(8);
507  tx1=new TLatex(10.8,-yMx*.8,"-1.3"); tx1->Draw(); tx1->SetTextColor(8);
508  tx1=new TLatex(17.1,-yMx*.8,"+1.3"); tx1->Draw();tx1->SetTextColor(8);
509  tx1=new TLatex(12.5,-yMx*.8,"polBeam #eta"); tx1->Draw();tx1->SetTextColor(8);
510  tx1=new TLatex(20.1,yMx*.8,"Barrel"); tx1->Draw();tx1->SetTextAngle(90);
511  tx1=new TLatex(19.1,yMx*.74,"+Endcap"); tx1->Draw();tx1->SetTextAngle(90);
512  tx1=new TLatex(10.1,yMx*.75,"-Endcap"); tx1->Draw();tx1->SetTextAngle(90);
513 
514  gPad->SetGridy();
515 
516  if(fdWAL) {// overlay model
517  TH1F *hWAL=fdWAL->Get("modW_AL_"+Q); assert(hWAL);
518  hWAL->Draw("same");
519  }
520 
521  cR->cd(); cR->Divide(1,2);
522 
523  //------- ALL, Anull plot -------
524  cR->cd(1);
525  TPad *cRL,*cRR; splitPadX(0.5,&cRL,&cRR);
526  x1=14.5;
527  for(int jj=0;jj<2;jj++) {
528  h=(TH1F*)hALL[iq]->Clone();
529  cRL->cd();
530  if(jj==1) { // changes for A_NULL
531  h=(TH1F*)hEpsNull[iq]->Clone();
532  cRR->cd();
533  yMx*=0.3;
534  }
535  gPad->SetLeftMargin(0.20);
536  ax=h->GetXaxis();
537  float txSz=0.08;
538  ax->SetTitle("polBeam #eta bins ");
539  ax->SetTitleSize(txSz); ax->SetTitleOffset(0.65);
540  ax->SetLabelSize(txSz); ax->SetLabelOffset(-0.01);
541  ax=h->GetYaxis();
542  ax->SetLabelSize(txSz); //ax->SetLabelOffset(-0.01);
543  ax->SetTitleSize(txSz); ax->SetTitleOffset(0.99);
544  h->Draw(); h->SetAxisRange(x1,x2);
545  // h->SetNdivisions(6,"y");
546  h->SetMinimum(-yMx); h->SetMaximum(yMx);
547  ln->Draw();
548 
549  bxEE->Draw();
550 
551  ar=new TArrow(14.6,yMx*1.05,18.3,yMx*1.05, 0.025); ar->Draw(); ar->SetLineColor(8);
552  tx1=new TLatex(14.8,yMx*.85,"0"); tx1->Draw(); tx1->SetTextColor(8);tx1->SetTextSize(txSz);
553  tx1=new TLatex(17.1,yMx*.85,"+1.3"); tx1->Draw();tx1->SetTextColor(8);tx1->SetTextSize(txSz);
554 
555  tx1=new TLatex(20.1,yMx*.8,"Barrel"); tx1->Draw();tx1->SetTextAngle(90);tx1->SetTextSize(txSz);
556  tx1=new TLatex(19.1,yMx*.74,"Endcap"); tx1->Draw();tx1->SetTextAngle(90);tx1->SetTextSize(txSz);
557  gPad->SetGridy();
558 
559  if(jj==0 && fdWAL) {// overlay model
560  TH1F *hWALL=fdWAL->Get("modW_ALL_"+Q); assert(hWALL);
561  hWALL->Draw("same");
562  }
563 
564  }
565  //------- yield plot -------
566  cR->cd(2);
567  gPad->SetLeftMargin(0.15);
568  h=hNused[iq];
569  h->Draw(); h->SetAxisRange(0,8.4);
570  ax=h->GetXaxis();
571  float txSz=0.06;
572  ax->SetTitle("STAR #eta bins ");
573  ax->SetTitleSize(txSz); ax->SetTitleOffset(0.9);
574  ax->SetLabelSize(txSz); ax->SetLabelOffset(0.01);
575  ax=h->GetYaxis();
576  ax->SetLabelSize(txSz);
577  ax->SetTitleSize(txSz); ax->SetTitleOffset(0.85);
578  h->Draw("h text"); h->SetMarkerSize(3.); //<-- larger text font
579  h->SetFillColor(kBlue);
580  h->SetMinimum(0.8);
581  bxEG=new TBox(6.5,0.1,8.5,1e6); bxEG->SetFillColor(15); bxEG->SetFillStyle(3944);
582  bxEG->Draw();
583  gPad->SetGridy();
584  gPad->SetLogy();
585 
586  ar=new TArrow(0.6,1.3,6.3,1.3, 0.025); ar->Draw(); ar->SetLineColor(8);
587  tx1=new TLatex(0.8,2,"-0.9"); tx1->Draw(); tx1->SetTextColor(8);
588  tx1=new TLatex(5.5,2,"+1.3"); tx1->Draw();tx1->SetTextColor(8);
589  tx1=new TLatex(8.1,3,"Barrel"); tx1->Draw();tx1->SetTextAngle(90); tx1->SetTextColor(8);
590  tx1=new TLatex(7.1,3,"Endcap"); tx1->Draw();tx1->SetTextAngle(90); tx1->SetTextColor(8);
591  tx1=new TLatex(3.1,2,"STAR #eta"); tx1->Draw();tx1->SetTextColor(8); tx1->SetTextColor(8);
592 
593 }
594 
595 
596 //======================================
597 //======================================
598 //======================================
599 void printNiceTable(TString inpCore, int iq, double pol1=0.5,double pol2=0.5 ) {
600  TString Q="P"; if (iq) Q="N";
601 
602  printf("\n******* W(eta) summary for charge=%s INPUT=%s *********\n",Q.Data(), inpCore.Data());
603 
604  printf("lumi-corrections: ");
605  for(int is=0; is<mxSS; is++)printf("%.3f, ",hLumi->GetBinContent(kA[is]+1));
606 
607  printf(" applied\nstar-bin, sum , yield ++ +- -+ -- , 1/sqrt(sum), beta\n");
608  for( int starPhysEtaBin=1; starPhysEtaBin<=8; starPhysEtaBin++) {
609  double val=hNused[iq]->GetBinContent(starPhysEtaBin);
610  if(val<=0 ) val=0.1;
611  printf("%d %6.0f, ", starPhysEtaBin,val);
612  for(int is=0; is<mxSS; is++) printf("%5.0f ",SumYieldPerSpin[iq][starPhysEtaBin-1][is]);
613  printf(", %.3f, %.2f\n", 1./sqrt(val), hBeta[iq]->GetBinContent(starPhysEtaBin));
614  }
615 
616  printf("Spin results: pol1=%.2f pol2=%.2f\npolBeam-bin, events, *** AL ***,sig*sqrt(M)\n",pol1, pol2);
617  for (int k=10;k <=20; k++) { // index polBeamEtaBin
618  double eve=hNused[iq]->GetBinContent(k);
619  val=hALfin[iq]->GetBinContent(k);
620  float err=hALfin[iq]->GetBinError(k);
621  printf("%d %6.0f %.3f +/- %.3f nSig=%.1f , %.2f\n",k, eve,val,err,fabs(val)/err , err*sqrt(eve));
622  }
623 
624  printf("polBeam-bin, events,*** ALL *** ,sig*sqrt(M)\n");
625  for (int k=15;k <=20; k++) { // index polBeamEtaBin
626  double eve=hNused[iq]->GetBinContent(k);
627  val=hALfin[iq]->GetBinContent(k);
628  val=hALL[iq]->GetBinContent(k);
629  err=hALL[iq]->GetBinError(k);
630  printf("%d %6.0f %.3f +/- %.3f nSig=%.1f , %.2f\n",k, eve,val,err,fabs(val)/err , err*sqrt(eve));
631  }
632  printf("polBeam-bin, events, *** NULL *** \n");
633  for (int k=15;k <=20; k++) { // index polBeamEtaBin
634  double eve=hNused[iq]->GetBinContent(k);
635  val=hEpsNull[iq]->GetBinContent(k);
636  err=hEpsNull[iq]->GetBinError(k);
637  err=hALL[iq]->GetBinError(k);
638  printf("%d %6.0f %.3f +/- %.3f nSig=%.1f \n",k, eve,val,err,fabs(val)/err );
639  }
640 
641  printf("\******* end ************** charge=%s ********\n",Q.Data());
642 }
643 
644 //========================================
645 void saveEtaSpectra(int iq) {
646  TString Q="P"; if (iq) Q="N";
647  c3=new TCanvas("etaYiled_"+Q,"etaYield_"+Q);
648  c3->Divide(3,3); //gStyle->SetOptStat(1111);
649 
650  TGraphErrors *gr= new TGraphErrors;
651  gr->SetName( "finEta_"+Q);
652  gr->SetTitle("final W Q=P+N yield ; average STAR #eta #\pm bin rms; events");
653 
654  //..... access eta spectra ....
655  float txSz=0.06;
656  double yyMx=70.;
657  if(iq) yyMx=26;
658  for( int starPhysEtaBin=1; starPhysEtaBin<=mxEta; starPhysEtaBin++) {
659 
660  TH1F *hEta=(TH1F *)fdData->Get(Form("Eta%d/AEta%dspinLepEta_",starPhysEtaBin,starPhysEtaBin)+Q); assert(hEta);
661  TH1F *hEtaE=(TH1F *)fdData->Get(Form("Eta%d/AEta%dspinELepEta_",starPhysEtaBin,starPhysEtaBin)+Q); assert(hEtaE);
662  hEtaE->SetFillColor(24);
663 
664  hEta->Rebin(4);hEtaE->Rebin(4);
665  c3->cd(starPhysEtaBin); gPad->SetBottomMargin(0.15);
666  if( starPhysEtaBin==6 || starPhysEtaBin ==7) hEta=hEtaE;
667  hEta->Draw(); // order matters!!!
668  if( starPhysEtaBin==5) { hEta->Add(hEtaE);hEtaE->Draw("same"); }
669  float eveEta= hEta->GetEntries();
670  float avrEta= hEta->GetMean();
671  float rmsEta= hEta->GetRMS();
672  printf("etaBin=%d eve=%d Eta=%.3f +/- %.3f\n",starPhysEtaBin,eveEta, avrEta, rmsEta);
673  ttx=new TLatex(-1,yyMx*0.85,Form("%.0f eve, #eta= %.2f",eveEta,avrEta));
674  ttx->Draw();ttx->SetTextSize(1.5*txSz);ttx->SetTextColor(kBlue);
675  gr->SetPoint(starPhysEtaBin-1,avrEta,eveEta);
676  gr->SetPointError(starPhysEtaBin-1,rmsEta,sqrt(eveEta));
677  hEta->SetMaximum(yyMx);
678  hEta->SetLineWidth(2);
679  ax=hEta->GetXaxis();
680 
681  ax->SetTitleSize(txSz); ax->SetTitleOffset(0.9);
682  ax->SetLabelSize(txSz); ax->SetLabelOffset(0.01);
683  ax=hEta->GetYaxis();
684  ax->SetLabelSize(txSz);
685  ax->SetTitleSize(txSz); ax->SetTitleOffset(0.85);
686 
687  gPad->SetGridy();
688  hEta->Write() ; hEta->SetAxisRange(-1.1,1.6);
689  }
690 
691  // show yield vs. float-eta
692  c3->cd(9); gr->Draw("A P");
693  gPad->SetLogy();
694  //gr->Print();
695  gr->Write();
696 
697 }
TDataSet::Clone
virtual TObject * Clone(const char *newname="") const
the custom implementation fo the TObject::Clone
Definition: TDataSet.cxx:308
Generated by   doxygen 1.8.5