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