StRoot: StEEmcPool/StEzSmdCal/macros/fitPrePost.C Source File

00001 TCanvas *c;
00002 double par[6];  
00003 double *epar;  
00004 TFile *fdA[12];
00005 FILE *wfd=stdout; // output web page
00006 FILE *gfd;
00007 
00008 //======================================
00009 fitPrePost() {
00010   gStyle->SetStatW(0.22);
00011   gStyle->SetStatH(0.22);
00012   
00013   //  f=new TFile("iter4-pp/sect05/sum-sect5.hist.root");
00014 
00015   TH1F *h1=new TH1F("mpv","MPV  MIP gated; MPV of ADC-ped",35,-5,65); 
00016   TH1F *h2=new TH1F("mpvE","relative error of MPV , MIP gated; err(MPV)/MPV ",50,0,0.5); 
00017   TH1F *h3=new TH1F("mpvS","relative width of L-peak,  MIP gated; sigma/MPV",25,0.,1.); 
00018  TH2F *h4=new TH2F("mpv2","MPV from ; gated w/ MIP ; inclusive spectrum;",25,0,50,25,0,25); 
00019 
00020 //===sww===Set layer to analyze in next line
00021   char cT='R';
00022   openAll(cT);
00023   
00024   const float feta[]=
00025     {1.95,1.855,1.765,1.675,1.59,1.51,1.435,1.365,1.3,1.235,1.17,1.115};
00026   
00027   int sec=5;
00028   char core[100];
00029   
00030   int eta;
00031   char sub='C';
00032   
00033 //===sww===Set eta range and sector range in next two "for" statements
00034   for(eta=1;eta<=12;eta++) {
00035     int nErr=0, nOK=0;
00036     float mpvL=999, mpvH=0; 
00037     h1->Reset(); h2->Reset(); h3->Reset(); h4->Reset(); 
00038     fprintf(wfd," <tr> <th> %d <td> \n",eta); 
00039     gStyle->SetOptStat(1001111);
00040 
00041     for(sec=1; sec<=12;sec++) {
00042       TFile *f=fdA[sec-1];
00043       for(sub='A';sub<='E';sub++)     {
00044         sprintf(core,"%02d%c%c%02d",sec,cT,sub,eta);
00045         TString coreT=core;
00046         ha=(TH1F*)f->Get("a"+coreT);
00047         hd=(TH1F*)f->Get("d"+coreT);
00048         c=new TCanvas("aa","aa",400,400);
00049         plotOne(ha,hd);
00050         //return;
00051         c->Print(coreT+".ps");
00052         TString errS=QaOne(ha,hd,cT);
00053         printf("errS=%s=%d\n",errS.Data(),errS.Sizeof());
00054         bool isBad=errS.Sizeof()>1;
00055         if(isBad) { // report error channel
00056           fprintf(wfd,"     %s ,\n",(coreT+"-"+errS).Data());
00057           nErr++;
00058           //    continue;
00059         }
00060         if(errS.Contains("mask")) continue;
00061         nOK++;
00062         //    return;
00063         //  
00064         //  c->Print(coreT+".gif");
00065         float mpv=par[4];
00066         float mpvEr=epar[4];
00067         if(mpvL>mpv) mpvL=mpv;
00068         if(mpvH<mpv) mpvH=mpv;
00069         
00070         int ieta=eta-1; 
00071 
00072         float fac=TMath::TanH(feta[ieta])/0.0009; // assumed 0.9 MeV per plastic
00073         float err=sqrt(mpvEr*mpvEr+1);
00074         float gain=mpv*fac; 
00075         float sig=err*fac;
00076         fprintf(gfd,"%s %.0f %.0f %.1f %.1f \n",core,gain,sig,mpv,mpvEr);
00077       
00078         h1->Fill(mpv);
00079         if(mpv>0) {
00080           h2->Fill(mpvEr/mpv);
00081           h3->Fill(par[5]/mpv);
00082         }
00083         //  h4->Fill(mpv,mpvInc);
00084         
00085       } // end of eta bin
00086       
00087     }// end of sector loop
00088     fprintf(wfd," <td> %d <td> %d\n",nOK,nErr);
00089 
00090     char txt[100],sumN[100], pdfN[100];
00091     sprintf(pdfN,"%cfitEta%02d.pdf",cT,eta);
00092     sprintf(txt,"cat *%02d.ps | ps2pdf - %s",eta,pdfN);
00093     printf("%s\n",txt);
00094     system(txt);
00095     sprintf(txt,"mv %s /star/u/wissink/cal2006/tmp/",pdfN);
00096     printf("%s\n",txt);
00097     system(txt);
00098     fprintf(wfd,"     <td> <a href=\"%s\"> PDF </a>\n",pdfN);    
00099     fprintf(wfd,"     <td> %.1f to %.1f \n",mpvL,mpvH);
00100     
00101     gStyle->SetOptStat(1111111);
00102 
00103     sprintf(sumN,"mpv%c-eta%02d",cT,eta);
00104     c=new TCanvas(sumN,sumN,600,600);
00105     c->Divide(1,3);
00106     c->cd(1); h1->Draw();
00107     c->cd(2); h2->Draw();
00108     c->cd(3); h3->Draw();
00109     // c->cd(4); h4->Draw("box");
00110     c->Print();
00111     
00112     sprintf(txt,"ps2pdf %s.ps %s.pdf",sumN,sumN);
00113     printf("%s\n",txt);
00114     system(txt);
00115     sprintf(txt,"mv %s.pdf /star/u/wissink/cal2006/tmp/",sumN);
00116     printf("%s\n",txt);
00117     system(txt);
00118     fprintf(wfd,"     <td> <a href=\"%s.pdf\"> PDF </a>\n",sumN);
00119       
00120     
00121 }// end of loop over eta bins
00122   
00123   if(wfd!=stdout)fclose(wfd);
00124   fclose(gfd);
00125 }
00126 
00127 //=================================
00128 plotOne(TH1F *ha, TH1F *hd) {
00129   assert(ha);
00130   assert(hd);
00131 
00132   float xMax=220;
00133 
00134   hd->SetAxisRange(4,80);
00135   float sum=hd->Integral();
00136   float hdY=hd->GetMaximum();
00137   printf("%s sum=%f\n",hd->GetName(),sum);
00138   if(sum>210) 
00139     hd->Rebin(2);
00140   else 
00141     hd->Rebin(4);
00142 
00143   //exceptions:
00144   if(strstr(hd->GetName(),"d12PB05")) hd->Rebin(2);
00145   if(strstr(hd->GetName(),"d12RC05")) hd->Rebin(2);
00146 
00147   ha->SetAxisRange(5,50);
00148   int kb=ha->GetMaximumBin();
00149   float yMax=ha->GetBinContent(kb);
00150   ha->SetAxisRange(-10,xMax); 
00151   hd->SetAxisRange(-10,xMax);
00152   
00153   c->Divide(1,2);
00154 
00155   c->cd(1);
00156   ha->Draw();
00157   //gPad->SetGrid();
00158   ha->Fit("gaus","R","",-2,3);
00159   
00160  
00161   ha->SetMaximum(yMax*1.5);
00162  
00163 
00164   c->cd(2);
00165   hd->SetMaximum(2.5*hdY);
00166   fitGausPlusLand(hd);
00167   // gPad->SetGrid();
00168 
00169 }
00170 
00171 
00172 //---------------------------------------
00173 float fitGausPlusLand(TH1F *h) {
00174 
00175   float xMax=160.;
00176 
00177   gF=new TF1("gF","gaus",-2.,2.);
00178   gF->SetLineWidth(2);
00179   gF->SetLineColor(kBlue);
00180   gF->SetLineStyle(2);
00181 
00182   lF=new TF1("lF","landau",3,xMax);
00183   lF->SetLineWidth(2);
00184   lF->SetLineStyle(2);
00185   lF->SetLineColor(kRed);
00186 
00187   glF=new TF1("glF","gaus(0)+landau(3)",-2,xMax);
00188   glF->SetLineWidth(1);
00189   glF->SetLineColor(kGreen);
00190   glF->SetParNames("ampl-G","mean-G","sig-G","const-L","MPV-L","sig-L");
00191 
00192 
00193   h->Fit("gF","R");
00194   // return 1;
00195   h->Fit("lF","R0Q+");
00196   
00197   
00198   // copy starting point
00199   gF->GetParameters(par+0);
00200   lF->GetParameters(par+3);//return 1;
00201   glF->SetParameters(par);
00202   //freez position & width of pedestal residuum
00203   glF->FixParameter(1,par[1]);
00204   glF->FixParameter(2,par[2]);
00205 
00206 
00207   h->Fit("glF","R");
00208 
00209   // retrieve components
00210   glF->GetParameters(par);
00211   epar= glF->GetParErrors();
00212 
00213   gF->SetParameters(par+0);
00214   lF->SetParameters(par+3);
00215 
00216   gF->SetRange(-5,5);
00217   gF->Draw("same");
00218 
00219   lF->SetRange(0,xMax);
00220   lF->Draw("same");
00221 
00222   float yMax=lF->GetMaximum();
00223   printf("max=%f\n",yMax);
00224   return yMax;
00225 }
00226 
00227 //=================================
00228 TString QaOne(TH1F *ha, TH1F *hd,char cT ) {
00229   assert(ha);
00230   assert(hd);
00231 
00232   //............. pedestal
00233   ha->SetAxisRange(-10,10);
00234   int kb=ha->GetMaximumBin();
00235   float xb=ha->GetBinCenter(kb);
00236   //  printf("ped xAdc=%f\n",xb);
00237   if(fabs(xb)>1.) return "ped";
00238   if(ha->GetEntries() <=0) return "empty";
00239   // change limit to 0.9999 for R analysis - sww
00240   if(ha->Integral()/ha->GetEntries()>0.9999) return "mask";
00241 
00242 
00243   //.............Landau 
00244   //printf("xxx=%f\n",epar[4]/par[4]);
00245 
00246   switch (cT) {
00247   case 'R': // post shower
00248     if(par[4]<6. || par[4]>30.)  return "valL";
00249     break;
00250   default:
00251     if(par[4]<10. || par[4]>75.)  return "valL";
00252   }
00253 
00254   if(fabs(epar[4])<0.35)  return "erL"; // single bin pathology
00255   if(epar[4]/par[4]>0.5)  return "sigL";
00256   
00257   hd->SetAxisRange(5,50);
00258   float sum=hd->Integral();
00259   if(sum<20)  return "lSum";
00260   //  if(sum>500)  return "hSum";
00261   return "";
00262 }
00263 
00264 
00265 
00266 //========================
00267 void openAll(char cT) {
00268   int i;
00269   char txt[200];
00270   for(i=0;i<12;i++) {
00271     //  sprintf(txt,"/star/data05/scratch/balewski/2005-eemcCal/day171-hist/iter4-outA/sum-sect%d.hist.root",i+1);
00272     // sprintf(txt,"./iter12-mc/sum-sect%d.hist.root",i+1);
00273     // sprintf(txt,"/star/data05/scratch/balewski/2005-eemcCal/mc-hist/iter8-mc/sum-sect%d.hist.root",i+1);
00274     if(i<9) sprintf(txt,"iter4-pp/sect0%d/sum-sect%d.hist.root",i+1,i+1);
00275     if(i>8) sprintf(txt,"iter4-pp/sect%d/sum-sect%d.hist.root",i+1,i+1);
00276     fdA[i]=new TFile(txt);
00277     assert(fdA[i]->IsOpen());
00278   }
00279   sprintf(txt,"/star/u/wissink/cal2006/iter4-pp/gains%c-allSect.dat",cT); 
00280   gfd=fopen(txt,"w");
00281   assert(gfd);
00282   fprintf(gfd,"# final %c-gains from MIP using UxV\n",cT);
00283   fprintf(gfd,"# format: gain (ch/GeV), errGain, LandauFit: MPV(ADC) sigMPV\n");
00284 
00285   
00286   sprintf(txt,"/star/u/wissink/cal2006/iter4-pp/gains%c.html",cT);
00287   
00288   wfd=fopen(txt,"w");
00289   assert(wfd);
00290   
00291   fprintf(wfd,"MIP -- > <a href=\"gains%c-allSect.dat\"> final %c-gains</a> table \n",cT,cT);
00292   fprintf(wfd,"<table border=1>\n");
00293   fprintf(wfd,"<tr> <th> etaBin <th> problems <th> # of good<th> # of errors <th> spectra <br> (60 tiles)<th> MPV range <br>(ADC) <th> summary <br> plot\n");
00294   
00295   
00296 }