Project:/TrackingHL/src/CotNegLogT0Fitter.cc

CDF KITS
[ source navigation ]
[ diff markup ]
[ identifier search ]
[ freetext search ]
[ file search ]
Architecture:	[ i386 ]
Version:	[ 4.10.4 ] [ 4.10.5 ] [ 4.8.4 ] [ 4.8.4l3s ] [ 4.8.5 ] [ 4.9.0 ] [ 4.9.1 ] [ 4.9.1.hpt3 ] [ 4.9.1hpt3 ] [ 4.9.1top1 ] [ 5.0.0 ] [ 5.1.0 ] [ 5.1.0beamonly ] [ 5.1.1 ] [ 5.2.0 ] [ 5.3.0 ] [ 5.3.1 ] [ 5.3.1dsp ] [ 5.3.3 ] [ 5.3.3_nt ] [ 5.3.4 ] [ 6.1.1 ] [ 6.1.1b ] [ 6.1.2 ] [ 6.1.3 ] [ 6.1.4 ] [ 6.1.4int3 ] [ 6.1.4mc ] [ 6.1.4mc_a ] [ 6.1.6 ] [ development ]
001 #include "TrackingHL/Utility/CotNegLogT0Fitter.hh"
002 #include "TMinuit.h"
003 #include <iostream>
004 #include <cassert>
005 #include <cerrno>
006 #include <limits>
007 #include <iomanip>
008 
009 
010 const double CotNegLogT0Fitter::TZERO_LIMIT = 10.0;
011 
012 
013 CotNegLogT0Fitter::CotNegLogT0Fitter() {
014   _input = 0;
015   _internal = 0;
016 }
017 
018 CotNegLogT0Fitter::~CotNegLogT0Fitter(){
019 }
020 
021 
022 CotT0Fitter::FitStatus 
023 CotNegLogT0Fitter::calculateT0 (double & tzero, double & etzero, 
024                                 CotT0Fitter::t0Input* input, 
025                                 CotT0Fitter::t0Internal* internal)
026 { 
027   if (input->trk_t0.size()==0)
028   {
029     if ( internal )
030     {
031       internal->ntrk  = 0;
032       internal->ndof  = 0;
033       internal->chi2  = -999;
034       internal->ncall = 0;
035     }
036     return UNDEFINED;
037   }
038 
039 // std::cout << "Max p = " << p << "  max prob = " << maxProb << std::endl
040 //           << "T0 guess:  " << tzero << " assuming PID " << maxPid << " "
041 //           << "for match " << mbest << std::endl;
042 
043   // Use production time with minimum neg log likelihood as seed t0.
044   //
045   double t0seed;
046   double minLLF = std::numeric_limits<double>::max();
047   for(int mi = 0; mi < input->trk_t0.size(); mi++) {
048           for (int m=PION; m<NPID; m++) {       // skip electrons and muons...
049                   if (m==PION) {              
050                           tzero = input->trk_t0[mi] ;
051                   } else if (m==KAON) {       
052                           tzero = input->trk_kt0[mi] ;
053                   } else if (m==PROTON) {             
054                           tzero = input->trk_pt0[mi] ;
055                   }
056                   double llf = getLikelihood( tzero, input );
057                   if ( llf < minLLF ) {
058                           minLLF = llf;
059                           t0seed = tzero;
060                   }
061           }
062         }
063 
064   tzero = t0seed;
065   
066   //then do the pid likelihood fit, without skipping any matches:
067   if ( internal )
068   {
069     internal->ncall = 0;
070   }
071   _input    = input;
072   _internal = internal;
073   
074   int status = do_pid_likelihood_fit(tzero,etzero);
075 
076   FitStatus fitStatus;
077   if ( status == 0 )
078   { 
079     // Check that fitted t0 is not close to the limits of
080     // the allowed range of t0 values. (Close means within
081     // about 3 sigma in order to give the fitter room to 
082     // operate properly.)
083     //
084     if ( fabs( tzero - TZERO_LIMIT ) < 0.2  
085          || 
086          fabs( tzero + TZERO_LIMIT ) < 0.2 )
087     {
088       fitStatus = CotT0Fitter::FIT_AT_LIMIT;
089     }
090     else
091     {
092       fitStatus = CotT0Fitter::OK;
093     }
094   }
095   else
096   {
097     fitStatus = CotT0Fitter::FIT_ERROR;
098   }
099 
100   // Fill internal information for the final state of the fit
101   // Sets ntrk, ndof, chi2 = likelihood at tzero
102   //
103   if ( internal )
104   {
105     getLikelihood( tzero, input, internal );
106   }
107 
108   return fitStatus;
109 }
110 
111 int CotNegLogT0Fitter::do_pid_likelihood_fit(double & tzero, double & etzero)
112 {
113   static TMinuit minuit(1);
114   double arglis[10];
115   int ierr = 0;
116   instance = this;
117 
118   minuit.SetFCN ( CotNegLogT0Fitter::pid_likelihood_fcn );
119   minuit.SetPrintLevel(-1);
120 
121   arglis[0] = 1;
122   minuit.mnexcm("SET SETNOWARNINGS",arglis,1,ierr);
123   minuit.mncler();
124   minuit.mnparm(0, "t0_mean", tzero, 0.01, -TZERO_LIMIT, TZERO_LIMIT, ierr);
125   arglis[0] = 1000;
126   minuit.mnexcm ("MIGRAD",arglis,1, ierr);
127 // std::cout << "Fitter return status = " << ierr << std::endl;
128   minuit.GetParameter(0, tzero, etzero );
129   return ierr;
130 }
131 
132 double CotNegLogT0Fitter::getLikelihood(double t0, CotT0Fitter::t0Input* input,
133                                         CotT0Fitter::t0Internal * internal) 
134 {
135   _input = input;
136   _internal = internal;
137   
138   double par[2];
139   par[0] = t0;
140   double llf;
141   if ( internal )
142   {
143     count_input( input, internal );
144     llf = pid_likelihood( par );
145     internal->chi2 = llf;
146   }
147   else
148   {
149     llf = pid_likelihood( par );
150   }
151   
152   return llf;
153 }
154 
155 double CotNegLogT0Fitter::pid_likelihood(double * par)
156 {
157   static const double norm = 2.5066283;
158   static const double minLog = log( std::numeric_limits<double>::min() );
159   double diff, sig, sigtail;
160   double outer, inner, prod;
161   
162 // Calculate the likelihood: 
163 //  
164 
165   outer = 0;
166 
167   for (int mi = 0; mi < _input->trk_t0.size(); mi++) { 
168       inner = 0;
169 
170       for (int m=PION; m<NPID; m++)     // skip electrons and muons...
171         {  
172                   prod = _input->trk_pidProb[mi][m];
173                   if (m==PION) {              
174                           diff  = _input->trk_t0[mi]  - par[0];
175                           sig   = _input->trk_et0[mi];
176                   } else if (m==KAON) {       
177                           diff  = _input->trk_kt0[mi]  - par[0];
178                           sig   = _input->trk_ekt0[mi];
179                   } else if (m==PROTON) {             
180                           diff  = _input->trk_pt0[mi]  - par[0];
181                           sig   = _input->trk_ept0[mi];
182                   }
183                           //          sigtail = 0.6;
184               if (sig!=0)
185                           prod *= exp(-diff*diff/(2*sig*sig))/(norm*sig)  
186 //                      * 0.9 + 0.1*exp(-diff*diff/(2*sigtail*sigtail))/(norm*sigtail)
187                                  ;
188                   inner += prod;
189         }
190 // std::cout << "  result inner = " 
191 //           << std::scientific << std::setprecision(20) << inner << std::endl;
192         if (inner>0.) {
193         double logInner = log( inner );
194         
195         // Don't allow log( inner ) values below that 
196         // defined by log( numeric_limits::min() ) in order
197         // to avoid discontinuities in likelihood function
198         //
199                 outer += ( logInner < minLog ? minLog : logInner );
200 
201 // std::cout << "  result log(inner) = "
202 //           << std::scientific << std::setprecision(20) 
203 //        << loginner << std::endl;
204 //         if ( loginner < minlog )
205 //      {
206 //        std::cout << "  underlog detected, setting log(inner) to "
207 //                  << minlog << std::endl;
208 //      }
209         } else {
210         // Add minimum possible log(double) value in order
211         // to maintain continuity of the likelihood function
212         //
213                 outer += minLog;
214 
215 //         double d = log( 1 / std::numeric_limits<double>::max() );
216 //         double d = log( std::numeric_limits<double>::min() );
217 //         std::cout << "zero inner. adding log(min) = log(" 
218 //                   << std::scientific << std::setprecision(20) 
219 //                << std::numeric_limits<double>::min() << ") = "
220 //                << d << std::endl;
221         }
222   }
223 // std::cout << "result final = " << -2*outer << std::endl;
224         return -2*outer;
225 }
226 
227 //global variables/functions for minuit:
228 CotNegLogT0Fitter* CotNegLogT0Fitter::instance = 0;
229 
230 void CotNegLogT0Fitter::pid_likelihood_fcn
231 (int &npar, double *gin, double &f, double *par, int iflag){
232   assert( instance );
233   f = instance->pid_likelihood(par);
234   if ( instance->_internal ) ++instance->_internal->ncall;
235 }
236 
237 
238 // Count tracks in t0Input and put results
239 // in t0Internal
240 // Is this necessary?
241 void CotNegLogT0Fitter::count_input( CotT0Fitter::t0Input    * input,
242                                      CotT0Fitter::t0Internal * internal )
243 {
244   int ntrk       = 0;
245  for ( int it = 0, itEnd = input->trk_t0.size() ; it < itEnd ; ++it )
246   {
247   //  if ( input->trk_used )
248     {
249       ++ntrk;
250     }
251   }
252   
253   internal->ntrk = ntrk;
254   internal->ndof = ntrk - 1;
255 }
[ source navigation ]
[ diff markup ]
[ identifier search ]
[ freetext search ]
[ file search ]
This page was automatically generated by the LXR engine.
The LXR team
Send problems or questions to cdfcode@fnal.gov