Main Page | Namespace List | Class Hierarchy | Alphabetical List | Class List | Directories | File List | Namespace Members | Class Members | File Members | Related Pages

trk::Circe Class Reference

#include <Circe.h>

List of all members.

Public Member Functions

 Circe ()
void SeedProngs (const Measurement *meas, unsigned int nprong, double *seed) const
double Fit (int nprog, const Measurement *m)
void MakeProngs (std::vector< recobase::Prong > &prong)
unsigned int NDim () const
 Return the number of fit parameters.
double DoEval (const double *p) const
ROOT::Math::IMultiGenFunction * Clone () const

Static Public Attributes

const unsigned int kNmax = 20
 Largest number of prongs that algorithm will find.

Private Attributes

const Measurement * fM
 Set of detector measurements.
unsigned int fN
 Number of prongs to fit.
double fVtx [3]
 Vertex location (xyz,cm).
double fTheta [kNmax]
 Polar angles of prongs (rad.).
double fPhi [kNmax]
 Azimuthal angles of prongs (rad.).
std::vector< int > fProngAssn
 Association between.

Detailed Description

Fit vertex and n-prongs to detector hits

This multi-prong fitter attempts to solve the same problem as the original Circe fitter for which its named (D.E. Fries, SLAC-99 UC-34, 1969), although the implementation is quite a bit different.

Definition at line 32 of file Circe.h.

Constructor & Destructor Documentation

trk::Circe::Circe (   )  [inline]
 

Definition at line 37 of file Circe.h.

00037 { }

Member Function Documentation

ROOT::Math::IMultiGenFunction * Circe::Clone (   )  const
 

Clone the function -- not completely sure if this is really needed. Return 0 so that errors are fatal.

Definition at line 133 of file Circe.cxx.

00133 { return 0; }

double Circe::DoEval (  const double *  par  )  const
 

Return the chi^2 comparison of the vertex and prong list (stored in "par") to the detector hits.

Parameters:
par - Vertex and direction parameters they are stored in order as par[0] = theta, par[1] = phi,... par[2*n+2]=vx, par[2*n+3]=vy, par[2*n+4]=vz
Returns:
- chi^2 measure of goodness of fit

Definition at line 30 of file Circe.cxx.

References geo::DsqrToLine(), trk::Measurement::fHits, fM, fN, trk::Measurement::fW, kNmax, and recobase::CellHit::View().

00031 {
00032   unsigned register int i;
00033   unsigned register int j;
00034   
00035   // Assert that limit on the number of prongs is observed
00036   if (fN>kNmax) {
00037     std::cerr << __FILE__ << ":" << __LINE__ 
00038               << " Can't fit more than " << kNmax << " prongs." 
00039               << std::endl; 
00040     abort();
00041   }
00042   
00043   // Unpack the input data. For convenience "track" the particles by
00044   // projecting them from the vertex 10 meters
00045   double vtx[3];
00046   double endpoint[kNmax][3];
00047   double theta;
00048   double phi;
00049   double ds = 10.0E2; // path length in cm
00050   std::vector<unsigned int> prongassn(fM->fHits.size());
00051 
00052   // Vertex location is last on the list
00053   vtx[0] = par[2*fN+0];
00054   vtx[1] = par[2*fN+1]; 
00055   vtx[2] = par[2*fN+2];
00056 
00057   // Direction angles are the first 2n variables
00058   for (i=0; i<fN; ++i) {
00059     theta = par[2*i+0];
00060     phi   = par[2*i+1];
00061     endpoint[i][0] = vtx[0]+ds*sin(theta)*cos(phi);
00062     endpoint[i][1] = vtx[1]+ds*sin(theta)*sin(phi);
00063     endpoint[i][2] = vtx[2]+ds*cos(theta);
00064   }
00065 
00066   double x0;          // Transverse hit        location (cm)
00067   double x1;          // Transverse vertex     location (cm)
00068   double x2;          // Transverse endpoint   location (cm)
00069   double y0;          // Longitudinal hit      location (cm)
00070   double y1;          // Longitudinal vertex   location (cm)
00071   double y2;          // Longitudinal endpoint location (cm)
00072   double dsqr;        // Distance^2 from prong to hit
00073   double dsqrmin;     // Shortest distance to any prong
00074   double dot;
00075   double score = 0.0; // Chi^2 score of comparison
00076   
00077   // Loop on all hits and find the prong they are closest to
00078   for (j=0; j<fM->fHits.size(); ++j) {
00079     const recobase::CellHit& h = (*fM->fHits[j]);
00080     
00081     x0 = h.ZCPos();
00082     y0 = h.TCPos();
00083 
00084     dsqrmin = 99.E99;
00085     for (i=0; i<fN; ++i) { 
00086       if (h.View()==geo::kX) {
00087         x1 = vtx[2];
00088         y1 = vtx[0];
00089         x2 = endpoint[i][2];
00090         y2 = endpoint[i][0];
00091       }
00092       else if (h.View()==geo::kY) {
00093         x1 = vtx[2];
00094         y1 = vtx[1];
00095         x2 = endpoint[i][2];
00096         y2 = endpoint[i][1];
00097       }
00098       else abort();
00099       // Check if the point is in the forward dirction. If not use the
00100       // distance to the vertex location as the score
00101       dot = (x0-x1)*(x2-x1) + (y0-y1)*(y2-y1);
00102       if (dot>=0.0) {
00103         dsqr = geo::DsqrToLine(x0, y0, x1, y1, x2, y2);
00104       }
00105       else {
00106         dsqr = pow(x0-x1,2)+pow(y0-y1,2);
00107       }
00108       if (dsqr<dsqrmin) dsqrmin = dsqr;
00109     }
00110     score += fM->fW[j]*dsqrmin;
00111   }
00112   /*
00113   std::cout << "chi^2=" << score << "\n";
00114   std::cout << "vtx=("
00115             << par[2*fN+0] << " " 
00116             << par[2*fN+1] 
00117             << par[2*fN+2] 
00118             << ")" << std::endl;
00119   for (unsigned int i=0; i<fN; ++i) {
00120     std::cout << "( " 
00121               << par[2*i+0] << ","
00122               << par[2*i+1] << ")" << std::endl;
00123   }
00124   */
00125   return score;
00126 }

