Circe.cxx

Go to the documentation of this file.

#include "TrackFit/Circe.h"
#include <iostream>
#include <cmath>
#include "Minuit2/Minuit2Minimizer.h"
#include "Geometry/geo.h"
#include "RecoBase/CellHit.h"
#include "RecoBase/Prong.h"
using namespace trk;

unsigned int Circe::NDim() const { return 2*fN+3; }

double Circe::DoEval(const double* par) const
{
  unsigned register int i;
  unsigned register int j;
  
  // Assert that limit on the number of prongs is observed
  if (fN>kNmax) {
    std::cerr << __FILE__ << ":" << __LINE__ 
              << " Can't fit more than " << kNmax << " prongs." 
              << std::endl; 
    abort();
  }
  
  // Unpack the input data. For convenience "track" the particles by
  // projecting them from the vertex 10 meters
  double vtx[3];
  double endpoint[kNmax][3];
  double theta;
  double phi;
  double ds = 10.0E2; // path length in cm
  std::vector<unsigned int> prongassn(fM->fHits.size());

  // Vertex location is last on the list
  vtx[0] = par[2*fN+0];
  vtx[1] = par[2*fN+1]; 
  vtx[2] = par[2*fN+2];

  // Direction angles are the first 2n variables
  for (i=0; i<fN; ++i) {
    theta = par[2*i+0];
    phi   = par[2*i+1];
    endpoint[i][0] = vtx[0]+ds*sin(theta)*cos(phi);
    endpoint[i][1] = vtx[1]+ds*sin(theta)*sin(phi);
    endpoint[i][2] = vtx[2]+ds*cos(theta);
  }

  double x0;          // Transverse hit        location (cm)
  double x1;          // Transverse vertex     location (cm)
  double x2;          // Transverse endpoint   location (cm)
  double y0;          // Longitudinal hit      location (cm)
  double y1;          // Longitudinal vertex   location (cm)
  double y2;          // Longitudinal endpoint location (cm)
  double dsqr;        // Distance^2 from prong to hit
  double dsqrmin;     // Shortest distance to any prong
  double dot;
  double score = 0.0; // Chi^2 score of comparison
  
  // Loop on all hits and find the prong they are closest to
  for (j=0; j<fM->fHits.size(); ++j) {
    const recobase::CellHit& h = (*fM->fHits[j]);
    
    x0 = h.ZCPos();
    y0 = h.TCPos();

    dsqrmin = 99.E99;
    for (i=0; i<fN; ++i) { 
      if (h.View()==geo::kX) {
        x1 = vtx[2];
        y1 = vtx[0];
        x2 = endpoint[i][2];
        y2 = endpoint[i][0];
      }
      else if (h.View()==geo::kY) {
        x1 = vtx[2];
        y1 = vtx[1];
        x2 = endpoint[i][2];
        y2 = endpoint[i][1];
      }
      else abort();
      // Check if the point is in the forward dirction. If not use the
      // distance to the vertex location as the score
      dot = (x0-x1)*(x2-x1) + (y0-y1)*(y2-y1);
      if (dot>=0.0) {
        dsqr = geo::DsqrToLine(x0, y0, x1, y1, x2, y2);
      }
      else {
        dsqr = pow(x0-x1,2)+pow(y0-y1,2);
      }
      if (dsqr<dsqrmin) dsqrmin = dsqr;
    }
    score += fM->fW[j]*dsqrmin;
  }
  /*
  std::cout << "chi^2=" << score << "\n";
  std::cout << "vtx=("
            << par[2*fN+0] << " " 
            << par[2*fN+1] 
            << par[2*fN+2] 
            << ")" << std::endl;
  for (unsigned int i=0; i<fN; ++i) {
    std::cout << "( " 
              << par[2*i+0] << ","
              << par[2*i+1] << ")" << std::endl;
  }
  */
  return score;
}

//......................................................................
ROOT::Math::IMultiGenFunction* Circe::Clone() const { return 0; }

