StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
rdN2AL.C
1 // converts raw inputs yields in to physics asymetries
2 // version 1.3
3 
4 /* mapping of spin4-index to helicity at STAR */
5 enum { ka=10, /* STAR pol B+ Y + */
6  kb=9, /* STAR pol B+ Y - */
7  kc=6, /* STAR pol B- Y + */
8  kd=5, /* STAR pol B- Y - */
9  mxSS=4,
10  mxQ=2 /*0=W+, 1=W- */
11 };
12 double pol1=0.376; // blue
13 double pol2=0.395; // yellow
14 // unpol background:
15 double betaVal[mxQ]={0.936,0.837}; //beta=S/(S+B) for W+ & W-, tau not included
16 double betaErr[mxQ]={0.017, 0.032};
17 // pol background:
18 double alphaVal[mxQ]={-0.0016, -0.0048};
19 double alphaErr[mxQ]={ 0.0008, 0.0024};
20 
21 
22 int kA[mxSS]={ka,kb,kc,kd};
23 TH1F *hData[mxQ];
24 TH1F *hAsy[mxQ];
25 TH1F *hLum=new TH1F("hLum","hLum", 50,0.5,50.5); // for lumi monitor and
26 
27 void rdN2AL(TString inpCore="run9setP1234") {
28  TString iPath="/star/data05/scratch/balewski/2009-Wana-SL09g-p/data/";
29  //iPath="sortMarch22/";
30  iPath="sortMay3/";
31 
32  TString oPath="out/";
33 
34  for(int iq=0;iq<mxQ;iq++) {
35  char *tt1="hDataP";
36  if(iq==1) tt1="hDataN";
37  hData[iq]=new TH1F(tt1,tt1,30,0.5,30.5);
38  tt1="hAsyP"; char *tt2="Positive charge";
39  if(iq==1){ tt1="hAsyN"; tt2="Negative charge";}
40  hAsy[iq]=new TH1F(tt1,tt2,7,0,7);
41  hAsy[iq]->GetXaxis()->SetLabelSize(0.065);
42  hAsy[iq]->SetLineWidth(2);
43  char *key[]={"AL blue","AL yell","AL avr","ALL","null test","AL* "};
44  for(int i=0;i<6;i++) hAsy[iq]->Fill(key[i],0.); // preset the order of keys
45  hAsy[iq]->SetMinimum(-0.95); hAsy[iq]->SetMaximum(0.95);
46  }
47 
48  TString fullOutName=oPath+inpCore+".wasy.hist.root";
49 
50  TString fullInpName=iPath; fullInpName+=inpCore;
51  fullInpName+=".wana.hist.root";
52  fd=new TFile(fullInpName);
53  if(! fd->IsOpen()) {
54  printf("EROR: input histo file not found, quit\n",fullInpName.Data());
55  return;
56  } else {
57  printf("Opened: %s\n",fullInpName.Data());
58  }
59 
60  getLumi( fd);
61  // hLum->Draw(); gPad->SetLogy();
62 
63  for(int iq=0;iq<mxQ;iq++) {// Q=+, -
64  printf("........processing iQ=%d .....\n",iq);
65  getSignal(fd,iq,"AspinY2"); //Aspin2=default, select here B,C for East/West
66  computeAsy(iq);
67  }
68  fd->Close();
69 
70  fdo=new TFile(fullOutName,"recreate"); assert(fdo->IsOpen());
71  hLum->Write();
72  hData[0]->Write();
73  hData[1]->Write();
74  hAsy[0]->Write();
75  hAsy[1]->Write();
76 
77  //........ nice plot with summary
78  can=new TCanvas(inpCore, inpCore,500,400);
79  gStyle->SetOptStat(0);
80  TPad *c=makeTitle(can,(iPath+inpCore).Data());
81  c->SetFillColor(kWhite);
82  c->Divide(2,1);
83  ln=new TLine(0,0,7,0); ln->SetLineStyle(2);
84 
85  for(int iq=0;iq<mxQ;iq++) {// Q=+, -
86  c->cd(1+iq); hAsy[iq]->Draw(); ln->Draw();
87  gPad->SetFillStyle(0);
88  TString tt=hAsy[iq]->GetTitle(); tt+=", unpol yield=";
89  float sum=hData[iq]->GetBinContent(5); tt+=(int)sum;
90  hAsy[iq]->SetTitle(tt);
91 
92  // add special markers for 3 physcis observables
93  int myMark[mxQ][5]={{21,28,25},{20,28,24}};
94  for(int kk=2; kk<=4;kk++) {
95  hx=(TH1*) hAsy[iq]->Clone();
96  hx->Draw("same"); hx->SetMarkerSize(2);
97  hx->SetAxisRange(kk,kk);
98 
99  hx->SetMarkerStyle(myMark[iq][kk-2]);
100  }
101 
102  }
103 
104 }
105 
106 //**********************************************************
107 //.................
108 void getLumi( TFile *fd) {
109  TH1F *h1=fd->Get("AspinY1"); assert(h1);// h1->Draw();
110 
111  double sum=0;
112  for( int k=1; k<=4; k++ ){
113  double y= h1->GetBinContent(kA[k-1]+1);
114  sum+=y;
115  hLum->SetBinContent(k,y);
116  }
117  hLum->SetBinContent(5,sum);
118 
119  printf(" lum total=%.0f\n", sum);
120  sum/=4.;
121  for( int k=1; k<=4; k++ ){
122  double y= hLum->GetBinContent(k);
123  printf("k=%d %.0f %.4f (spin4=%d)\n",k,y, y/sum,kA[k-1]);
124  hLum->SetBinContent(k+10,y/sum);
125  }
126  hLum->SetBinContent(6,pol1);
127  hLum->SetBinContent(7,pol2);
128  printf("Polarization beam_1=%.3f beam_2=%.3f\n", pol1, pol2);
129 
130  // fd->ls();
131 }
132 
133 //**********************************************************
134 
135 //.................
136 void getSignal( TFile *fd, int iq, char *hcore) {
137  TString hname=hcore;
138  if(iq==0) hname+="_P";
139  if(iq==1) hname+="_N";
140  TH1F *h1=fd->Get(hname); assert(h1); //h1->Draw();
141 
142  double sum=0, sumL=0;
143  for( int k=1; k<=4; k++ ){
144  double y= h1->GetBinContent(kA[k-1]+1);
145 #if 0 // add ET>50
146  if(k==1 && iq==0) y+=4;
147  if(k==2 && iq==0) y+=3;
148  if(k==3 && iq==0) y+=1;
149  if(k==4 && iq==0) y+=2;
150  if(k==2 && iq==1) y+=2;
151  if(k==3 && iq==1) y+=1;
152 #endif
153 
154  printf("k=%d yield=%.0f\n",k,y);
155  sum+=y;
156 
157  hData[iq]->SetBinContent(k,y);
158  double L= hLum->GetBinContent(k+10);
159  double yL=y/L;
160  double syL=sqrt(y)/L;
161  sumL+=yL;
162  hData[iq]->SetBinContent(k+10,yL);
163  hData[iq]->SetBinError(k+10,syL);
164  }
165  hData[iq]->SetBinContent(5,sum);
166  hData[iq]->SetBinContent(15,sumL);
167 
168  // hData[iq]->Draw("e text");
169 }
170 
171 //**********************************************************
172 //---------------------
173 void computeAsy(int iq) {
174  TH1F *hD= hData[iq];
175  TH1F *hA= hAsy[iq];
176 
177  double P1=hLum->GetBinContent(6);
178  double P2=hLum->GetBinContent(7);
179 
180  double beta=betaVal[iq], sBeta=betaErr[iq];
181  double alpha=alphaVal[iq], sAlpha=alphaErr[iq];
182 
183 
184  double Ma=hD->GetBinContent(11), sMa=hD->GetBinError(11), VMa=sMa*sMa;
185  double Mb=hD->GetBinContent(12), sMb=hD->GetBinError(12), VMb=sMb*sMb;
186  double Mc=hD->GetBinContent(13), sMc=hD->GetBinError(13), VMc=sMc*sMc;
187  double Md=hD->GetBinContent(14), sMd=hD->GetBinError(14), VMd=sMd*sMd;
188 
189  printf(" yield total=%.0f S/(S+B)=%.2f +/-%f alpha=%f +/-%f\n", Ma+Mb+Mc+Md,beta,sBeta,alpha,sAlpha);
190 
191  //..................... AL for beam 1
192  double eps, sEps, AL,sAL;
193 
194  doEps_I( Ma+ Mb, Mc+ Md,
195  VMa+VMb, VMc+VMd, eps, sEps );
196 
197  doPolBckgCorr(eps, sEps, P1, alpha, sAlpha, beta, sBeta,
198  AL, sAL);
199  hD->SetBinContent(21,eps); hA->SetBinContent(1,AL);
200  hD->SetBinError (21,sEps); hA->SetBinError (1,sAL);
201  // printf("eps1 %f %f %f\n", eps, sEps, 1./sqrt(Ma+Mb+Mc+Md));
202  printf("AL(1) %.3f %.3f nSig=%.2f\n", AL, sAL, AL/sAL);
203 
204  //...................... AL for beam 2
205  doEps_I( Ma+ Mc, Mb +Md,
206  VMa+VMc, VMb+VMd, eps, sEps );
207 
208  doPolBckgCorr(eps, sEps, P2, alpha, sAlpha, beta, sBeta,
209  AL, sAL);
210 
211  hD->SetBinContent(22,eps); hA->SetBinContent(2,AL);
212  hD->SetBinError (22,sEps); hA->SetBinError (2,sAL);
213  // printf("eps2 %f %f %f\n", eps, sEps, 1./sqrt(Ma+Mb+Mc+Md));
214  printf("AL(2) %.3f %.3f nSig=%.2f\n", AL, sAL, AL/sAL);
215 
216  //...................... average AL for beam 1+2
217  doEps_II( Ma, Md, Mb +Mc,
218  VMa, VMd, VMb+VMc, eps, sEps );
219  doPolBckgCorr(eps, sEps, (P1+P2)/2., alpha, sAlpha, beta, sBeta,
220  AL, sAL);
221 
222  hD->SetBinContent(23,eps); hA->SetBinContent(3,AL);
223  hD->SetBinError (23,sEps); hA->SetBinError (3,sAL);
224  // printf("eps3 %f %f %f\n", eps, sEps, 1./sqrt(Ma+Mb+Mc+Md)/sqrt(2.));
225  printf("AL(1+2) %.3f $\\pm$ %.4f nSig=%.2f\n", AL, sAL, AL/sAL);
226 
227  //...................... DSS ALL -regular
228  doEps_I( Ma+ Md, Mb +Mc,
229  VMa+VMd, VMb+VMc, eps, sEps );
230  double ALL, sALL;
231  doPolBckgCorr(eps, sEps, P2*P2, alpha, sAlpha, beta, sBeta,
232  ALL, sALL);
233 
234 
235  hD->SetBinContent(24,eps); hA->SetBinContent(4,ALL);
236  hD->SetBinError (24,sEps); hA->SetBinError (4,sALL);
237  // printf("eps4 %f %f %f\n", eps, sEps, 1./sqrt(Ma+Mb+Mc+Md));
238  printf("ALL(reg) %.2f %.2f nSig=%.1f\n", ALL, sALL, ALL/sALL);
239 
240  //...................... ALL : PV barr = null test
241  doEps_I( Mb, Mc,
242  VMb, VMc, eps, sEps );
243 
244  hD->SetBinContent(25,eps); hA->SetBinContent(5,eps);
245  hD->SetBinError (25,sEps); hA->SetBinError (5,sEps);
246  printf("eps5= Null %.3f $\\pm$ %.3f nSig=%.2f\n", eps, sEps, eps/ sEps );
247 
248  //...................... PV ALL, both beams
249  doEps_I( Ma , Md,
250  VMa, VMd, eps, sEps );
251  doPolBckgCorr(eps, sEps, (P1+P2), alpha, sAlpha, beta, sBeta,
252  AL, sAL);
253 
254  hD->SetBinContent(26,eps); hA->SetBinContent(6,AL);
255  hD->SetBinError (26,sEps); hA->SetBinError (6,sAL);
256  // printf("eps6 %f %f %f\n", eps, sEps, 1./sqrt(Ma+Mb+Mc+Md));
257  printf("AL*(eps6) %.3f %.3f nSig=%.1f\n", AL, sAL, AL/sAL);
258 
259 }
260 
261 //------------
262 void doEps_I( double a, double b, double va, double vb,
263  double &eps, double &sEps) {
264  double sum=a+b;
265  eps= (a-b)/sum;
266  double xx=b*b*va+ a*a*vb;
267  sEps= sqrt( 4 * xx/sum/sum/sum/sum);
268 }
269 
270 //------------
271 void doEps_II( double a, double b,double c, double Va, double Vb, double Vc,
272  double &eps, double &sEps) {
273  double sum=a+b+c;
274  eps= (a-b)/sum;
275  double vEps=(4*a*c*Vb + c*c*(Va + Vb) + b*b*(4*Va + Vc) +
276  a*a*(4*Vb + Vc) + b*(4*c*Va - 2*a*Vc))/sum/sum/sum/sum;
277  sEps= sqrt(vEps);
278 }
279 
280 //------------
281 void doPolBckgCorr( double eps, double sEps,double P, double a, double sa, double b, double sb,
282  double &AL, double &sAL) {
283  AL=(eps/P-a)/b;
284  double v1=sEps*sEps/P/P;
285  double v2=sa*sa;
286  double v3=AL*AL*sb*sb;
287  sAL=sqrt(v1+v2+v3)/b;
288 }
289 
290 
291 //------------------------
292 Void splitPadX(float x, TPad **cL, TPad **cR) {
293  (*cL) = new TPad("padL", "apdL",0.0,0.,x,0.95);
294  (*cL)->Draw();
295  (*cR) = new TPad("padL", "apdL",x+0.005,0.,1.0,0.95);
296  (*cR)->Draw();
297 }
298 
299 //------------------------
300 TPad *makeTitle(TCanvas *c,char *core) {
301  c->Range(0,0,1,1);
302  TPad *pad0 = new TPad("pad0", "apd0",0.0,0.95,1.,1.);
303  pad0->Draw();
304  pad0->cd();
305 
306  TPaveText *pt = new TPaveText(0,0.,1,1,"br");
307  pt->Draw();
308  TDatime dt;
309  TString txt2=core;
310 
311  txt2+=", ";
312  txt2+=dt.AsString();
313  pt->AddText(txt2);
314  txt2="--";
315  pt->AddText(txt2);
316 
317  c->cd();
318  pad = new TPad("pad1", "apd1",0.0,0.0,1,.95);
319  pad->Draw();
320  return pad;
321 }
Generated by   doxygen 1.8.5