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

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