//......................................................................
void Circe::SeedProngs(const Measurement* meas,
                       unsigned int       nprong,
                       double*            seed) const
{
  // Choose vertex position at most upstream end in z
  // direction. That's where the neutrinos come from...
  double sumx=0.0;
  double sumy=0.0;
  double sumz=0.0;
  double sumwx=0.0;
  double sumwy=0.0;
  double sumwz=0.0;
  for (unsigned int i=0; i<meas->fHits.size(); ++i) {
    const recobase::CellHit& h = (*meas->fHits[i]);
    if (h.View()==geo::kX) {
      sumx  += meas->fW[i]*h.TCPos();
      sumwx += meas->fW[i];
    }
    if (h.View()==geo::kY) {
      sumy  += meas->fW[i]*h.TCPos();
      sumwy += meas->fW[i];
    }
    sumz  += meas->fW[i]*h.ZCPos();
    sumwz += meas->fW[i];
  }
  seed[2*fN+0] = sumx/sumwx;
  seed[2*fN+1] = sumy/sumwy;
  seed[2*fN+2] = sumz/sumwz;
  
  // See the directions so that the first goes along the beam
  // direction and the others are distributed around in the x-y plane
  for (unsigned int i=0; i<nprong; ++i) {
    if (i==0) {
      seed[2*i+0] = 0.0;
      seed[2*i+1] = 0.0;
    }
    else {
      seed[2*i+0] = 0.0;
      seed[2*i+1] = M_PI*(float)i/(float)(nprong-1);
    }
  }
}

//......................................................................
double Circe::Fit(int nprong, const Measurement* meas) 
{
  double seed[2*kNmax+3];

  // Setup the fitter
  fN = nprong;
  fM = meas;

  // Seed the fit
  this->SeedProngs(meas, nprong, seed);

  // Setup one of the root minimizers
  // ROOT::Minuit2::Minuit2Minimizer mini(ROOT::Minuit2::kSimplex);
  ROOT::Minuit2::Minuit2Minimizer mini;
  mini.SetPrintLevel(0);
  mini.SetTolerance(0.1);
  mini.SetStrategy(2);
  mini.SetMaxFunctionCalls(10000);

  for (unsigned int i=0; i<fN; ++i) {
    char buff1[16];
    char buff2[16];
    sprintf(buff1, "theta_%d",i);
    sprintf(buff2, "phi_%d",  i);
    mini.SetVariable(2*i+0, buff1, seed[2*i+0], 0.1);
    mini.SetVariable(2*i+1, buff2, seed[2*i+1], 0.1);
  }
  mini.SetVariable(2*fN+0, "x0", seed[2*fN+0], 10.0);
  mini.SetVariable(2*fN+1, "y0", seed[2*fN+1], 10.0);
  mini.SetVariable(2*fN+2, "z0", seed[2*fN+2], 10.0);

  mini.SetFunction(*this);
  mini.Minimize();
  
  // Re-do the fit to guard against local minima
  for (unsigned int i=0; i<fN; ++i) {
    char buff1[16];
    char buff2[16];
    sprintf(buff1, "theta_%d",i);
    sprintf(buff2, "phi_%d",  i);
    mini.SetVariable(2*i+0, buff1, mini.X()[2*i+0], 0.1);
    mini.SetVariable(2*i+1, buff2, mini.X()[2*i+1], 0.1);
  }
  mini.SetVariable(2*fN+0, "x0", mini.X()[2*fN+0], 10.0);
  mini.SetVariable(2*fN+1, "y0", mini.X()[2*fN+1], 10.0);
  mini.SetVariable(2*fN+2, "z0", mini.X()[2*fN+2], 10.0);
  mini.Minimize();

  // Unpack the results of the fit
  for (unsigned int i=0; i<fN; ++i) {
    fTheta[i] = mini.X()[2*i+0];
    fPhi[i]   = mini.X()[2*i+1];
  }
  fVtx[0] = mini.X()[2*fN+0];
  fVtx[1] = mini.X()[2*fN+1];
  fVtx[2] = mini.X()[2*fN+2];
  
  // Return the standard deviation as a measure of goodness and
  // doneness
  double ndof = (float)((2*meas->fHits.size())-(2*fN+3));
  if (ndof<1.0) ndof = 1.0;
  return sqrt(mini.MinValue()/ndof);
}

//......................................................................

void Circe::MakeProngs(std::vector<recobase::Prong>& prong) 
{
  for (unsigned int i=0; i<fN; ++i) {
    recobase::Prong p;

    p.SetStart(fVtx);
    
    double ds = 1.E2; // track length in meters
    double xend[3];
    xend[0] = fVtx[0] + ds*sin(fTheta[i])*cos(fPhi[i]);
    xend[1] = fVtx[1] + ds*sin(fTheta[i])*sin(fPhi[i]);
    xend[2] = fVtx[2] + ds*cos(fTheta[i]);
    p.SetStop(xend);

    prong.push_back(p);
  }
}

---