StRoot: StEmcPool/StBarrelCalib2008/macros/doAvrPed.C Source File

00001 /* notes, R9067013 OLD
00002    BPRS:
00003    kBad=0x1:  yield in ADC range [100,300] 100 counts, input was ~220 events
00004           7 tiles:  softID=3301-4,3322-4
00005    kBad=0x2: yiled @ [mean +/- 9 ADC] below 70% 
00006           0 tiles
00007    kBad=0x4: bad chi2/dof >10          
00008           42-7=35 tiles, softID=3863 ...3940, about 50%
00009    kBad=x8 = ped low & high, set wide [105,295]
00010           0 tiles
00011    kBad=0x10  sticky low bit @ pedestal area
00012           85-42=43 , tiles softID=3862 ...3940, about 50%
00013    kBad=0x20  sig ped >2.5 ADC
00014           132-85=47 tiles, softID=3021-3032, 3041-3048, 3052,3103,3621-3652
00015 
00016 
00017 */
00018 
00019 const int  mxTP=2, mxBTiles=4800;
00020 char  cTile[mxTP]={'T','P'};
00021 char  *cTile4[mxTP]={"BTOW","BPRS"};
00022 TH1F * hPed, *sigP, *hChi,*hSig;
00023 TH1I * hStat;
00024 TGraphErrors *grP;
00025 TFile *fd2=0, *fd1=0;
00026 TH1F *hPeak; // yield of pedestal
00027 TCanvas *c=0;
00028 int isRawAdc=0; // 1=raw, 0=pedSub
00029 
00030 doAvrPed( int k=10,int run=9067013) {
00031   gStyle->SetPalette(1,0);
00032   gStyle->SetOptFit(1);
00033   gStyle->SetOptStat(1001110);
00034 
00035   char *pathIn="outA3/";
00036   char *pathOut="iter4/";
00037   char txt[1000], txt2[1000];
00038 
00039   sprintf(txt,"%s/R%dp.barCal.-1.hist.root",pathIn,run);
00040 
00041   fd1=new TFile(txt);  assert(fd1->IsOpen());
00042   //  fd1.ls(); return;
00043   sprintf(txt,"%s/pedBprsR%davr.hist.root",pathOut,run);
00044   fd2=new TFile(txt,"recreate"); 
00045   assert(fd2->IsOpen());
00046   int i;
00047    grP=new TGraphErrors;
00048    //  gr->SetMarkerStyle(20+m);
00049    grP->SetMarkerSize(0.5);
00050 
00051   TH2F * h2T= (TH2F *)fd1->Get("BTOW_c0"); assert(h2T);
00052   TH2F * h2P= (TH2F *)fd1->Get("BPRS_c0"); assert(h2P);
00053   h2P->GetXaxis()->SetLabelSize(0.08);
00054   h2P->GetYaxis()->SetLabelSize(0.08);
00055   h2P->Draw("colz");
00056   //  Printf("XXX %f\n", h2P->GetYaxis()->SetTitleSize());
00057 
00058   // return;
00059  int id1=1, id2=id1+k-1;
00060   if(k<0) { id1=1; id2=mxBTiles; }
00061   c=new TCanvas("aa","aa",400,300);
00062   // id1=id2=k;//tmp
00063   fitPedBT(1,run,h2P,id1,id2 );
00064   gPad->SetLogy();
00065 
00066   // save output histo
00067   fd2->cd();  hPed->Write();
00068   hStat->Write();
00069   hChi->Write();
00070   hSig->Write();
00071   hPeak->Write();
00072   
00073   // return;
00074   c=new TCanvas(fd1->GetName(),fd1->GetName(),1000,800);
00075   c->Divide(1,5);
00076   gStyle->SetOptStat(10);
00077   c->cd(1);  h2P->Draw("colz");  h2P->SetAxisRange(id1,id2); 
00078   hPed->Draw("same");
00079 
00080   c->cd(2); hStat->Draw(); hStat->SetAxisRange(id1,id2); 
00081    hStat->SetMinimum(.5);
00082    if(hStat->GetMaximum()>10) gPad->SetLogy();
00083    gPad->SetGrid();
00084  
00085    c->cd(3);
00086    hChi->Draw();  hChi->SetAxisRange(id1,id2);// hChi->SetMaximum(20);
00087    gPad->SetGrid();
00088    
00089   c->cd(4);    
00090   hSig->Draw();  hSig->SetAxisRange(id1,id2); hSig->SetMaximum(5);
00091   gPad->SetGrid();
00092 
00093   c->cd(5);
00094   hPeak->Draw();  hPeak->SetAxisRange(id1,id2); hPeak->SetMinimum(0.7);
00095   gPad->SetGrid();
00096 
00097   //  gPad->SetLogy();
00098   // c->cd(1); hPed->Draw(); hPed->SetAxisRange(id1,id2); 
00099   
00100 
00101   return;
00102   fd2->cd();
00103   hres->Write();
00104   grPedR[0]->Write();
00105 
00106 }
00107 
00108 
00109 
00110 //==============================================
00111 //==============================================
00112 //==============================================
00113 void  fitPedBT( int itp,int run,TH2F *h2, int id1, int id2) {
00114   if(id2>18000) id2=18000;
00115     char txt[1000], txt2[1000];
00116     
00117     sprintf(txt,"ped%s",cTile4[itp]);
00118     sprintf(txt2,"ped %s R%d; %s soft ID; pedestal +/- fit err (ADC)",cTile4[itp],run,cTile4[itp]);
00119     hPed=new TH1F(txt,txt2,mxBTiles,0.5,mxBTiles+0.5);
00120     hPed->Sumw2();
00121     hPed->GetXaxis()->SetLabelSize(0.08);
00122     hPed->GetYaxis()->SetLabelSize(0.08);
00123 
00124     sprintf(txt,"stat%s",cTile4[itp]);
00125     sprintf(txt2,"status %s R%d; %s soft ID; jan status",cTile4[itp],run,cTile4[itp]);
00126     hStat=new TH1I(txt,txt2,mxBTiles,0.5,mxBTiles+0.5); hStat->Reset(-1); // default is bad
00127     hStat->GetXaxis()->SetLabelSize(0.08);
00128     hStat->GetYaxis()->SetLabelSize(0.08);
00129 
00130     sprintf(txt,"chi%s",cTile4[itp]);
00131     sprintf(txt2,"chi2/DOF %s R%d; %s soft ID; ped Chi2/DOF",cTile4[itp],run,cTile4[itp]);
00132     hChi=new TH1F(txt,txt2,mxBTiles,0.5,mxBTiles+0.5); hStat->Reset(-1); // default is bad
00133     hChi->GetXaxis()->SetLabelSize(0.08);
00134     hChi->GetYaxis()->SetLabelSize(0.08);
00135     hChi->GetYaxis()->SetTitleSize(0.2);
00136 
00137     sprintf(txt,"sigPed%s",cTile4[itp]);
00138     sprintf(txt2,"sigma(ped) %s R%d; %s soft ID; sig(ped) ADC",cTile4[itp],run,cTile4[itp]);
00139     hSig=new TH1F(txt,txt2,mxBTiles,0.5,mxBTiles+0.5); hStat->Reset(-1); // default is bad
00140     hSig->GetXaxis()->SetLabelSize(0.08);
00141     hSig->GetYaxis()->SetLabelSize(0.08);
00142 
00143     hPeak=new TH1F("hpedPeak", "integral of pedestal peak;soft ID", mxBTiles,0.5,mxBTiles+0.5);
00144     hPeak->GetXaxis()->SetLabelSize(0.08);
00145     hPeak->GetYaxis()->SetLabelSize(0.08);
00146 
00147   
00148   float par_nsig=3; // integration range for QA
00149   float par_rms=2; // to approximate fit range
00150   float cut_yield0=100; 
00151   float cut_yieldR1=0.7;
00152   float cut_ch2min=0.1; 
00153   float cut_ch2ndf=200.; 
00154   float cut_pedL=105;
00155   float cut_pedH=295;
00156   float cut_pedH=295;
00157   float cut_sigPed=2.7;
00158 
00159   //  float cut_minYiled=50;
00160 
00161   axX=h2->GetXaxis();
00162   float x1=axX->GetXmin();
00163   float x2=axX->GetXmax();
00164   int nbX=axX->GetNbins();
00165   printf("X-axis range  --> [%.1f, %.1f], nb=%d %s\n",x1,x2,nbX,axX->GetTitle());
00166 
00167   axY=h2->GetYaxis();
00168   float y1=axY->GetXmin();
00169   float y2=axY->GetXmax();
00170   int nbY=axY->GetNbins();
00171   printf("Y-axis range  --> [%.1f, %.1f], nb=%d\n",y1,y2,nbY);
00172 
00173   TH1F*h=new TH1F("h1","h1",nbY,y1,y2); // working histo for 1-D spectrum
00174   // do projections
00175   int ih;
00176   for(ih=id1; ih<=id2; ih++) {
00177     char txt1[100], txt2[1000];
00178     sprintf(txt1,"id%d",ih);
00179     sprintf(txt2,"%s soft id=%d;%s ",cTile4[itp],ih,axY->GetTitle());
00180     h->SetTitle(txt2);
00181     h->Reset();  h->SetAxisRange(y1,y2);
00182     int i;
00183     int kBad=1;
00184     for(i=1;i<=nbY;i++) h->SetBinContent(i,h2->GetBinContent(ih,i));
00185     //for(i=1;i<=nbY;i++) printf("%d %f \n",i,h2->GetBinContent(ih,i));
00186     float mean=h->GetMean();
00187     float yield=h->Integral();
00188     h->SetEntries(yield);
00189     printf("******** work on %s  mean=%.1f rms=%.1f yield=%.1f\n",txt1,mean,h->GetRMS(),yield);
00190 
00191     if(yield<cut_yield0) {  hStat->SetBinContent(ih,kBad); continue; }
00192     kBad<<=1;
00193 
00194     if(isSticky(h,mean)) {  hStat->SetBinContent(ih,kBad); continue; }
00195     kBad<<=1;
00196 
00197     // if(ih>=12367) return;
00198     float adc1=mean-par_rms*par_nsig;
00199     float adc2=mean+par_rms*par_nsig-1;
00200     h->SetAxisRange( adc1,adc2);
00201     float yield2=h->Integral();
00202     float r1=yield2/yield;
00203     printf(" range=%d,%d  r1=%.3f yield2=%.1f\n",adc1,adc2,r1,yield2);
00204     
00205     if(r1< cut_yieldR1) {  hStat->SetBinContent(ih,kBad); continue; }
00206     kBad<<=1;
00207 
00208     h->Fit("gaus","Q","Rh", adc1, adc2);
00209     TF1 *ff=h->GetFunction("gaus"); assert(ff);
00210     ff->SetLineColor(kRed);
00211     ff->SetLineWidth(1.);
00212     h->Draw();
00213     
00214     float ped=ff->GetParameter(1);
00215     float pedErr=ff->GetParError(1);
00216     float sig=ff->GetParameter(2);
00217     float chi2=ff->GetChisquare();
00218     float ndf=ff->GetNDF();
00219   
00220     printf("chi2=%f ndf=%f\n", chi2,ndf);
00221     if(chi2<cut_ch2min)  {  hStat->SetBinContent(ih,kBad); continue; }
00222     // keep the same flag
00223     hChi->SetBinContent(ih, chi2/ndf);
00224     if(chi2/ndf> cut_ch2ndf)  {  hStat->SetBinContent(ih,kBad); continue; }
00225     kBad<<=1;
00226  
00227     adc1=ped-sig*par_nsig;
00228     adc2=ped+sig*par_nsig;
00229     h->SetAxisRange( adc1,adc2);
00230     float yield3=h->Integral();
00231     float r2=yield2/yield;
00232     printf("ped=%f sig=%f range=%d,%d r2=%.3f\n",ped,sig,adc1,adc2,r2);
00233 
00234     hSig->SetBinContent(ih,sig);
00235     if(sig>cut_sigPed || sig<0) {  hStat->SetBinContent(ih,kBad); continue; }
00236     kBad<<=1;
00237 
00238     if(isRawAdc){
00239       if(ped< cut_pedL)  {  hStat->SetBinContent(ih,kBad); continue; }
00240       if(ped> cut_pedH)  {  hStat->SetBinContent(ih,kBad); continue; }
00241       kBad<<=1;
00242     }
00243      
00244     hPed->SetBinContent(ih,ped);
00245     hPed->SetBinError(ih,pedErr);
00246 
00247     int np=grP->GetN();
00248     grP->SetPoint(np,ih,ped);
00249     grP->SetPointError(np,0,pedErr);
00250     
00251     hPeak->SetBinContent(ih,r2);
00252     hStat->SetBinContent(ih,0); // good  
00253   }
00254   h->SetAxisRange(y1,y2);
00255   
00256   printStat(id1,id2);
00257   fd2->cd();
00258   h->Write();
00259 }
00260 
00261 //================
00262 void isSticky(TH1F *h, float ped) {
00263   assert(h);
00264   axX=h->GetXaxis();
00265   int kPed=axX->FindBin(ped);
00266   printf("%s ped=%f kPed=%d\n",h->GetName(),ped,kPed);
00267   float sum1=0,sum2=0;
00268   for(int i=kPed-10;i<=kPed+10;i++) {
00269     float val=h->GetBinContent(i);
00270     // printf("i=%d val=%f\n",i,val);
00271     if(i%2==0) sum1+=val;
00272     if(i%2==1) sum2+=val;
00273   }
00274   float r=-1;
00275   if(sum1>sum2) r=sum2/sum1;
00276   else r=sum1/sum2;
00277   printf("sum1=%f sum2=%f, r=%f\n",sum1,sum2,r);
00278   if(r<0.1) return true; // is sticky bit
00279   return false; // good spectrum
00280 
00281 }
00282 
00283 
00284 //================
00285 void printStat(int id1,int id2) {
00286   int nBad=0, nGood=0;
00287   printf("softID stat \n",);
00288   
00289   for(int i=id1; i<=id2;i++) {
00290     int val=hStat->GetBinContent(i);
00291     if(val<0.5) { 
00292       if(val==0)nGood++; 
00293       continue;
00294     }
00295     printf("%d 0x%0x\n",i,val);
00296   //  if(val>0x10)  printf("softID=%d stat=0x%0x\n",i,val);
00297     nBad++;
00298   }
00299   printf("\nstat summary: tot=%d good=%d bad=%d\n",id2-id1+1,nGood, nBad);
00300 }