StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
fitPrePost.C
1 TCanvas *c;
2 double par[6];
3 double *epar;
4 TFile *fdA[12];
5 FILE *wfd=stdout; // output web page
6 FILE *gfd;
7 
8 //======================================
9 fitPrePost() {
10  gStyle->SetStatW(0.22);
11  gStyle->SetStatH(0.22);
12 
13  // f=new TFile("iter4-pp/sect05/sum-sect5.hist.root");
14 
15  TH1F *h1=new TH1F("mpv","MPV MIP gated; MPV of ADC-ped",35,-5,65);
16  TH1F *h2=new TH1F("mpvE","relative error of MPV , MIP gated; err(MPV)/MPV ",50,0,0.5);
17  TH1F *h3=new TH1F("mpvS","relative width of L-peak, MIP gated; sigma/MPV",25,0.,1.);
18  TH2F *h4=new TH2F("mpv2","MPV from ; gated w/ MIP ; inclusive spectrum;",25,0,50,25,0,25);
19 
20 //===sww===Set layer to analyze in next line
21  char cT='R';
22  openAll(cT);
23 
24  const float feta[]=
25  {1.95,1.855,1.765,1.675,1.59,1.51,1.435,1.365,1.3,1.235,1.17,1.115};
26 
27  int sec=5;
28  char core[100];
29 
30  int eta;
31  char sub='C';
32 
33 //===sww===Set eta range and sector range in next two "for" statements
34  for(eta=1;eta<=12;eta++) {
35  int nErr=0, nOK=0;
36  float mpvL=999, mpvH=0;
37  h1->Reset(); h2->Reset(); h3->Reset(); h4->Reset();
38  fprintf(wfd," <tr> <th> %d <td> \n",eta);
39  gStyle->SetOptStat(1001111);
40 
41  for(sec=1; sec<=12;sec++) {
42  TFile *f=fdA[sec-1];
43  for(sub='A';sub<='E';sub++) {
44  sprintf(core,"%02d%c%c%02d",sec,cT,sub,eta);
45  TString coreT=core;
46  ha=(TH1F*)f->Get("a"+coreT);
47  hd=(TH1F*)f->Get("d"+coreT);
48  c=new TCanvas("aa","aa",400,400);
49  plotOne(ha,hd);
50  //return;
51  c->Print(coreT+".ps");
52  TString errS=QaOne(ha,hd,cT);
53  printf("errS=%s=%d\n",errS.Data(),errS.Sizeof());
54  bool isBad=errS.Sizeof()>1;
55  if(isBad) { // report error channel
56  fprintf(wfd," %s ,\n",(coreT+"-"+errS).Data());
57  nErr++;
58  // continue;
59  }
60  if(errS.Contains("mask")) continue;
61  nOK++;
62  // return;
63  //
64  // c->Print(coreT+".gif");
65  float mpv=par[4];
66  float mpvEr=epar[4];
67  if(mpvL>mpv) mpvL=mpv;
68  if(mpvH<mpv) mpvH=mpv;
69 
70  int ieta=eta-1;
71 
72  float fac=TMath::TanH(feta[ieta])/0.0009; // assumed 0.9 MeV per plastic
73  float err=sqrt(mpvEr*mpvEr+1);
74  float gain=mpv*fac;
75  float sig=err*fac;
76  fprintf(gfd,"%s %.0f %.0f %.1f %.1f \n",core,gain,sig,mpv,mpvEr);
77 
78  h1->Fill(mpv);
79  if(mpv>0) {
80  h2->Fill(mpvEr/mpv);
81  h3->Fill(par[5]/mpv);
82  }
83  // h4->Fill(mpv,mpvInc);
84 
85  } // end of eta bin
86 
87  }// end of sector loop
88  fprintf(wfd," <td> %d <td> %d\n",nOK,nErr);
89 
90  char txt[100],sumN[100], pdfN[100];
91  sprintf(pdfN,"%cfitEta%02d.pdf",cT,eta);
92  sprintf(txt,"cat *%02d.ps | ps2pdf - %s",eta,pdfN);
93  printf("%s\n",txt);
94  system(txt);
95  sprintf(txt,"mv %s /star/u/wissink/cal2006/tmp/",pdfN);
96  printf("%s\n",txt);
97  system(txt);
98  fprintf(wfd," <td> <a href=\"%s\"> PDF </a>\n",pdfN);
99  fprintf(wfd," <td> %.1f to %.1f \n",mpvL,mpvH);
100 
101  gStyle->SetOptStat(1111111);
102 
103  sprintf(sumN,"mpv%c-eta%02d",cT,eta);
104  c=new TCanvas(sumN,sumN,600,600);
105  c->Divide(1,3);
106  c->cd(1); h1->Draw();
107  c->cd(2); h2->Draw();
108  c->cd(3); h3->Draw();
109  // c->cd(4); h4->Draw("box");
110  c->Print();
111 
112  sprintf(txt,"ps2pdf %s.ps %s.pdf",sumN,sumN);
113  printf("%s\n",txt);
114  system(txt);
115  sprintf(txt,"mv %s.pdf /star/u/wissink/cal2006/tmp/",sumN);
116  printf("%s\n",txt);
117  system(txt);
118  fprintf(wfd," <td> <a href=\"%s.pdf\"> PDF </a>\n",sumN);
119 
120 
121 }// end of loop over eta bins
122 
123  if(wfd!=stdout)fclose(wfd);
124  fclose(gfd);
125 }
126 
127 //=================================
128 plotOne(TH1F *ha, TH1F *hd) {
129  assert(ha);
130  assert(hd);
131 
132  float xMax=220;
133 
134  hd->SetAxisRange(4,80);
135  float sum=hd->Integral();
136  float hdY=hd->GetMaximum();
137  printf("%s sum=%f\n",hd->GetName(),sum);
138  if(sum>210)
139  hd->Rebin(2);
140  else
141  hd->Rebin(4);
142 
143  //exceptions:
144  if(strstr(hd->GetName(),"d12PB05")) hd->Rebin(2);
145  if(strstr(hd->GetName(),"d12RC05")) hd->Rebin(2);
146 
147  ha->SetAxisRange(5,50);
148  int kb=ha->GetMaximumBin();
149  float yMax=ha->GetBinContent(kb);
150  ha->SetAxisRange(-10,xMax);
151  hd->SetAxisRange(-10,xMax);
152 
153  c->Divide(1,2);
154 
155  c->cd(1);
156  ha->Draw();
157  //gPad->SetGrid();
158  ha->Fit("gaus","R","",-2,3);
159 
160 
161  ha->SetMaximum(yMax*1.5);
162 
163 
164  c->cd(2);
165  hd->SetMaximum(2.5*hdY);
166  fitGausPlusLand(hd);
167  // gPad->SetGrid();
168 
169 }
170 
171 
172 //---------------------------------------
173 float fitGausPlusLand(TH1F *h) {
174 
175  float xMax=160.;
176 
177  gF=new TF1("gF","gaus",-2.,2.);
178  gF->SetLineWidth(2);
179  gF->SetLineColor(kBlue);
180  gF->SetLineStyle(2);
181 
182  lF=new TF1("lF","landau",3,xMax);
183  lF->SetLineWidth(2);
184  lF->SetLineStyle(2);
185  lF->SetLineColor(kRed);
186 
187  glF=new TF1("glF","gaus(0)+landau(3)",-2,xMax);
188  glF->SetLineWidth(1);
189  glF->SetLineColor(kGreen);
190  glF->SetParNames("ampl-G","mean-G","sig-G","const-L","MPV-L","sig-L");
191 
192 
193  h->Fit("gF","R");
194  // return 1;
195  h->Fit("lF","R0Q+");
196 
197 
198  // copy starting point
199  gF->GetParameters(par+0);
200  lF->GetParameters(par+3);//return 1;
201  glF->SetParameters(par);
202  //freez position & width of pedestal residuum
203  glF->FixParameter(1,par[1]);
204  glF->FixParameter(2,par[2]);
205 
206 
207  h->Fit("glF","R");
208 
209  // retrieve components
210  glF->GetParameters(par);
211  epar= glF->GetParErrors();
212 
213  gF->SetParameters(par+0);
214  lF->SetParameters(par+3);
215 
216  gF->SetRange(-5,5);
217  gF->Draw("same");
218 
219  lF->SetRange(0,xMax);
220  lF->Draw("same");
221 
222  float yMax=lF->GetMaximum();
223  printf("max=%f\n",yMax);
224  return yMax;
225 }
226 
227 //=================================
228 TString QaOne(TH1F *ha, TH1F *hd,char cT ) {
229  assert(ha);
230  assert(hd);
231 
232  //............. pedestal
233  ha->SetAxisRange(-10,10);
234  int kb=ha->GetMaximumBin();
235  float xb=ha->GetBinCenter(kb);
236  // printf("ped xAdc=%f\n",xb);
237  if(fabs(xb)>1.) return "ped";
238  if(ha->GetEntries() <=0) return "empty";
239  // change limit to 0.9999 for R analysis - sww
240  if(ha->Integral()/ha->GetEntries()>0.9999) return "mask";
241 
242 
243  //.............Landau
244  //printf("xxx=%f\n",epar[4]/par[4]);
245 
246  switch (cT) {
247  case 'R': // post shower
248  if(par[4]<6. || par[4]>30.) return "valL";
249  break;
250  default:
251  if(par[4]<10. || par[4]>75.) return "valL";
252  }
253 
254  if(fabs(epar[4])<0.35) return "erL"; // single bin pathology
255  if(epar[4]/par[4]>0.5) return "sigL";
256 
257  hd->SetAxisRange(5,50);
258  float sum=hd->Integral();
259  if(sum<20) return "lSum";
260  // if(sum>500) return "hSum";
261  return "";
262 }
263 
264 
265 
266 //========================
267 void openAll(char cT) {
268  int i;
269  char txt[200];
270  for(i=0;i<12;i++) {
271  // sprintf(txt,"/star/data05/scratch/balewski/2005-eemcCal/day171-hist/iter4-outA/sum-sect%d.hist.root",i+1);
272  // sprintf(txt,"./iter12-mc/sum-sect%d.hist.root",i+1);
273  // sprintf(txt,"/star/data05/scratch/balewski/2005-eemcCal/mc-hist/iter8-mc/sum-sect%d.hist.root",i+1);
274  if(i<9) sprintf(txt,"iter4-pp/sect0%d/sum-sect%d.hist.root",i+1,i+1);
275  if(i>8) sprintf(txt,"iter4-pp/sect%d/sum-sect%d.hist.root",i+1,i+1);
276  fdA[i]=new TFile(txt);
277  assert(fdA[i]->IsOpen());
278  }
279  sprintf(txt,"/star/u/wissink/cal2006/iter4-pp/gains%c-allSect.dat",cT);
280  gfd=fopen(txt,"w");
281  assert(gfd);
282  fprintf(gfd,"# final %c-gains from MIP using UxV\n",cT);
283  fprintf(gfd,"# format: gain (ch/GeV), errGain, LandauFit: MPV(ADC) sigMPV\n");
284 
285 
286  sprintf(txt,"/star/u/wissink/cal2006/iter4-pp/gains%c.html",cT);
287 
288  wfd=fopen(txt,"w");
289  assert(wfd);
290 
291  fprintf(wfd,"MIP -- > <a href=\"gains%c-allSect.dat\"> final %c-gains</a> table \n",cT,cT);
292  fprintf(wfd,"<table border=1>\n");
293  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");
294 
295 
296 }
Generated by   doxygen 1.8.5