StRoot: StEEmcPool/StEzSmdCal/SmdGains.cxx Source File

00001 // $Id: SmdGains.cxx,v 1.11 2009/12/03 22:35:03 ogrebeny Exp $
00002  
00003 #include <assert.h>
00004 #include <stdlib.h>
00005 #include <string.h>
00006 
00007 #include <TClonesArray.h>
00008 #include <TObjArray.h> 
00009 #include <TString.h> 
00010 #include <TCanvas.h> 
00011 #include <TGraphErrors.h> 
00012 #include <TLine.h> 
00013 #include <TF1.h> 
00014 #include <TH1.h> 
00015 #include <TH2.h> 
00016 #include <TFile.h> 
00017 #include "TMath.h"
00018 #include "SmdGains.h"
00019 
00020 
00021 ClassImp(SmdGains)
00022 
00023 //--------------------------------------------------
00024 //--------------------------------------------------
00025 SmdGains::SmdGains(){
00026   HList=0;
00027   sectID=0;
00028   fdIn=0;
00029   c1=c2=0;
00030   gnCorFn=0;
00031   planeUV='X';
00032   memset(hA,0,sizeof(hA));
00033   memset(grA,0,sizeof(grA));
00034 
00035   idealMipEne=1.3; // (MeV) in 7mm of plastic
00036   // limits for expo fit
00037   // sww 1/17/07 - change adcMax from 100 -> 140
00038   adcMin=40;  adcMax=140; // SMD
00039    minSum=50;
00040   // sww 1/17/07 - change maxRelEr from 0.4 -> 0.5
00041   maxRelEr=0.5; // maximal relative error at any stage of calculations
00042   minMipEne=0.4;// (MeV) lower/upper thres for Landau Ene fit
00043   maxMipEne=4;
00044 }
00045 
00046 //--------------------------------------------------
00047 //--------------------------------------------------
00048 void SmdGains::init(){
00049   int i;
00050   assert(HList);
00051   
00052   char tt1[100];
00053   sprintf(tt1,"%02d%c",sectID,planeUV);
00054   plCore=tt1;
00055 
00056   //............... histos ............
00057   hA[0]=new TH1F("sum"+plCore,plCore+" Integral from raw spectra; strip ID; total counts",290,0.5,290.5);
00058    hA[1]=new TH1F("Lm"+plCore,plCore+" Mean of Landau fit to average MIP ene (2strip); MPV (MeV)",30,.5,2);
00059   hA[2]=new TH1F("Lw"+plCore,plCore+" Width of Landau fit to average MIP ene (2strip); Width=sigma (MeV)",20,0,1);
00060 
00061  
00062   for(i=0;i<mxH;i++) 
00063     if(hA[i]) HList->Add(hA[i]);
00064  
00065   //................ TGRaphs
00066   TGraphErrors*  gr=new TGraphErrors;
00067 
00068   grA[0]=0;
00069   
00070   gr=new TGraphErrors;
00071   gr->SetMarkerStyle(24);
00072   gr->SetName("mpvN"+plCore); // MPV of Landau fit
00073   gr->SetTitle(plCore+ " MPV of Landau fit to N 2-strips; strip ID; MIP energy (MeV)");
00074   grA[1]=gr;
00075 
00076 
00077   for(i=0;i<mxH;i++) 
00078     if(grA[i]) HList->Add(grA[i]);
00079   
00080   //............... other initializations ...........
00081   printf("cuts for %s : adcMin=%d ,adcMax=%d minSum=%d maxRelEr=%f\n",plCore.Data(),adcMin,adcMax,minSum, maxRelEr);
00082 
00083   for(i=0;i<mxS;i++) str[i].id=i+1;
00084 
00085 }
00086 
00087 //--------------------------------------------------
00088 //--------------------------------------------------
00089 TFile* SmdGains::open(TString fn) {
00090   fdIn=new TFile(fn);
00091   assert(fdIn->IsOpen());
00092   return fdIn;
00093 }
00094 
00095 
00096 
00097 //-------------------------------------------------
00098 //-------------------------------------------------
00099 void SmdGains::doGainCorr(int str1, int str2, int ns, int pl){
00100 
00101   TGraphErrors  *gr= grA[1];
00102   printf("doGainCorr() for %s, average over %d-strips pl=%d\n",gr->GetName(),ns,pl);
00103   TString nn=gr->GetName();
00104 
00105   if(pl==0) {
00106     c2=new TCanvas(nn,nn,300,400);
00107   } else {
00108     c2=new TCanvas(nn,nn,600,800);
00109   }
00110 
00111   c2->Divide(5,6);
00112   int i,k;
00113   for(i=str1,k=1; i<mxS; i+=ns,k++) {
00114     if(i>=str2) break;
00115     c2->cd(k);
00116      avrMipNEne(i,ns);
00117   }
00118   c2->cd(29);  hA[1]->Draw();
00119   c2->cd(30);  hA[2]->Draw();
00120 
00121     if( pl&1) c2->Print(nn+".ps");
00122   // if( pl&2) c2->Print(nn+".gif");
00123 
00124 }
00125 
00126 
00127 
00128 
00129 //-----------------------------------------
00130 //-----------------------------------------
00131 void SmdGains:: avrMipNEne( int str1,int ns){
00132   // sum MIP ene (from 2 strips) of ns strips
00133   char tit[100];
00134   assert(str1>0 && ns>0);
00135   int i;
00136 
00137   TH1F *hs=0;
00138   for(i=str1;i<str1+ns;i++) {
00139     if(i>mxS) break;
00140     sprintf(tit,"e%s%03d",plCore.Data(),i);
00141     TH1F* h= (TH1F*)fdIn->Get(tit); assert(h);
00142     HList->Add(h);
00143     //  c1->cd(k+1);
00144     //  h->Draw(); gPad->SetLogy(0);
00145     
00146     if(hs==0) {
00147       hs=(TH1F*) h->Clone();
00148     } else {
00149       hs->Add(h);
00150     }
00151   }
00152 
00153   sprintf(tit,"avr%s%03d",plCore.Data(),str1);
00154   hs->SetName(tit);
00155   TString nn=hs->GetTitle(); nn="avr"+nn;
00156   hs->SetTitle(nn);
00157   hs->Draw();
00158   HList->Add(hs);
00159   float sum=hs->Integral();
00160   hs->Rebin(3);
00161   if(sum<150)  hs->Rebin();
00162 
00163   // printf("%s --> sum=%.1f\n",hs->GetName(),sum );
00164   if(sum<50) return ;
00165 
00166   hs->Fit("landau","RQ","", minMipEne, maxMipEne);
00167   TF1* f=hs->GetFunction("landau");
00168   assert(f);
00169   f->SetLineColor(kRed);
00170   f->SetLineWidth(1);
00171   double *par=f->GetParameters();
00172   double *epar=f->GetParErrors();
00173   float   mpv=par[1];
00174   float   empv=epar[1];
00175   float rerr=0;
00176   
00177   hA[1]->Fill(mpv) ;
00178   hA[2]->Fill(par[2]) ;
00179 
00180   if(par[1]>0) rerr=empv/mpv;
00181   hs->SetAxisRange(0,5.);
00182  
00183   printf("%s MPV=%.2f +/- %.2f (%.1f%c) \n",hs->GetName(),mpv,empv,rerr*100.,37);
00184   
00185   if(rerr>maxRelEr)  {
00186     for(i=str1;i<str1+ns;i++) {
00187       if(i>mxS) break;    
00188       str[i-1].flag+=8;
00189     }
00190   }
00191 
00192   float x=str1+ns/2.; 
00193   TGraphErrors*  gr=grA[1];
00194   int n=gr->GetN();
00195   gr->SetPoint(n,x,mpv);
00196   gr->SetPointError(n,ns/2.,empv);
00197 
00198   return ;
00199 }
00200 
00201 
00202 //-----------------------------------------
00203 //-----------------------------------------
00204 void SmdGains:: plTGraph(const Char_t *shpFunc,int ig, int pl){
00205   // sum MIP ene (from 2 strips) of k strips
00206 
00207   TGraphErrors  *gr= grA[ig];
00208   assert(gr);
00209   printf("plot %s\n",gr->GetName());
00210   TString nn=gr->GetName();
00211   nn="C"+nn;  
00212   c2=new TCanvas(nn,nn,300,250);
00213   gr->Draw("AP");
00214   gr->Fit(shpFunc);
00215   gr->SetMinimum(0.7);
00216   gr->SetMaximum(1.8);
00217 
00218   gnCorFn=gr->GetFunction(shpFunc); assert(gnCorFn);
00219   gnCorFn->SetLineColor(kRed);
00220   gnCorFn ->SetLineWidth(2);
00221 
00222   if(ig==1 ) {
00223     TList *Lx;  TLine *ln;
00224     Lx=gr->GetListOfFunctions();    assert(Lx);
00225     ln=new TLine(1,idealMipEne, mxS,idealMipEne); ln->SetLineColor(kBlue); Lx->Add(ln);
00226   }
00227 
00228 
00229   if( pl&1) c2->Print(nn+".ps");
00230   if( pl&2) c2->Print(nn+".gif");
00231 
00232  }
00233   
00234 //-----------------------------------------
00235 //-----------------------------------------
00236 void SmdGains:: plFGC(){
00237   TString nn="fgc"+plCore;
00238   c2=new TCanvas(nn,nn,500,400);
00239 
00240   c2->Divide(2,2);
00241   c2->cd(1);  hA[1]->Draw();
00242   c2->cd(2);  hA[2]->Draw();
00243   c2->cd(3);  hA[3]->Draw();
00244   c2->cd(4);  hA[4]->Draw();
00245 
00246  }
00247   
00248 
00249 //-----------------------------------------
00250 //-----------------------------------------
00251 
00252 void SmdGains::fitSlopesSmd(int str1, int str2, int pl) {
00253   char tit[100];
00254   int ns=str2-str1+1;
00255   assert(str1>0 && ns>0);
00256   sprintf(tit,"%02d%c%03d+%d",sectID,planeUV,str1,ns);
00257   printf("fitSlopes() %s\n",tit);
00258   TString ttC=tit;
00259 
00260   TLine *ln0=new TLine(0,0,0,5000); // memory leak
00261   ln0->SetLineColor(kGreen);
00262 
00263   if(c1==0) c1=new TCanvas("aa1","aa1",800,700);// big
00264   c1->Clear();
00265   if(pl) {
00266     c1->Divide(5,6);
00267   } else {
00268     c1->Divide(2,2);
00269   }
00270   c1->SetName(ttC);
00271   c1->SetTitle(ttC);
00272   
00273   int k=0;
00274   int i;
00275   TList *Lx;  TLine *ln;
00276   for(i=str1;i<=str2;i++,k++) {
00277     if(k>mxS) break;
00278     StripG *s=&str[i-1];
00279     sprintf(tit,"a%s%03d",plCore.Data(),i);
00280     TH1F* h= (TH1F*)fdIn->Get(tit); assert(h);
00281     h->Rebin(2);
00282     HList->Add(h);
00283     c1->cd(k+1);
00284     h->SetAxisRange(adcMin,adcMax);
00285     float sum=h->Integral();
00286     h->SetAxisRange(-20,1.5*adcMax);
00287     //printf(" %s sum=%f\n",h->GetName(),sum);
00288     
00289     hA[0]->Fill(i,sum);
00290     s->sum1=(int)sum;
00291     h->Draw(); h->SetLineColor(kBlue);gPad->SetLogy();
00292     gPad->SetGridx(0);      gPad->SetGridy(0);
00293     float ymax=6000-str1*10;    h->SetMaximum(ymax);
00294     float ymin=100.9-str1/2.7;   
00295    // sww 1/17/07 change test from 30 ->300
00296     if(i<=300) ymin=0.9;
00297     h->SetMinimum(ymin);    
00298 
00299     Lx=h->GetListOfFunctions();    assert(Lx);
00300     ln=new TLine(adcMin,0,adcMin,30000); ln->SetLineColor(kMagenta);ln->SetLineStyle(2);
00301     Lx->Add(ln);
00302     ln=new TLine(adcMax,0,adcMax,30000);  ln->SetLineColor(kMagenta); ln->SetLineStyle(2);
00303     Lx->Add(ln);
00304     if(sum<minSum) { s->flag+=1; continue;}
00305     h->Fit("expo","RQ","",adcMin,adcMax);
00306 
00307     // ....exceptions in 2005 data, day 49
00308     if(strstr(h->GetName(),"a02V269")) h->Fit("expo","RQ","",60,130);
00309     if(strstr(h->GetName(),"a06U077")) h->Fit("expo","RQ","",80,150); // was wrong [800,150]
00310        //... end
00311     ln0->Draw();
00312     TF1* f=h->GetFunction("expo");
00313     f->SetLineColor(kRed);
00314     f->SetLineWidth(1);
00315     double *par=f->GetParameters();
00316     double *epar=f->GetParErrors();
00317 
00318     s->sl=par[1];
00319     s->esl=epar[1];
00320     if(epar[1]>maxRelEr *TMath::Abs(par[1])) {
00321       s->flag+=2;
00322     } 
00323   }
00324 
00325   if(pl&1) c1->Print(ttC+".ps");
00326   if(pl&2) c1->Print(ttC+".gif");
00327  
00328 }
00329 
00330 
00331 //-----------------------------------------
00332 //-----------------------------------------
00333 
00334 void SmdGains::fitSlopesTile(int eta1, int nEta, char cT, int pl) {
00335   char tit[100];
00336   assert(eta1>0 && nEta>0);
00337   sprintf(tit,"%02d%cA-E%02d+%d",sectID,cT,eta1,nEta);
00338   printf("fitSlopes() %s\n",tit);
00339   TString ttC=tit;
00340   float xLow=-50;
00341 
00342   if(cT=='T') { adcMin=15;  adcMax=45; } // towers
00343   if(c1==0) c1=new TCanvas("aa2","aa2",800,700);// big
00344 
00345   c1->Clear();
00346   if(pl)
00347     c1->Divide(5,6);
00348   else
00349     c1->Divide(2,2);
00350   c1->SetName(ttC);
00351   c1->SetTitle(ttC);
00352   
00353   int k=0;
00354   int i;
00355   TList *Lx;  TLine *ln;
00356   for(i=eta1;i<eta1+nEta;i++) {
00357     char sub='A';
00358     for(sub='A';sub<='E';sub++) {
00359       k++;
00360       sprintf(tit,"a%02d%c%c%02d",sectID,cT,sub,i); 
00361       TH1F* h= (TH1F*)fdIn->Get(tit); assert(h);
00362       // h->Rebin(2);
00363       HList->Add(h);
00364       c1->cd(k);
00365       h->SetAxisRange(adcMin,adcMax);
00366       float sum=h->Integral();
00367       h->SetAxisRange(xLow,1.5*adcMax);
00368       if(cT=='T') h->SetAxisRange(-20,100);
00369       // printf(" %s sum=%f\n",h->GetName(),sum);
00370       
00371       h->Draw(); h->SetLineColor(kBlue);
00372       gPad->SetLogy();
00373       gPad->SetGridx(0);      gPad->SetGridy(0);
00374       h->SetMinimum(.9);    h->SetMaximum(9000);
00375       //if(sum<4*minSum)  h->Rebin();
00376       
00377       Lx=h->GetListOfFunctions();    assert(Lx);
00378       ln=new TLine(adcMin,0,adcMin,30000); ln->SetLineColor(kMagenta); Lx->Add(ln);
00379       ln=new TLine(adcMax,0,adcMax,30000);  ln->SetLineColor(kMagenta);  Lx->Add(ln); 
00380       float sl=0,esl=0;
00381       if(sum>minSum) { //do fit
00382         h->Fit("expo","RQ","",adcMin,adcMax);
00383         TF1* f=h->GetFunction("expo");
00384         f->SetLineColor(kRed);
00385         f->SetLineWidth(1);
00386         double *par=f->GetParameters();
00387         double *epar=f->GetParErrors();
00388         sl=par[1];
00389         esl=epar[1];
00390       }
00391 
00392       printf("# %s %.4f %.4f %f\n",h->GetName(),sl,esl,sum);
00393 
00394     }
00395   }
00396 
00397   if(pl&1) c1->Print(ttC+".ps");
00398   if(pl&2) c1->Print(ttC+".gif");
00399  
00400 }
00401 
00402 
00403 
00404 //-------------------------------------------------
00405 //-------------------------------------------------
00406 void SmdGains:: doSlopesOnly(float fac){
00407   printf("working on gains from slopes only, fac=%f\n",fac);
00408 
00409   int i;
00410   // ...... relative gains
00411   for(i=0;i<mxS;i++) {
00412     StripG *s=str+i;
00413     if(s->sl<0){ // calculate is possible
00414       s->gc=-fac/s->sl;
00415       s->egc=s->esl*s->gc/TMath::Abs(s->sl);
00416     }
00417     if(s->egc>maxRelEr *s->gc) {
00418       s->flag+=4;
00419       s->gc=s->egc=0;
00420     }
00421 
00422   }
00423   
00424 }
00425 
00426 #if 0
00427 //-------------------------------------------------
00428 //-------------------------------------------------
00429 void SmdGains:: finish(int k) {
00430   int i;
00431   for(i=0;i<mxS;i++) {
00432     StripG *s=str+i;
00433     s->print();
00434   }
00435   
00436   
00437 }
00438 
00439 #endif
00440 
00441 //-------------------------------------------------
00442 //-------------------------------------------------
00443 void SmdGains::saveGains(FILE *fd) {
00444   if(fd==0) fd=stdout;
00445   int i,nBad=0;
00446   for(i=0;i<mxS;i++) {
00447     StripG *s=str+i;
00448     float rer=-1;
00449     if(s->gc>0) {
00450       rer=100.*s->egc/s->gc;
00451     } else {
00452       nBad++;
00453     }
00454     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);
00455   }
00456   fprintf(fd,"#nBad =%d\n",nBad);
00457 }
00458 
00459 //-------------------------------------------------
00460 //-------------------------------------------------
00461 void SmdGains::saveHisto(const Char_t *fname){
00462   TString outName;
00463   if(fname){
00464     outName=fname;
00465   } else {
00466     outName="smd"+plCore;
00467   }
00468   outName+=".hist.root";
00469   TFile f( outName,"recreate");
00470   assert(f.IsOpen());
00471   printf("%d histos are written  to '%s' ...\n",HList->GetEntries(),outName.Data());
00472   HList->Write();
00473   f.Close();
00474 }
00475 
00476 //=============================================
00477 //=============================================
00478 //=============================================
00479 //=============================================
00480 
00481 StripG::StripG(){
00482 clear();
00483 }
00484 
00485 //-------------------------
00486 void StripG::clear(){
00487   id=0;
00488   sl=esl=999;
00489   sum1=0; 
00490   gc=egc=0;
00491   flag=0;
00492   mpv1=empv1=0;
00493 }
00494 
00495 //-------------------------
00496 void StripG::print(){
00497   if(flag) printf("*");
00498   if (flag==0) {
00499     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);
00500   } else {
00501     printf("strip id=%d  --- sum1=%d\n",id,sum1);
00502   }
00503 }
00504 
00505 /*****************************************************************
00506  * $Log: SmdGains.cxx,v $
00507  * Revision 1.11  2009/12/03 22:35:03  ogrebeny
00508  * Fixed compiler warnings, mostly char* -> const char*
00509  *
00510  * Revision 1.10  2009/01/26 14:37:42  fisyak
00511  * Add missing (in ROOT 5.22) includes
00512  *
00513  * Revision 1.9  2007/08/21 13:10:04  balewski
00514  * final, used in 2006 offline calibration by soScott
00515  *
00516  *
00517  * VS: ----------------------------------------------------------------------
00518  *
00519  * Revision 1.7  2005/09/29 13:57:57  balewski
00520  * after SMD gains were rescaled
00521  *
00522  * Revision 1.6  2005/08/09 18:46:31  balewski
00523  * after smd calib in 2005
00524  *
00525  * Revision 1.5  2004/11/02 21:29:08  balewski
00526  * ange hist range
00527  *
00528  * Revision 1.4  2004/10/08 14:34:50  balewski
00529  * as used for PQRUV calib for pp200, 2004
00530  *
00531  * Revision 1.3  2004/09/22 00:45:52  balewski
00532  * ready for calib of smd
00533  *
00534  * Revision 1.2  2004/09/14 19:38:43  balewski
00535  * new version, SMD calib is now too complicated
00536  *
00537  * Revision 1.1  2004/09/11 04:57:34  balewski
00538  * cleanup
00539  *
00540  *
00541  *
00542  ********************************************************************/
00543 
00544