StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
SmdGains.cxx
1 // $Id: SmdGains.cxx,v 1.11 2009/12/03 22:35:03 ogrebeny Exp $
2 
3 #include <assert.h>
4 #include <stdlib.h>
5 #include <string.h>
6 
7 #include <TClonesArray.h>
8 #include <TObjArray.h>
9 #include <TString.h>
10 #include <TCanvas.h>
11 #include <TGraphErrors.h>
12 #include <TLine.h>
13 #include <TF1.h>
14 #include <TH1.h>
15 #include <TH2.h>
16 #include <TFile.h>
17 #include "TMath.h"
18 #include "SmdGains.h"
19 
20 
21 ClassImp(SmdGains)
22 
23 //--------------------------------------------------
24 //--------------------------------------------------
25 SmdGains::SmdGains(){
26  HList=0;
27  sectID=0;
28  fdIn=0;
29  c1=c2=0;
30  gnCorFn=0;
31  planeUV='X';
32  memset(hA,0,sizeof(hA));
33  memset(grA,0,sizeof(grA));
34 
35  idealMipEne=1.3; // (MeV) in 7mm of plastic
36  // limits for expo fit
37  // sww 1/17/07 - change adcMax from 100 -> 140
38  adcMin=40; adcMax=140; // SMD
39  minSum=50;
40  // sww 1/17/07 - change maxRelEr from 0.4 -> 0.5
41  maxRelEr=0.5; // maximal relative error at any stage of calculations
42  minMipEne=0.4;// (MeV) lower/upper thres for Landau Ene fit
43  maxMipEne=4;
44 }
45 
46 //--------------------------------------------------
47 //--------------------------------------------------
48 void SmdGains::init(){
49  int i;
50  assert(HList);
51 
52  char tt1[100];
53  sprintf(tt1,"%02d%c",sectID,planeUV);
54  plCore=tt1;
55 
56  //............... histos ............
57  hA[0]=new TH1F("sum"+plCore,plCore+" Integral from raw spectra; strip ID; total counts",290,0.5,290.5);
58  hA[1]=new TH1F("Lm"+plCore,plCore+" Mean of Landau fit to average MIP ene (2strip); MPV (MeV)",30,.5,2);
59  hA[2]=new TH1F("Lw"+plCore,plCore+" Width of Landau fit to average MIP ene (2strip); Width=sigma (MeV)",20,0,1);
60 
61 
62  for(i=0;i<mxH;i++)
63  if(hA[i]) HList->Add(hA[i]);
64 
65  //................ TGRaphs
66  TGraphErrors* gr=new TGraphErrors;
67 
68  grA[0]=0;
69 
70  gr=new TGraphErrors;
71  gr->SetMarkerStyle(24);
72  gr->SetName("mpvN"+plCore); // MPV of Landau fit
73  gr->SetTitle(plCore+ " MPV of Landau fit to N 2-strips; strip ID; MIP energy (MeV)");
74  grA[1]=gr;
75 
76 
77  for(i=0;i<mxH;i++)
78  if(grA[i]) HList->Add(grA[i]);
79 
80  //............... other initializations ...........
81  printf("cuts for %s : adcMin=%d ,adcMax=%d minSum=%d maxRelEr=%f\n",plCore.Data(),adcMin,adcMax,minSum, maxRelEr);
82 
83  for(i=0;i<mxS;i++) str[i].id=i+1;
84 
85 }
86 
87 //--------------------------------------------------
88 //--------------------------------------------------
89 TFile* SmdGains::open(TString fn) {
90  fdIn=new TFile(fn);
91  assert(fdIn->IsOpen());
92  return fdIn;
93 }
94 
95 
96 
97 //-------------------------------------------------
98 //-------------------------------------------------
99 void SmdGains::doGainCorr(int str1, int str2, int ns, int pl){
100 
101  TGraphErrors *gr= grA[1];
102  printf("doGainCorr() for %s, average over %d-strips pl=%d\n",gr->GetName(),ns,pl);
103  TString nn=gr->GetName();
104 
105  if(pl==0) {
106  c2=new TCanvas(nn,nn,300,400);
107  } else {
108  c2=new TCanvas(nn,nn,600,800);
109  }
110 
111  c2->Divide(5,6);
112  int i,k;
113  for(i=str1,k=1; i<mxS; i+=ns,k++) {
114  if(i>=str2) break;
115  c2->cd(k);
116  avrMipNEne(i,ns);
117  }
118  c2->cd(29); hA[1]->Draw();
119  c2->cd(30); hA[2]->Draw();
120 
121  if( pl&1) c2->Print(nn+".ps");
122  // if( pl&2) c2->Print(nn+".gif");
123 
124 }
125 
126 
127 
128 
129 //-----------------------------------------
130 //-----------------------------------------
131 void SmdGains:: avrMipNEne( int str1,int ns){
132  // sum MIP ene (from 2 strips) of ns strips
133  char tit[100];
134  assert(str1>0 && ns>0);
135  int i;
136 
137  TH1F *hs=0;
138  for(i=str1;i<str1+ns;i++) {
139  if(i>mxS) break;
140  sprintf(tit,"e%s%03d",plCore.Data(),i);
141  TH1F* h= (TH1F*)fdIn->Get(tit); assert(h);
142  HList->Add(h);
143  // c1->cd(k+1);
144  // h->Draw(); gPad->SetLogy(0);
145 
146  if(hs==0) {
147  hs=(TH1F*) h->Clone();
148  } else {
149  hs->Add(h);
150  }
151  }
152 
153  sprintf(tit,"avr%s%03d",plCore.Data(),str1);
154  hs->SetName(tit);
155  TString nn=hs->GetTitle(); nn="avr"+nn;
156  hs->SetTitle(nn);
157  hs->Draw();
158  HList->Add(hs);
159  float sum=hs->Integral();
160  hs->Rebin(3);
161  if(sum<150) hs->Rebin();
162 
163  // printf("%s --> sum=%.1f\n",hs->GetName(),sum );
164  if(sum<50) return ;
165 
166  hs->Fit("landau","RQ","", minMipEne, maxMipEne);
167  TF1* f=hs->GetFunction("landau");
168  assert(f);
169  f->SetLineColor(kRed);
170  f->SetLineWidth(1);
171  double *par=f->GetParameters();
172  const double *epar=f->GetParErrors();
173  float mpv=par[1];
174  float empv=epar[1];
175  float rerr=0;
176 
177  hA[1]->Fill(mpv) ;
178  hA[2]->Fill(par[2]) ;
179 
180  if(par[1]>0) rerr=empv/mpv;
181  hs->SetAxisRange(0,5.);
182 
183  printf("%s MPV=%.2f +/- %.2f (%.1f%c) \n",hs->GetName(),mpv,empv,rerr*100.,37);
184 
185  if(rerr>maxRelEr) {
186  for(i=str1;i<str1+ns;i++) {
187  if(i>mxS) break;
188  str[i-1].flag+=8;
189  }
190  }
191 
192  float x=str1+ns/2.;
193  TGraphErrors* gr=grA[1];
194  int n=gr->GetN();
195  gr->SetPoint(n,x,mpv);
196  gr->SetPointError(n,ns/2.,empv);
197 
198  return ;
199 }
200 
201 
202 //-----------------------------------------
203 //-----------------------------------------
204 void SmdGains:: plTGraph(const Char_t *shpFunc,int ig, int pl){
205  // sum MIP ene (from 2 strips) of k strips
206 
207  TGraphErrors *gr= grA[ig];
208  assert(gr);
209  printf("plot %s\n",gr->GetName());
210  TString nn=gr->GetName();
211  nn="C"+nn;
212  c2=new TCanvas(nn,nn,300,250);
213  gr->Draw("AP");
214  gr->Fit(shpFunc);
215  gr->SetMinimum(0.7);
216  gr->SetMaximum(1.8);
217 
218  gnCorFn=gr->GetFunction(shpFunc); assert(gnCorFn);
219  gnCorFn->SetLineColor(kRed);
220  gnCorFn ->SetLineWidth(2);
221 
222  if(ig==1 ) {
223  TList *Lx; TLine *ln;
224  Lx=gr->GetListOfFunctions(); assert(Lx);
225  ln=new TLine(1,idealMipEne, mxS,idealMipEne); ln->SetLineColor(kBlue); Lx->Add(ln);
226  }
227 
228 
229  if( pl&1) c2->Print(nn+".ps");
230  if( pl&2) c2->Print(nn+".gif");
231 
232  }
233 
234 //-----------------------------------------
235 //-----------------------------------------
236 void SmdGains:: plFGC(){
237  TString nn="fgc"+plCore;
238  c2=new TCanvas(nn,nn,500,400);
239 
240  c2->Divide(2,2);
241  c2->cd(1); hA[1]->Draw();
242  c2->cd(2); hA[2]->Draw();
243  c2->cd(3); hA[3]->Draw();
244  c2->cd(4); hA[4]->Draw();
245 
246  }
247 
248 
249 //-----------------------------------------
250 //-----------------------------------------
251 
252 void SmdGains::fitSlopesSmd(int str1, int str2, int pl) {
253  char tit[100];
254  int ns=str2-str1+1;
255  assert(str1>0 && ns>0);
256  sprintf(tit,"%02d%c%03d+%d",sectID,planeUV,str1,ns);
257  printf("fitSlopes() %s\n",tit);
258  TString ttC=tit;
259 
260  TLine *ln0=new TLine(0,0,0,5000); // memory leak
261  ln0->SetLineColor(kGreen);
262 
263  if(c1==0) c1=new TCanvas("aa1","aa1",800,700);// big
264  c1->Clear();
265  if(pl) {
266  c1->Divide(5,6);
267  } else {
268  c1->Divide(2,2);
269  }
270  c1->SetName(ttC);
271  c1->SetTitle(ttC);
272 
273  int k=0;
274  int i;
275  TList *Lx; TLine *ln;
276  for(i=str1;i<=str2;i++,k++) {
277  if(k>mxS) break;
278  StripG *s=&str[i-1];
279  sprintf(tit,"a%s%03d",plCore.Data(),i);
280  TH1F* h= (TH1F*)fdIn->Get(tit); assert(h);
281  h->Rebin(2);
282  HList->Add(h);
283  c1->cd(k+1);
284  h->SetAxisRange(adcMin,adcMax);
285  float sum=h->Integral();
286  h->SetAxisRange(-20,1.5*adcMax);
287  //printf(" %s sum=%f\n",h->GetName(),sum);
288 
289  hA[0]->Fill(i,sum);
290  s->sum1=(int)sum;
291  h->Draw(); h->SetLineColor(kBlue);gPad->SetLogy();
292  gPad->SetGridx(0); gPad->SetGridy(0);
293  float ymax=6000-str1*10; h->SetMaximum(ymax);
294  float ymin=100.9-str1/2.7;
295  // sww 1/17/07 change test from 30 ->300
296  if(i<=300) ymin=0.9;
297  h->SetMinimum(ymin);
298 
299  Lx=h->GetListOfFunctions(); assert(Lx);
300  ln=new TLine(adcMin,0,adcMin,30000); ln->SetLineColor(kMagenta);ln->SetLineStyle(2);
301  Lx->Add(ln);
302  ln=new TLine(adcMax,0,adcMax,30000); ln->SetLineColor(kMagenta); ln->SetLineStyle(2);
303  Lx->Add(ln);
304  if(sum<minSum) { s->flag+=1; continue;}
305  h->Fit("expo","RQ","",adcMin,adcMax);
306 
307  // ....exceptions in 2005 data, day 49
308  if(strstr(h->GetName(),"a02V269")) h->Fit("expo","RQ","",60,130);
309  if(strstr(h->GetName(),"a06U077")) h->Fit("expo","RQ","",80,150); // was wrong [800,150]
310  //... end
311  ln0->Draw();
312  TF1* f=h->GetFunction("expo");
313  f->SetLineColor(kRed);
314  f->SetLineWidth(1);
315  double *par=f->GetParameters();
316  const double *epar=f->GetParErrors();
317 
318  s->sl=par[1];
319  s->esl=epar[1];
320  if(epar[1]>maxRelEr *TMath::Abs(par[1])) {
321  s->flag+=2;
322  }
323  }
324 
325  if(pl&1) c1->Print(ttC+".ps");
326  if(pl&2) c1->Print(ttC+".gif");
327 
328 }
329 
330 
331 //-----------------------------------------
332 //-----------------------------------------
333 
334 void SmdGains::fitSlopesTile(int eta1, int nEta, char cT, int pl) {
335  char tit[100];
336  assert(eta1>0 && nEta>0);
337  sprintf(tit,"%02d%cA-E%02d+%d",sectID,cT,eta1,nEta);
338  printf("fitSlopes() %s\n",tit);
339  TString ttC=tit;
340  float xLow=-50;
341 
342  if(cT=='T') { adcMin=15; adcMax=45; } // towers
343  if(c1==0) c1=new TCanvas("aa2","aa2",800,700);// big
344 
345  c1->Clear();
346  if(pl)
347  c1->Divide(5,6);
348  else
349  c1->Divide(2,2);
350  c1->SetName(ttC);
351  c1->SetTitle(ttC);
352 
353  int k=0;
354  int i;
355  TList *Lx; TLine *ln;
356  for(i=eta1;i<eta1+nEta;i++) {
357  char sub='A';
358  for(sub='A';sub<='E';sub++) {
359  k++;
360  sprintf(tit,"a%02d%c%c%02d",sectID,cT,sub,i);
361  TH1F* h= (TH1F*)fdIn->Get(tit); assert(h);
362  // h->Rebin(2);
363  HList->Add(h);
364  c1->cd(k);
365  h->SetAxisRange(adcMin,adcMax);
366  float sum=h->Integral();
367  h->SetAxisRange(xLow,1.5*adcMax);
368  if(cT=='T') h->SetAxisRange(-20,100);
369  // printf(" %s sum=%f\n",h->GetName(),sum);
370 
371  h->Draw(); h->SetLineColor(kBlue);
372  gPad->SetLogy();
373  gPad->SetGridx(0); gPad->SetGridy(0);
374  h->SetMinimum(.9); h->SetMaximum(9000);
375  //if(sum<4*minSum) h->Rebin();
376 
377  Lx=h->GetListOfFunctions(); assert(Lx);
378  ln=new TLine(adcMin,0,adcMin,30000); ln->SetLineColor(kMagenta); Lx->Add(ln);
379  ln=new TLine(adcMax,0,adcMax,30000); ln->SetLineColor(kMagenta); Lx->Add(ln);
380  float sl=0,esl=0;
381  if(sum>minSum) { //do fit
382  h->Fit("expo","RQ","",adcMin,adcMax);
383  TF1* f=h->GetFunction("expo");
384  f->SetLineColor(kRed);
385  f->SetLineWidth(1);
386  double *par=f->GetParameters();
387  const double *epar=f->GetParErrors();
388  sl=par[1];
389  esl=epar[1];
390  }
391 
392  printf("# %s %.4f %.4f %f\n",h->GetName(),sl,esl,sum);
393 
394  }
395  }
396 
397  if(pl&1) c1->Print(ttC+".ps");
398  if(pl&2) c1->Print(ttC+".gif");
399 
400 }
401 
402 
403 
404 //-------------------------------------------------
405 //-------------------------------------------------
406 void SmdGains:: doSlopesOnly(float fac){
407  printf("working on gains from slopes only, fac=%f\n",fac);
408 
409  int i;
410  // ...... relative gains
411  for(i=0;i<mxS;i++) {
412  StripG *s=str+i;
413  if(s->sl<0){ // calculate is possible
414  s->gc=-fac/s->sl;
415  s->egc=s->esl*s->gc/TMath::Abs(s->sl);
416  }
417  if(s->egc>maxRelEr *s->gc) {
418  s->flag+=4;
419  s->gc=s->egc=0;
420  }
421 
422  }
423 
424 }
425 
426 #if 0
427 //-------------------------------------------------
428 //-------------------------------------------------
429 void SmdGains:: finish(int k) {
430  int i;
431  for(i=0;i<mxS;i++) {
432  StripG *s=str+i;
433  s->print();
434  }
435 
436 
437 }
438 
439 #endif
440 
441 //-------------------------------------------------
442 //-------------------------------------------------
443 void SmdGains::saveGains(FILE *fd) {
444  if(fd==0) fd=stdout;
445  int i,nBad=0;
446  for(i=0;i<mxS;i++) {
447  StripG *s=str+i;
448  float rer=-1;
449  if(s->gc>0) {
450  rer=100.*s->egc/s->gc;
451  } else {
452  nBad++;
453  }
454  fprintf(fd,"%s%03d %.1f %.1f (%.1f%c) sum1=%d flag=0x%0x mpv1=%.1f %.1f\n",plCore.Data(),s->id,s->gc,s->egc,rer,37,s->sum1,s->flag,s->mpv1,s->empv1);
455  }
456  fprintf(fd,"#nBad =%d\n",nBad);
457 }
458 
459 //-------------------------------------------------
460 //-------------------------------------------------
461 void SmdGains::saveHisto(const Char_t *fname){
462  TString outName;
463  if(fname){
464  outName=fname;
465  } else {
466  outName="smd"+plCore;
467  }
468  outName+=".hist.root";
469  TFile f( outName,"recreate");
470  assert(f.IsOpen());
471  printf("%d histos are written to '%s' ...\n",HList->GetEntries(),outName.Data());
472  HList->Write();
473  f.Close();
474 }
475 
476 //=============================================
477 //=============================================
478 //=============================================
479 //=============================================
480 
481 StripG::StripG(){
482 clear();
483 }
484 
485 //-------------------------
486 void StripG::clear(){
487  id=0;
488  sl=esl=999;
489  sum1=0;
490  gc=egc=0;
491  flag=0;
492  mpv1=empv1=0;
493 }
494 
495 //-------------------------
496 void StripG::print(){
497  if(flag) printf("*");
498  if (flag==0) {
499  printf("strip id=%d sl=%f esl=%f gc=%f egc=%f sum1=%d mpv1=%f empv1=%f\n",id,sl,esl,gc,egc,sum1,mpv1,empv1);
500  } else {
501  printf("strip id=%d --- sum1=%d\n",id,sum1);
502  }
503 }
504 
505 /*****************************************************************
506  * $Log: SmdGains.cxx,v $
507  * Revision 1.11 2009/12/03 22:35:03 ogrebeny
508  * Fixed compiler warnings, mostly char* -> const char*
509  *
510  * Revision 1.10 2009/01/26 14:37:42 fisyak
511  * Add missing (in ROOT 5.22) includes
512  *
513  * Revision 1.9 2007/08/21 13:10:04 balewski
514  * final, used in 2006 offline calibration by soScott
515  *
516  *
517  * VS: ----------------------------------------------------------------------
518  *
519  * Revision 1.7 2005/09/29 13:57:57 balewski
520  * after SMD gains were rescaled
521  *
522  * Revision 1.6 2005/08/09 18:46:31 balewski
523  * after smd calib in 2005
524  *
525  * Revision 1.5 2004/11/02 21:29:08 balewski
526  * ange hist range
527  *
528  * Revision 1.4 2004/10/08 14:34:50 balewski
529  * as used for PQRUV calib for pp200, 2004
530  *
531  * Revision 1.3 2004/09/22 00:45:52 balewski
532  * ready for calib of smd
533  *
534  * Revision 1.2 2004/09/14 19:38:43 balewski
535  * new version, SMD calib is now too complicated
536  *
537  * Revision 1.1 2004/09/11 04:57:34 balewski
538  * cleanup
539  *
540  *
541  *
542  ********************************************************************/
543 
544 
Generated by   doxygen 1.8.5