StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
doAvrPed.C
1 /* notes, R9067013 OLD
2  BPRS:
3  kBad=0x1: yield in ADC range [100,300] 100 counts, input was ~220 events
4  7 tiles: softID=3301-4,3322-4
5  kBad=0x2: yiled @ [mean +/- 9 ADC] below 70%
6  0 tiles
7  kBad=0x4: bad chi2/dof >10
8  42-7=35 tiles, softID=3863 ...3940, about 50%
9  kBad=x8 = ped low & high, set wide [105,295]
10  0 tiles
11  kBad=0x10 sticky low bit @ pedestal area
12  85-42=43 , tiles softID=3862 ...3940, about 50%
13  kBad=0x20 sig ped >2.5 ADC
14  132-85=47 tiles, softID=3021-3032, 3041-3048, 3052,3103,3621-3652
15 
16 
17 */
18 
19 const int mxTP=2, mxBTiles=4800;
20 char cTile[mxTP]={'T','P'};
21 char *cTile4[mxTP]={"BTOW","BPRS"};
22 TH1F * hPed, *sigP, *hChi,*hSig;
23 TH1I * hStat;
24 TGraphErrors *grP;
25 TFile *fd2=0, *fd1=0;
26 TH1F *hPeak; // yield of pedestal
27 TCanvas *c=0;
28 int isRawAdc=0; // 1=raw, 0=pedSub
29 
30 doAvrPed( int k=10,int run=9067013) {
31  gStyle->SetPalette(1,0);
32  gStyle->SetOptFit(1);
33  gStyle->SetOptStat(1001110);
34 
35  char *pathIn="outA3/";
36  char *pathOut="iter4/";
37  char txt[1000], txt2[1000];
38 
39  sprintf(txt,"%s/R%dp.barCal.-1.hist.root",pathIn,run);
40 
41  fd1=new TFile(txt); assert(fd1->IsOpen());
42  // fd1.ls(); return;
43  sprintf(txt,"%s/pedBprsR%davr.hist.root",pathOut,run);
44  fd2=new TFile(txt,"recreate");
45  assert(fd2->IsOpen());
46  int i;
47  grP=new TGraphErrors;
48  // gr->SetMarkerStyle(20+m);
49  grP->SetMarkerSize(0.5);
50 
51  TH2F * h2T= (TH2F *)fd1->Get("BTOW_c0"); assert(h2T);
52  TH2F * h2P= (TH2F *)fd1->Get("BPRS_c0"); assert(h2P);
53  h2P->GetXaxis()->SetLabelSize(0.08);
54  h2P->GetYaxis()->SetLabelSize(0.08);
55  h2P->Draw("colz");
56  // Printf("XXX %f\n", h2P->GetYaxis()->SetTitleSize());
57 
58  // return;
59  int id1=1, id2=id1+k-1;
60  if(k<0) { id1=1; id2=mxBTiles; }
61  c=new TCanvas("aa","aa",400,300);
62  // id1=id2=k;//tmp
63  fitPedBT(1,run,h2P,id1,id2 );
64  gPad->SetLogy();
65 
66  // save output histo
67  fd2->cd(); hPed->Write();
68  hStat->Write();
69  hChi->Write();
70  hSig->Write();
71  hPeak->Write();
72 
73  // return;
74  c=new TCanvas(fd1->GetName(),fd1->GetName(),1000,800);
75  c->Divide(1,5);
76  gStyle->SetOptStat(10);
77  c->cd(1); h2P->Draw("colz"); h2P->SetAxisRange(id1,id2);
78  hPed->Draw("same");
79 
80  c->cd(2); hStat->Draw(); hStat->SetAxisRange(id1,id2);
81  hStat->SetMinimum(.5);
82  if(hStat->GetMaximum()>10) gPad->SetLogy();
83  gPad->SetGrid();
84 
85  c->cd(3);
86  hChi->Draw(); hChi->SetAxisRange(id1,id2);// hChi->SetMaximum(20);
87  gPad->SetGrid();
88 
89  c->cd(4);
90  hSig->Draw(); hSig->SetAxisRange(id1,id2); hSig->SetMaximum(5);
91  gPad->SetGrid();
92 
93  c->cd(5);
94  hPeak->Draw(); hPeak->SetAxisRange(id1,id2); hPeak->SetMinimum(0.7);
95  gPad->SetGrid();
96 
97  // gPad->SetLogy();
98  // c->cd(1); hPed->Draw(); hPed->SetAxisRange(id1,id2);
99 
100 
101  return;
102  fd2->cd();
103  hres->Write();
104  grPedR[0]->Write();
105 
106 }
107 
108 
109 
110 //==============================================
111 //==============================================
112 //==============================================
113 void fitPedBT( int itp,int run,TH2F *h2, int id1, int id2) {
114  if(id2>18000) id2=18000;
115  char txt[1000], txt2[1000];
116 
117  sprintf(txt,"ped%s",cTile4[itp]);
118  sprintf(txt2,"ped %s R%d; %s soft ID; pedestal +/- fit err (ADC)",cTile4[itp],run,cTile4[itp]);
119  hPed=new TH1F(txt,txt2,mxBTiles,0.5,mxBTiles+0.5);
120  hPed->Sumw2();
121  hPed->GetXaxis()->SetLabelSize(0.08);
122  hPed->GetYaxis()->SetLabelSize(0.08);
123 
124  sprintf(txt,"stat%s",cTile4[itp]);
125  sprintf(txt2,"status %s R%d; %s soft ID; jan status",cTile4[itp],run,cTile4[itp]);
126  hStat=new TH1I(txt,txt2,mxBTiles,0.5,mxBTiles+0.5); hStat->Reset(-1); // default is bad
127  hStat->GetXaxis()->SetLabelSize(0.08);
128  hStat->GetYaxis()->SetLabelSize(0.08);
129 
130  sprintf(txt,"chi%s",cTile4[itp]);
131  sprintf(txt2,"chi2/DOF %s R%d; %s soft ID; ped Chi2/DOF",cTile4[itp],run,cTile4[itp]);
132  hChi=new TH1F(txt,txt2,mxBTiles,0.5,mxBTiles+0.5); hStat->Reset(-1); // default is bad
133  hChi->GetXaxis()->SetLabelSize(0.08);
134  hChi->GetYaxis()->SetLabelSize(0.08);
135  hChi->GetYaxis()->SetTitleSize(0.2);
136 
137  sprintf(txt,"sigPed%s",cTile4[itp]);
138  sprintf(txt2,"sigma(ped) %s R%d; %s soft ID; sig(ped) ADC",cTile4[itp],run,cTile4[itp]);
139  hSig=new TH1F(txt,txt2,mxBTiles,0.5,mxBTiles+0.5); hStat->Reset(-1); // default is bad
140  hSig->GetXaxis()->SetLabelSize(0.08);
141  hSig->GetYaxis()->SetLabelSize(0.08);
142 
143  hPeak=new TH1F("hpedPeak", "integral of pedestal peak;soft ID", mxBTiles,0.5,mxBTiles+0.5);
144  hPeak->GetXaxis()->SetLabelSize(0.08);
145  hPeak->GetYaxis()->SetLabelSize(0.08);
146 
147 
148  float par_nsig=3; // integration range for QA
149  float par_rms=2; // to approximate fit range
150  float cut_yield0=100;
151  float cut_yieldR1=0.7;
152  float cut_ch2min=0.1;
153  float cut_ch2ndf=200.;
154  float cut_pedL=105;
155  float cut_pedH=295;
156  float cut_pedH=295;
157  float cut_sigPed=2.7;
158 
159  // float cut_minYiled=50;
160 
161  axX=h2->GetXaxis();
162  float x1=axX->GetXmin();
163  float x2=axX->GetXmax();
164  int nbX=axX->GetNbins();
165  printf("X-axis range --> [%.1f, %.1f], nb=%d %s\n",x1,x2,nbX,axX->GetTitle());
166 
167  axY=h2->GetYaxis();
168  float y1=axY->GetXmin();
169  float y2=axY->GetXmax();
170  int nbY=axY->GetNbins();
171  printf("Y-axis range --> [%.1f, %.1f], nb=%d\n",y1,y2,nbY);
172 
173  TH1F*h=new TH1F("h1","h1",nbY,y1,y2); // working histo for 1-D spectrum
174  // do projections
175  int ih;
176  for(ih=id1; ih<=id2; ih++) {
177  char txt1[100], txt2[1000];
178  sprintf(txt1,"id%d",ih);
179  sprintf(txt2,"%s soft id=%d;%s ",cTile4[itp],ih,axY->GetTitle());
180  h->SetTitle(txt2);
181  h->Reset(); h->SetAxisRange(y1,y2);
182  int i;
183  int kBad=1;
184  for(i=1;i<=nbY;i++) h->SetBinContent(i,h2->GetBinContent(ih,i));
185  //for(i=1;i<=nbY;i++) printf("%d %f \n",i,h2->GetBinContent(ih,i));
186  float mean=h->GetMean();
187  float yield=h->Integral();
188  h->SetEntries(yield);
189  printf("******** work on %s mean=%.1f rms=%.1f yield=%.1f\n",txt1,mean,h->GetRMS(),yield);
190 
191  if(yield<cut_yield0) { hStat->SetBinContent(ih,kBad); continue; }
192  kBad<<=1;
193 
194  if(isSticky(h,mean)) { hStat->SetBinContent(ih,kBad); continue; }
195  kBad<<=1;
196 
197  // if(ih>=12367) return;
198  float adc1=mean-par_rms*par_nsig;
199  float adc2=mean+par_rms*par_nsig-1;
200  h->SetAxisRange( adc1,adc2);
201  float yield2=h->Integral();
202  float r1=yield2/yield;
203  printf(" range=%d,%d r1=%.3f yield2=%.1f\n",adc1,adc2,r1,yield2);
204 
205  if(r1< cut_yieldR1) { hStat->SetBinContent(ih,kBad); continue; }
206  kBad<<=1;
207 
208  h->Fit("gaus","Q","Rh", adc1, adc2);
209  TF1 *ff=h->GetFunction("gaus"); assert(ff);
210  ff->SetLineColor(kRed);
211  ff->SetLineWidth(1.);
212  h->Draw();
213 
214  float ped=ff->GetParameter(1);
215  float pedErr=ff->GetParError(1);
216  float sig=ff->GetParameter(2);
217  float chi2=ff->GetChisquare();
218  float ndf=ff->GetNDF();
219 
220  printf("chi2=%f ndf=%f\n", chi2,ndf);
221  if(chi2<cut_ch2min) { hStat->SetBinContent(ih,kBad); continue; }
222  // keep the same flag
223  hChi->SetBinContent(ih, chi2/ndf);
224  if(chi2/ndf> cut_ch2ndf) { hStat->SetBinContent(ih,kBad); continue; }
225  kBad<<=1;
226 
227  adc1=ped-sig*par_nsig;
228  adc2=ped+sig*par_nsig;
229  h->SetAxisRange( adc1,adc2);
230  float yield3=h->Integral();
231  float r2=yield2/yield;
232  printf("ped=%f sig=%f range=%d,%d r2=%.3f\n",ped,sig,adc1,adc2,r2);
233 
234  hSig->SetBinContent(ih,sig);
235  if(sig>cut_sigPed || sig<0) { hStat->SetBinContent(ih,kBad); continue; }
236  kBad<<=1;
237 
238  if(isRawAdc){
239  if(ped< cut_pedL) { hStat->SetBinContent(ih,kBad); continue; }
240  if(ped> cut_pedH) { hStat->SetBinContent(ih,kBad); continue; }
241  kBad<<=1;
242  }
243 
244  hPed->SetBinContent(ih,ped);
245  hPed->SetBinError(ih,pedErr);
246 
247  int np=grP->GetN();
248  grP->SetPoint(np,ih,ped);
249  grP->SetPointError(np,0,pedErr);
250 
251  hPeak->SetBinContent(ih,r2);
252  hStat->SetBinContent(ih,0); // good
253  }
254  h->SetAxisRange(y1,y2);
255 
256  printStat(id1,id2);
257  fd2->cd();
258  h->Write();
259 }
260 
261 //================
262 void isSticky(TH1F *h, float ped) {
263  assert(h);
264  axX=h->GetXaxis();
265  int kPed=axX->FindBin(ped);
266  printf("%s ped=%f kPed=%d\n",h->GetName(),ped,kPed);
267  float sum1=0,sum2=0;
268  for(int i=kPed-10;i<=kPed+10;i++) {
269  float val=h->GetBinContent(i);
270  // printf("i=%d val=%f\n",i,val);
271  if(i%2==0) sum1+=val;
272  if(i%2==1) sum2+=val;
273  }
274  float r=-1;
275  if(sum1>sum2) r=sum2/sum1;
276  else r=sum1/sum2;
277  printf("sum1=%f sum2=%f, r=%f\n",sum1,sum2,r);
278  if(r<0.1) return true; // is sticky bit
279  return false; // good spectrum
280 
281 }
282 
283 
284 //================
285 void printStat(int id1,int id2) {
286  int nBad=0, nGood=0;
287  printf("softID stat \n",);
288 
289  for(int i=id1; i<=id2;i++) {
290  int val=hStat->GetBinContent(i);
291  if(val<0.5) {
292  if(val==0)nGood++;
293  continue;
294  }
295  printf("%d 0x%0x\n",i,val);
296  // if(val>0x10) printf("softID=%d stat=0x%0x\n",i,val);
297  nBad++;
298  }
299  printf("\nstat summary: tot=%d good=%d bad=%d\n",id2-id1+1,nGood, nBad);
300 }
Generated by   doxygen 1.8.5