double Circe::Fit (  int  nprong,
const Measurement *  meas
) 
 

Find the best fit for nprong prongs to the measurements in "meas"

Parameters:
nprong - Number of prongs to assume in the fit
meas - Collection of hits and weights to fit
Returns:
chi2 - Of best fit

Definition at line 195 of file Circe.cxx.

References trk::Measurement::fHits, fM, fN, fPhi, fTheta, fVtx, kNmax, and SeedProngs().

Referenced by trk::CirceFit::Reco().

00196 {
00197   double seed[2*kNmax+3];
00198 
00199   // Setup the fitter
00200   fN = nprong;
00201   fM = meas;
00202 
00203   // Seed the fit
00204   this->SeedProngs(meas, nprong, seed);
00205 
00206   // Setup one of the root minimizers
00207   // ROOT::Minuit2::Minuit2Minimizer mini(ROOT::Minuit2::kSimplex);
00208   ROOT::Minuit2::Minuit2Minimizer mini;
00209   mini.SetPrintLevel(0);
00210   mini.SetTolerance(0.1);
00211   mini.SetStrategy(2);
00212   mini.SetMaxFunctionCalls(10000);
00213 
00214   for (unsigned int i=0; i<fN; ++i) {
00215     char buff1[16];
00216     char buff2[16];
00217     sprintf(buff1, "theta_%d",i);
00218     sprintf(buff2, "phi_%d",  i);
00219     mini.SetVariable(2*i+0, buff1, seed[2*i+0], 0.1);
00220     mini.SetVariable(2*i+1, buff2, seed[2*i+1], 0.1);
00221   }
00222   mini.SetVariable(2*fN+0, "x0", seed[2*fN+0], 10.0);
00223   mini.SetVariable(2*fN+1, "y0", seed[2*fN+1], 10.0);
00224   mini.SetVariable(2*fN+2, "z0", seed[2*fN+2], 10.0);
00225 
00226   mini.SetFunction(*this);
00227   mini.Minimize();
00228   
00229   // Re-do the fit to guard against local minima
00230   for (unsigned int i=0; i<fN; ++i) {
00231     char buff1[16];
00232     char buff2[16];
00233     sprintf(buff1, "theta_%d",i);
00234     sprintf(buff2, "phi_%d",  i);
00235     mini.SetVariable(2*i+0, buff1, mini.X()[2*i+0], 0.1);
00236     mini.SetVariable(2*i+1, buff2, mini.X()[2*i+1], 0.1);
00237   }
00238   mini.SetVariable(2*fN+0, "x0", mini.X()[2*fN+0], 10.0);
00239   mini.SetVariable(2*fN+1, "y0", mini.X()[2*fN+1], 10.0);
00240   mini.SetVariable(2*fN+2, "z0", mini.X()[2*fN+2], 10.0);
00241   mini.Minimize();
00242 
00243   // Unpack the results of the fit
00244   for (unsigned int i=0; i<fN; ++i) {
00245     fTheta[i] = mini.X()[2*i+0];
00246     fPhi[i]   = mini.X()[2*i+1];
00247   }
00248   fVtx[0] = mini.X()[2*fN+0];
00249   fVtx[1] = mini.X()[2*fN+1];
00250   fVtx[2] = mini.X()[2*fN+2];
00251   
00252   // Return the standard deviation as a measure of goodness and
00253   // doneness
00254   double ndof = (float)((2*meas->fHits.size())-(2*fN+3));
00255   if (ndof<1.0) ndof = 1.0;
00256   return sqrt(mini.MinValue()/ndof);
00257 }

