StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
addRuns.C
1 const int mxBx=120, mxPol=4;
2 
3 char *cpolBY[mxPol]={"B+Y+","B+Y-","B-Y+","B-Y-"};
4 char bXselection[mxBx];
5 int polPattBY[mxBx];
6 double Lum[mxPol],Lum2[mxPol];
7 TFile *outH=0;
8 TString cutL=" ";
9 
10 
11 //-------------------------------------
12 //-------------------------------------
13 void updatePattern(char *run) {
14  int irun=atoi(run+1);
15  assert(irun>0);
16  printf("Update pattern for run=%d\n",irun);
17 
18  char * polPattY,* polPattB,*dum;
19 
20  char * //yellow
21  selPatt="a.........k....BBBBBa.........k....pppppa.........k....YYYYY";
22  dum ="| | | | | | |";
23  dum ="1 11 21 31 41 51 60";
24 
25  if(irun<3017036) { // (STAR) Pattern1
26  polPattY= "--++--++--++--++--++--++--++--++--++--++--++--++--++--+aaaaa";
27  polPattB= "-+-+-+-+-+-+-+-aaaaa-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+";
28  dum = "| | | | | | |";
29  dum = "1 11 21 31 41 51 60";
30  } else { // (STAR) Pattern2 with 3 unpolarized bXings
31  polPattY= "0-+-+-+-+-+-+-+-+-+-0+-+-+-+-+-+-+-+-+-+0-+-+-+-+-+-+-+aaaaa";
32  polPattB= "0++--++--++--++aaaaa0--++--++--++--++--+0+--++--++--++--++--";
33  dum = "| | | | | | |";
34  dum = "1 11 21 31 41 51 60";
35  }
36 
37  int i;
38  for(i=0;i<60;i++) {
39  bXselection[2*i]=selPatt[i];
40  bXselection[2*i+1]='i'; // ignore
41  }
42 
43  for(i=0;i<60;i++) {
44  int k=2*i;
45  polPattBY[k]=-2; //default
46  if (polPattB[i]=='+' && polPattY[i]=='+' ) polPattBY[k]=0;
47  else if(polPattB[i]=='+' && polPattY[i]=='-' ) polPattBY[k]=1;
48  else if(polPattB[i]=='-' && polPattY[i]=='+' ) polPattBY[k]=2;
49  else if(polPattB[i]=='-' && polPattY[i]=='-' ) polPattBY[k]=3;
50  }
51 
52 }
53 
54 
55 
56 //-------------------------------------
57 // M A I N
58 //-------------------------------------
59 addRuns(TString fill=0, TString runL=0, TString wrkDir=0) {
60  gStyle->SetPalette(1,0);
61  gStyle->SetOptStat(1111111);
62 
63  //fill="F2201x"; runL="R2362128 R777";
64  //wrkDir="/star/data04/sim/balewski/LcpRun2/setEta10/";
65 
66  // create output histo
67  TString outF=wrkDir+fill+".hist.root";
68  //TString outF=fill+".hist.root";
69  outH=new TFile(outF,"RECREATE");
70  assert(outH->IsOpen());
71  printf("save outH -->%s\n", outF.Data());
72 
73 
74  // ----- access scaler data
75  TString inpScal="JoannaScal/scalVsRun.hist.root";
76  scal=new TFile(inpScal);
77  assert(scal->IsOpen());
78 
79  char *run=strtok(runL.Data()," "); // init 'strtok'
80 
81  float sumIn=0;
82  int nRun=0, nLum=0; // counts runs & scaler data for them
83  do { // loop over runs
84  printf("\n\n process run %02d '%s' \n",nRun,run);
85  updatePattern(run);
86 
87  //.................. access LCP histos
88  TString lcpHist=wrkDir+run+".hist.root";
89  TFile * lcp=new TFile(lcpHist);
90  assert(lcp->IsOpen());
91  if(!lcp->IsOpen()) {
92  printf("#fail %s does not open, skip\n",lcpHist.Data());
93  continue;
94  }
95  nRun++;
96 
97 
98  // ................ access scale data (if exist)
99  TH1F *hSc=(TH1F*)scal->Get(run);
100  if(hSc) {
101  printf(" use absolute normalization\n");
102  sumScal(hSc);
103  calcNorm(run);
104  nLum++;
105  } else { // do not scal by lumin
106  printf(" NO Scaler data, --> no renormalization\n");
107  for(int i=0;i<mxPol;i++) {Lum[i]=Lum2[i]=1.;}
108  }
109 
110  //lcp->ls();
111 
112  TList *top = lcp->GetListOfKeys();
113  TIter iter(top);
114  TObject*ob=0;
115  // ........... loop over histo with LCP
116  while( ob=iter()) {
117  char *name=ob->GetName();
118  if(strstr(name,"PhiBx")==0) continue;
119  // if(strstr(name,"PhiBxPt1")==0 && strstr(name,"PhiBxAll")==0) continue;
120 
121  TH2F *hIn=(TH2F*)lcp->Get(name);
122  printf(" %10s integral=%.1f\n",name, hIn->Integral());
123  if(strstr(name,"PhiBxAll"))sumIn+=hIn->Integral();
124 
125  if(nRun==1) createHisto1D(hIn);
126  TH1F * h[mxPol];
127  TString cut=name+5;
128  fetch1D(cut,h);
129  accumulateLcp(hIn,h);
130 
131  } // end of loop over histos in one run
132  lcp->Close();
133  delete lcp;
134 
135 
136  } while(run=strtok(0," ")); // advance by one nam
137 
138  printf("acumulation done nRun=%d nLum=%d\n", nRun, nLum);
139  assert(nLum==0 || nRun==nLum);
140 
141  if(nLum==0) autoLum("All");
142  printf("left2=%s=\n",cutL.Data());
143  printf("#fill %s , inp= %d , ",fill.Data(),sumIn);
144  float sumAll=summaryLcp();
145  printf(" AllEff= %.2f\n",sumAll/sumIn );
146  outH->Write();
147 
148  //printLcp("PhiBxAll");
149 
150 
151  // printLcp("PhiBxPtM");
152 
153  // outH->ls();
154  (outH->Get("PhiB+Y-All"))->Draw();
155  TH1F*h1=(TH1F*) outH->Get("PhiB+Y+All");
156  h1->Draw("same"); h1->SetLineColor(kRed);
157 
158 
159  // hsc->Draw();
160 
161 }
162 
163 
164 
165 //-------------------------------------
166 //-------------------------------------
167 int bx2pol(int bx1) {
168  if(bXselection[bx1-1]!='.') return -3;
169  int pol= polPattBY[bx1-1];
170  assert(pol<mxPol);
171  return pol;
172 }
173 
174 //-------------------------------------
175 //-------------------------------------
176 void createHisto1D(TH2F *hIn){
177  TString name=hIn->GetName();
178  TString cut=name.Data()+5;
179  cutL+=" "+cut;
180  printf("create outHist for %s\n",name.Data());
181  TAxis *phiAx=hIn->GetYaxis();
182 
183  outH->cd();
184  int k;
185  for(k=0;k<mxPol;k++) {
186  TString name1="Phi";
187  name1+=cpolBY[k];
188  name1+=cut;
189  char tit2[100]={"aaa bbb"};
190  sprintf(tit2,"LCP vs. Phi/rad, pol=%s, cut=%s",cpolBY[k],cut.Data());
191  printf("k=%d %s '%s'\n",k,name1.Data(), tit2);
192  TH1F *h=new TH1F(name1,tit2,phiAx->GetNbins(),phiAx->GetXmin(), phiAx->GetXmax());
193  h->Sumw2(); // be carefull, not use Fill(x,w) !
194  }
195 }
196 
197 
198 //-------------------------------------
199 //-------------------------------------
200 void accumulateLcp(TH2F *hIn, TH1F**h){
201  TString name=hIn->GetName();
202  // printf("accu histo %s\n",name.Data());
203 
204  // fetch proper 1D histo
205  int k;
206  for(k=0;k<mxPol;k++) {
207  TString name1=hIn->GetName();
208  name1.ReplaceAll("Bx",cpolBY[k]);
209  h[k]=(TH1F *)outH->Get(name1);
210  assert(h[k]);
211  }
212 
213  assert(hIn);
214  TAxis *X=hIn->GetXaxis();
215  int nx=X->GetNbins();
216  TAxis *Yax=hIn->GetYaxis();
217  assert(Yax);
218  int ny=Yax->GetNbins();
219  int ix,iy;
220 
221  double nSel=0;
222  for(iy=1; iy<=ny; iy++) {
223  float phi=Yax->GetBinCenter(iy);
224  //intf("iy=%d phi/deg=%f\n", iy,phi/3.1416*180);
225  for(ix=1;ix<=nx;ix++) {
226  int pol=bx2pol(ix);
227  if(pol<0) continue;
228  double val=hIn->GetBinContent(ix,iy);
229  if (val==0) continue;
230  nSel+=val;
231  TH1F * h1=h[pol];
232  h1->AddBinContent(iy,val/Lum[pol]);
233  // do error propagation by hand
234  double err=h1->GetBinError(iy);
235  h1->SetBinError(iy,sqrt(err*err+val/Lum2[pol]));
236  // printf("add %d %f %d\n",pol,val, h[pol]->Integral());
237  }
238  }
239 
240  // cross checks
241 
242  return;
243  printf("%s --> %.1f nSel=%.1f\n",hIn->GetName(), hIn->Integral(),nSel);
244  for(k=0;k<mxPol;k++) {
245  TH1F *h1=h[k];
246  printf("%s --> %.0f\n",h1->GetName(), h1->Integral());
247  }
248 
249 }
250 
251 
252 //-------------------------------------
253 //-------------------------------------
254 void printLcp(char *cut){
255  printf("print histo %s, format: kPhi polBY ++ +- -+ --\n#tot: ",cut);
256 
257  TH1F * h[mxPol];
258 
259  // fetch proper 1D histo
260  int k;
261  double sum=0;
262  for(k=0;k<mxPol;k++) {
263  TString name1=cut;
264  name1.ReplaceAll("Bx",cpolBY[k]);
265  h[k]=(TH1F *)outH->Get(name1);
266  assert(h[k]);
267  sum+=h[k]->Integral();
268  printf(" %.1f ", h[k]->Integral());
269  }
270 
271  printf("\n#rel: ",cut);
272  for(k=0;k<mxPol;k++) {
273  printf(" %.3f ", h[k]->Integral()*4/sum);
274  }
275  printf("\n");
276  TAxis *X=h[0]->GetXaxis();
277  int nx=X->GetNbins();
278  int ix;
279 
280  for(ix=1; ix<=nx; ix++) {
281  printf("%2d ",ix);
282  for(k=0;k<mxPol;k++) {
283  float val=h[k]->GetBinContent(ix);
284  float err=h[k]->GetBinError(ix);
285  printf("%8.1f +/- %5.1f, ",val,err);
286  }
287  printf("\n");
288  }
289 }
290 
291 //-------------------------------------
292 //-------------------------------------
293 void fetch1D(TString cut,TH1F **h) {// fetch proper 1D histo
294  int k;
295  for(k=0;k<mxPol;k++) {
296  TString name1="Phi";
297  name1+=cpolBY[k];
298  name1+=cut;
299  h[k]=(TH1F *)outH->Get(name1);
300  assert(h[k]);
301  }
302 }
303 
304 //-------------------------------------
305 //-------------------------------------
306 void sumScal(TH1F *h) {
307  memset(Lum,0,sizeof(Lum)); // clear normalization
308  memset(Lum2,0,sizeof(Lum2));
309 
310  assert(h);
311  TAxis *X=h->GetXaxis();
312  int nx=X->GetNbins();
313  int ix;
314  int nBx=0;
315  for(ix=1;ix<=nx;ix++) {
316  int pol=bx2pol(ix);
317  if(pol<0) continue;
318  assert(pol<mxPol);
319  Lum[pol]+=h->GetBinContent(ix);
320  nBx++;
321  // printf("%3d %c %2d %8.1f %d \n",ix,bXselection[ix-1],polPattBY[ix-1], h->GetBinContent(ix),nBx);
322  }
323 
324  int i;
325  printf("Lum: ");
326  for(i=0;i<mxPol;i++)
327  printf("%d ",Lum[i]);
328 
329  printf(" usedBx=%d\n",nBx);
330 
331 }
332 
333 //-------------------------------------
334 //-------------------------------------
335 void calcNorm(char *name){
336  int i;
337  float sum=0;
338  printf("\nyield ");
339  for(i=0;i<mxPol;i++){
340  sum+=Lum[i];
341  printf("%d ",Lum[i]);
342  }
343 
344  int id=atoi(name+1);
345  //float day=(id-3000000.)/1000.;
346 
347  printf("\n#norm %s %d ",name , id);
348  sum/=4.;
349  // renormalize
350  for(i=0;i<mxPol;i++){
351  Lum[i]/=sum;
352  Lum2[i]=Lum[i]*Lum[i];
353  printf("%.4f ",Lum[i]);
354  }
355  printf("\n");
356 }
357 
358 
359 
360 //-------------------------------------
361 //-------------------------------------
362 void autoLum( char *cut0) {
363  printf("autoLuminosiy using cut=%s\n",cut0);
364 
365  TH1F * h[mxPol];
366  fetch1D(cut0,h);
367  int k;
368  for(k=0;k<mxPol;k++) {
369  Lum[k]=h[k]->Integral();
370  }// end of loop over pol
371 
372  calcNorm("AutoLum7777");
373 
374  // scale all histos for all cuts
375  TString cut1=" "+cutL;
376  char *cut=strtok(cut1.Data()," "); // init 'strtok'
377  do {
378  TH1F * h[mxPol];
379  fetch1D(cut,h);
380  int k;
381  for(k=0;k<mxPol;k++) {
382  h[k]->Scale(1./Lum[k]);
383  }// end of loop over pol
384  } while(cut=strtok(0," ")); // advance by one nam
385 }
386 
387 
388 //-------------------------------------
389 //-------------------------------------
390 float summaryLcp(){
391  TString cut1=cutL;
392  //printf("summary for cutL=%s\n",cut1.Data());
393  char *cut=strtok(cut1.Data()," "); // init 'strtok'
394  int sumAll=-1;
395  do {
396  // printf("cut: %s\n",cut);
397  float sum[mxPol];
398  float sum0=0;
399  TH1F * h[mxPol];
400  fetch1D(cut,h);
401  int k;
402  for(k=0;k<mxPol;k++) {
403  // printf(" %s %.1f -->",name1.Data(),h->Integral());
404  sum[k]=h[k]->Integral();
405  sum0+=sum[k];
406  }// end of loop over pol
407 
408  printf("%s= %.1f , ",cut,sum0);
409  if(strstr(cut,"All"))sumAll=sum0;
410 
411  } while(cut=strtok(0," ")); // advance by one nam
412 
413  return sumAll;
414 }
415 
416 
417 
418 
419 
420 
421 
422 
Generated by   doxygen 1.8.5