StRoot: StEEmcPool/StEEmcPointMaker/EEmcSectorFit.cxx Source File

00001 
00015 #include <cassert>
00016 
00017 #include "EEmcSectorFit.h"
00018 #include "eeSinglePeak.h"
00019 #include "TH1F.h"
00020 #include "TCanvas.h"
00021 #include "TF1.h"
00022 #include "TString.h"
00023 #include <iostream>
00024 #include <algorithm>
00025 #include "TLine.h"
00026 #include "TList.h"
00027 ClassImp(EEmcSectorFit);
00028 
00029 // ----------------------------------------------------------------------------
00030 // default constructor has enough parameters for 10 gammas per plane
00031 EEmcSectorFit::EEmcSectorFit(Int_t max) : TMinuit(4*max)
00032 {
00033   fLwarn=false; 
00034   SetPrintLevel(-1); 
00035   Command("SET ERR 1");
00036   mSMD[0]=0;
00037   mSMD[1]=0;
00038   doPermutations = true;
00039   mNDF=0;
00040   mChi2=0.; 
00041 }
00042 
00043 
00044 EEmcSectorFit::~EEmcSectorFit()
00045 {
00046 }
00047 
00048 // ----------------------------------------------------------------------------
00049 Double_t EEmcSectorFit::FitFunc( Double_t x, Int_t plane ) const
00050 {
00051   // loop over all gammas
00052   Double_t sum = 0.;
00053   for ( UInt_t i=0;i<yield.size();i++ )
00054     {
00055       Double_t a=yield[i];
00056       Double_t s=sigma[i];
00057       Double_t m=(plane==0)?umean[i]:vmean[i];
00058       Double_t X[]={x,0.,0.};
00059       Double_t P[]={a,m,s,0.2,3.5}; // norm, mean, sigma, frac2, rel width
00060       sum += eeSinglePeak(X,P);
00061     }
00062   return sum;
00063 }
00064 // ----------------------------------------------------------------------------
00065 Int_t EEmcSectorFit::Eval(Int_t np,Double_t* gr,Double_t& chi2,Double_t* par,Int_t flg)
00066 {
00067  
00068   // initialize chi2 to zero
00069   Double_t mychi2 = 0.;
00070   mChi2 = 0.;
00071   mNDF  = 0; 
00072 
00073   // pass paramerters to the arrays containing parameters
00074   for ( UInt_t i=0;i<yield.size();i++ )
00075     {
00076       yield[i] = par[4*i+0];
00077       sigma[i] = par[4*i+1];
00078       umean[i] = par[4*i+2];
00079       vmean[i] = par[4*i+3];
00080     }
00081 
00082   // loop over both planes
00083   for ( Int_t plane=0;plane<2; plane++ )
00084     {
00085 
00086       // get histogram for this plane
00087       TH1F *histo = mSMD[plane];
00088       // loop over all smd strips
00089       for ( Int_t strip=0;strip<288;strip++ )
00090         {
00091           // get n mips in this plane/strip
00092           Double_t nmip = histo -> GetBinContent(strip+1);
00093           // evaluate the fit function
00094           Double_t x    = (Double_t)strip + 0.5;
00095           Double_t f    = FitFunc(x, plane);
00096           // add to chi2
00097           if ( nmip > 0. ) {
00098             mychi2 += (nmip-f)*(nmip-f)/nmip;
00099             mNDF++; 
00100           }
00101           // kludge to count 0's near mean against chi2
00102           // need to go to a max likelyhood method
00103           else if ( f > 0.499 ) {
00104             mNDF++; 
00105             mychi2 += (nmip-f)*(nmip-f)/f;
00106           }
00107         }
00108     }
00109 
00110   chi2 = mychi2;
00111   mChi2 = mychi2;
00112   mNDF -= 4*(int)yield.size(); 
00113 
00114   return 0;
00115 
00116 }
00117 // ----------------------------------------------------------------------------
00118 Double_t EEmcSectorFit::Residual( Int_t x, Int_t plane ) const
00119 {
00120   // zero residual if x is w/in +/- 2 strips of another gamma
00121   for ( UInt_t i=0;i<yield.size();i++ )
00122     {
00123       Double_t xx=(plane==0)?umean[i]:vmean[i];
00124       if ( TMath::Abs(xx-(Double_t)x-0.5) <= 1.0 ) return 0.;
00125     }
00126   TH1F *histo=mSMD[plane];
00127   Int_t bin = x + 1;  
00128   Double_t r = histo->GetBinContent(bin) - FitFunc( (Double_t)x+0.5, plane );
00129   return r;
00130 }
00131 // ----------------------------------------------------------------------------
00132 Double_t EEmcSectorFit::Residual( Int_t x, Int_t plane, Int_t dx, Int_t side ) const
00133 {
00134 
00135 
00136   Int_t min;//=TMath::Max(x-dx,0);
00137   Int_t max;//=TMath::Min(x+dx,288);
00138 
00139   
00140   min=TMath::Max(x-dx,0);
00141   max=TMath::Min(x+dx,288);
00142   
00143   if ( side==1 ) 
00144     {
00145       min=TMath::Max(x,0);
00146     }
00147   else if ( side==2 )
00148     {
00149       max=TMath::Min(x,288);
00150     }
00151 
00152     
00153   Double_t r=0.; 
00154   for ( Int_t xx=min;xx<=max;xx++ )
00155     {
00156       r += Residual( xx, plane );
00157     }
00158   return r;
00159 } 
00160 // ----------------------------------------------------------------------------
00161 Int_t EEmcSectorFit::MaxStrip(Int_t plane) const
00162 {
00163   Double_t max=0.;
00164   Int_t imax=-1;
00165   for ( Int_t i=0;i<288;i++ )
00166     {
00167       Double_t r=Residual(i,plane);
00168       if ( r>max )
00169         {
00170           max=r;
00171           imax=i;
00172         }
00173     }
00174   return imax;
00175 }
00176 // ----------------------------------------------------------------------------
00177 void EEmcSectorFit::AddCandidate(Double_t y, Double_t s, Double_t u, Double_t v)
00178 {
00179   //  std::cout << "Adding candidate" << std::endl;
00180   yield.push_back(y);
00181   sigma.push_back(s);
00182   umean.push_back(u);
00183   vmean.push_back(v);
00184 #if 0
00185   TLine *lU = new TLine(u,0.,u,y);lU->SetLineColor(2);
00186   TLine *lV = new TLine(v,0.,v,y);lV->SetLineColor(2);
00187   mSMD[0]->GetListOfFunctions()->Add(lU);
00188   mSMD[1]->GetListOfFunctions()->Add(lV);
00189 #endif
00191   InitParameters();
00192 }
00193 // ----------------------------------------------------------------------------
00194 // initialzes parameters for the N-gamma fit
00195 void EEmcSectorFit::InitParameters()
00196 {
00197   Int_t flg=0;
00198   Command("CLEAR");
00199   Command("SET ERR 1");
00200   for ( Int_t i=0;i<(Int_t)yield.size();i++ )
00201     {
00202       mnparm(4*i+0, TString("yield"), yield[i], 0.1, 0.25*yield[i], 2.0*yield[i], flg );
00203       mnparm(4*i+1, TString("sigma"), sigma[i], 0.1, 0.5, 2.5,          flg );
00204       mnparm(4*i+2, TString("umean"), umean[i], 0.1, 0.0, 288.0,        flg );
00205       mnparm(4*i+3, TString("vmean"), vmean[i], 0.1, 0.0, 288.0,        flg );
00206     }
00207 }
00208 // ----------------------------------------------------------------------------
00209 void EEmcSectorFit::Draw(Option_t *opts)
00210 {
00211   //  TCanvas *canvas = new TCanvas("c1","c1",800,400);
00212   /*
00213   TPad *canvas=gPad;
00214   if ( !gPad ) {
00215     gPad=(TPad *)new TCanvas("c1","c1",500,500);
00216     gPad->Divide(2);
00217     canvas=gPad;
00218   }
00219   */
00220   //  canvas->Divide(2);
00221 
00222   /*
00223   mSMD[0]->GetListOfFunctions()->Delete();
00224   mSMD[1]->GetListOfFunctions()->Delete();
00225   */
00226 
00227   assert(mSMD[0]);
00228   assert(mSMD[1]);
00229 
00230   Double_t umin=999.;
00231   Double_t vmin=888.;
00232   Double_t umax=-1.;
00233   Double_t vmax=-1.;
00234 
00236   for ( UInt_t i=0;i<yield.size();i++ )
00237     {
00238       TString uname="fitu";uname+=i;      
00239       TString vname="fitv";vname+=i;
00240       TF1 *fitu=new TF1(uname,eeSinglePeak,0.,288.,5);
00241       TF1 *fitv=new TF1(vname,eeSinglePeak,0.,288.,5);
00242       fitu->SetLineColor(4);
00243       fitv->SetLineColor(4);
00244       fitu->SetLineWidth(1);
00245       fitv->SetLineWidth(1);
00246       fitu->SetParameter(0, yield[i] );
00247       fitu->SetParameter(1, umean[i] );
00248       fitu->SetParameter(2, sigma[i] );
00249       fitu->SetParameter(3, 0.2 );
00250       fitu->SetParameter(4, 3.5 );
00251       fitv->SetParameter(0, yield[i] );
00252       fitv->SetParameter(1, vmean[i] );
00253       fitv->SetParameter(2, sigma[i] );
00254       fitv->SetParameter(3, 0.2 );
00255       fitv->SetParameter(4, 3.5 );
00256       mSMD[0]->GetListOfFunctions()->Add( fitu );
00257       mSMD[1]->GetListOfFunctions()->Add( fitv );
00258       umin=TMath::Min( fitu->GetParameter(1), umin );
00259       umax=TMath::Max( fitu->GetParameter(1), umax );
00260       vmin=TMath::Min( fitv->GetParameter(1), vmin );
00261       vmax=TMath::Max( fitv->GetParameter(1), vmax );      
00262     }
00263 
00264   umin=TMath::Max(umin,1.0);
00265   umax=TMath::Min(umax,287.0);
00266   vmin=TMath::Max(vmin,1.0);
00267   vmax=TMath::Min(vmax,287.0);
00268 
00269   if ( umin < umax )
00270     mSMD[0]->GetXaxis()->SetRangeUser(umin-20.,umax+20.);
00271   if ( vmin < vmax )
00272     mSMD[1]->GetXaxis()->SetRangeUser(vmin-20.,vmax+20.);
00273 
00274 
00275   //  canvas->cd(1);
00276   //  mSMD[0]->Draw(opts);
00277   //  canvas->cd(2);
00278   //  mSMD[1]->Draw(opts);
00279 
00280 }
00281 // ----------------------------------------------------------------------------
00282 void EEmcSectorFit::TryPermutations()
00283 {
00284 
00285   if ( !doPermutations ) return;
00286   
00287   std::vector<Double_t> chi2vec;
00288 
00290   std::vector< std::vector<Double_t> > yields;
00291   std::vector< std::vector<Double_t> > sigmas;
00292   std::vector< std::vector<Double_t> > umeans;
00293   std::vector< std::vector<Double_t> > vmeans;
00294   std::vector< Double_t > mychi2;
00295   
00296   std::vector< Double_t > old_yield = yield;
00297   std::vector< Double_t > old_sigma = sigma;
00298   std::vector< Double_t > old_umean = umean;
00299   std::vector< Double_t > old_vmean = vmean;
00300 
00302   while ( std::next_permutation( old_vmean.begin(), old_vmean.end() ) )
00303     {
00304 
00306       vmean = old_vmean;
00307       
00309       InitParameters();
00310 
00312       for ( UInt_t i=0;i<sigma.size();i++ ) {
00313         FixParameter(4*i+1);
00314       }
00316       Migrad();
00318       mychi2.push_back( mChi2 );
00319       yields.push_back( yield );
00320       sigmas.push_back( sigma );
00321       umeans.push_back( umean );
00322       vmeans.push_back( vmean );
00324       yield = old_yield;
00325       sigma = old_sigma;
00326       umean = old_umean;
00327 
00328     }
00329 
00331     Double_t min = 9.0E12;
00332     for ( UInt_t i=0;i<mychi2.size();i++ )
00333       {
00334         if ( mychi2[i]< min )
00335           {
00336             min=mychi2[i];
00337             yield=yields[i];
00338             sigma=sigmas[i];
00339             umean=umeans[i];
00340             vmean=vmeans[i];
00341           }
00342       }
00343 
00345     InitParameters();
00347     for ( UInt_t i=0;i<sigma.size();i++ ) {
00348       Release(4*i+1);
00349     }
00350     Migrad();
00351     
00352 }
00353 // ----------------------------------------------------------------------------
00354 void EEmcSectorFit::print() const
00355 {
00356     std::cout << "doPermutations=" << doPermutations << std::endl;
00357     std::cout << "chi^2         =" << mChi2 << std::endl; 
00358     std::cout << "ndf           =" << mNDF << std::endl; 
00359 }
00360 // ----------------------------------------------------------------------------
00361 void EEmcSectorFit::Clear(Option_t *opts)
00362 {
00363   yield.clear();
00364   sigma.clear();
00365   umean.clear();
00366   vmean.clear();
00367   Command("CLEAR");
00368 }
00369 // ----------------------------------------------------------------------------
00370 void EEmcSectorFit::AddFits( TH1F *uu, TH1F *vv )
00371 {
00372 
00373   assert(uu);
00374   assert(vv);
00375 
00376   Double_t umin=999.;
00377   Double_t vmin=888.;
00378   Double_t umax=-1.;
00379   Double_t vmax=-1.;
00380 
00382   for ( UInt_t i=0;i<yield.size();i++ )
00383     {
00384       TString uname="fitu";uname+=i;      
00385       TString vname="fitv";vname+=i;
00386       TF1 *fitu=new TF1(uname,eeSinglePeak,0.,288.,5);
00387       TF1 *fitv=new TF1(vname,eeSinglePeak,0.,288.,5);
00388       fitu->SetLineColor(4);
00389       fitv->SetLineColor(4);
00390       fitu->SetLineWidth(1);
00391       fitv->SetLineWidth(1);
00392       fitu->SetParameter(0, yield[i] );
00393       fitu->SetParameter(1, umean[i] );
00394       fitu->SetParameter(2, sigma[i] );
00395       fitu->SetParameter(3, 0.2 );
00396       fitu->SetParameter(4, 3.5 );
00397       fitv->SetParameter(0, yield[i] );
00398       fitv->SetParameter(1, vmean[i] );
00399       fitv->SetParameter(2, sigma[i] );
00400       fitv->SetParameter(3, 0.2 );
00401       fitv->SetParameter(4, 3.5 );
00402       uu->GetListOfFunctions()->Add( fitu );
00403       vv->GetListOfFunctions()->Add( fitv );
00404       umin=TMath::Min( fitu->GetParameter(1), umin );
00405       umax=TMath::Max( fitu->GetParameter(1), umax );
00406       vmin=TMath::Min( fitv->GetParameter(1), vmin );
00407       vmax=TMath::Max( fitv->GetParameter(1), vmax );      
00408     }
00409 
00410   umin=TMath::Max(umin,1.0);
00411   umax=TMath::Min(umax,287.0);
00412   vmin=TMath::Max(vmin,1.0);
00413   vmax=TMath::Min(vmax,287.0);
00414 
00415   if ( umin < umax )
00416     uu->GetXaxis()->SetRangeUser(umin-20.,umax+20.);
00417   if ( vmin < vmax )
00418     vv->GetXaxis()->SetRangeUser(vmin-20.,vmax+20.);
00419 
00420 
00421 }