void Circe::MakeProngs (  std::vector< recobase::Prong > &  prong  ) 
 

Definition at line 261 of file Circe.cxx.

References fPhi, fTheta, fVtx, recobase::Prong::SetStart(), and recobase::Prong::SetStop().

Referenced by trk::CirceFit::Reco().

00262 {
00263   for (unsigned int i=0; i<fN; ++i) {
00264     recobase::Prong p;
00265 
00266     p.SetStart(fVtx);
00267     
00268     double ds = 1.E2; // track length in meters
00269     double xend[3];
00270     xend[0] = fVtx[0] + ds*sin(fTheta[i])*cos(fPhi[i]);
00271     xend[1] = fVtx[1] + ds*sin(fTheta[i])*sin(fPhi[i]);
00272     xend[2] = fVtx[2] + ds*cos(fTheta[i]);
00273     p.SetStop(xend);
00274 
00275     prong.push_back(p);
00276   }
00277 }

unsigned int Circe::NDim (   )  const
 

Return the number of fit parameters.

Definition at line 17 of file Circe.cxx.

References fN.

00017 { return 2*fN+3; }

void Circe::SeedProngs (  const Measurement *  meas,
unsigned int  nprong,
double *  seed
)  const
 

Pick a reasonable starting point for the fit

Parameters:
meas - the set of detector measurements
nprong - the number of prongs to be fit
seed - on return the seed parameters for the fit

Definition at line 143 of file Circe.cxx.

References trk::Measurement::fHits, fN, trk::Measurement::fW, and recobase::CellHit::View().

Referenced by Fit().

00146 {
00147   // Choose vertex position at most upstream end in z
00148   // direction. That's where the neutrinos come from...
00149   double sumx=0.0;
00150   double sumy=0.0;
00151   double sumz=0.0;
00152   double sumwx=0.0;
00153   double sumwy=0.0;
00154   double sumwz=0.0;
00155   for (unsigned int i=0; i<meas->fHits.size(); ++i) {
00156     const recobase::CellHit& h = (*meas->fHits[i]);
00157     if (h.View()==geo::kX) {
00158       sumx  += meas->fW[i]*h.TCPos();
00159       sumwx += meas->fW[i];
00160     }
00161     if (h.View()==geo::kY) {
00162       sumy  += meas->fW[i]*h.TCPos();
00163       sumwy += meas->fW[i];
00164     }
00165     sumz  += meas->fW[i]*h.ZCPos();
00166     sumwz += meas->fW[i];
00167   }
00168   seed[2*fN+0] = sumx/sumwx;
00169   seed[2*fN+1] = sumy/sumwy;
00170   seed[2*fN+2] = sumz/sumwz;
00171   
00172   // See the directions so that the first goes along the beam
00173   // direction and the others are distributed around in the x-y plane
00174   for (unsigned int i=0; i<nprong; ++i) {
00175     if (i==0) {
00176       seed[2*i+0] = 0.0;
00177       seed[2*i+1] = 0.0;
00178     }
00179     else {
00180       seed[2*i+0] = 0.0;
00181       seed[2*i+1] = M_PI*(float)i/(float)(nprong-1);
00182     }
00183   }
00184 }

Member Data Documentation

const Measurement* trk::Circe::fM [private]
 

Set of detector measurements.

Definition at line 51 of file Circe.h.

Referenced by DoEval(), and Fit().

unsigned int trk::Circe::fN [private]
 

Number of prongs to fit.

Definition at line 52 of file Circe.h.

Referenced by DoEval(), Fit(), NDim(), and SeedProngs().

double trk::Circe::fPhi[kNmax] [private]
 

Azimuthal angles of prongs (rad.).

Definition at line 55 of file Circe.h.

Referenced by Fit(), and MakeProngs().

std::vector<int> trk::Circe::fProngAssn [private]
 

Association between.

Definition at line 56 of file Circe.h.

double trk::Circe::fTheta[kNmax] [private]
 

Polar angles of prongs (rad.).

Definition at line 54 of file Circe.h.

Referenced by Fit(), and MakeProngs().

double trk::Circe::fVtx[3] [private]
 

Vertex location (xyz,cm).

Definition at line 53 of file Circe.h.

Referenced by Fit(), and MakeProngs().

const unsigned int trk::Circe::kNmax = 20 [static]
 

Largest number of prongs that algorithm will find.

Definition at line 35 of file Circe.h.

Referenced by DoEval(), and Fit().

The documentation for this class was generated from the following files:
Generated on Sun Nov 22 04:45:33 2009 for NOvA Offline by  doxygen 1.3.9.1