DimuonFitter_module.cc

Go to the documentation of this file.

///
/// \file    DimuonFitter_module.cc
/// \brief   Run the break-point fitter
/// \author  messier@indiana.edu
/// \version $Id:$
///
#include <algorithm>
//
#include "TH2F.h"
#include "TNtuple.h"
//
#include "art/Framework/Core/ModuleMacros.h"
#include "art/Framework/Core/EDProducer.h"
#include "art/Framework/Services/Optional/TFileService.h"
#include "canvas/Persistency/Common/FindManyP.h"
//
#include "Geometry/Geometry.h"
#include "GeometryObjects/CellGeo.h"
#include "RecoBase/Vertex.h"
#include "RecoBase/Cluster.h"
#include "RecoBase/Prong.h"
#include "RecoBase/Track.h"
#include "RecoBase/FitSum.h"
//
#include "Utilities/AssociationUtil.h"
//
#include "BreakPointFitter/func/HitList3D.h"
#include "BreakPointFitter/func/TrackBasis.h"
#include "BreakPointFitter/func/Path.h"
#include "BreakPointFitter/func/ScatteringSurfaces.h"
#include "BreakPointFitter/func/Lutz.h"
#include "BreakPointFitter/func/BPException.h"

namespace bpfit { class DimuonFitter; }

class bpfit::DimuonFitter : public art::EDProducer
{
public:
  const art::ProductToken<std::vector<rb::Cluster>> fSliceToken;

public:
  // Member data
  int          fRun;       ///< Slice identifiers (run, subrun, event, slice)
  int          fSubrun;    ///< are useful to have as class scope for logging
  int          fEvent;     ///< debugging information in each of the methods
  int          fSlice;     ///< used during the reconstruction.
  int          fVerbosity; ///< 0-100, how much text output to generate
  TNtuple*     fPass1Nt;   ///< pass one ntuple 
public:
  //......................................................................
  
  void reconfigure(fhicl::ParameterSet const& p)
  {
    fVerbosity = p.get<int>("Verbosity");
  }
  
  //......................................................................

  explicit DimuonFitter(fhicl::ParameterSet const& p) :
    fSliceToken(consumes<std::vector<rb::Cluster>>(p.get<std::string>("SlicerLabel"))),
    fRun(0),
    fSubrun(0),
    fEvent(0),
    fSlice(0)
  {
    this->produces< std::vector<rb::Vertex> >();
    this->produces< std::vector<rb::Track> >();
    this->reconfigure(p);
  }

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

  virtual ~DimuonFitter() { }

  void beginJob(){
  art::ServiceHandle<art::TFileService> f;
  fPass1Nt = f->make<TNtuple>("p1", 
                              "pass 1 ntuple",
                              "run:event:slice:"
                              "np:nhitx:nhity:"
                              "np5:"
                              "nc5l:nc5r:"
                              "nc5t:nc5b");
  
  }
  //......................................................................

  void FindSlices(const art::Event& evt, 
                  art::Handle< std::vector<rb::Cluster> >& h,
                  art::PtrVector<rb::Cluster>& s) 
  {
    unsigned int i;
    evt.getByToken(fSliceToken, h);
    for (i=0; i<h->size(); ++i) {
      s.push_back(art::Ptr<rb::Cluster>(h, i));
    }
    if (fVerbosity>10) {
      std::cout << __FILE__ << ":" << __LINE__ 
                << " Found " << s.size() << " slices."
                << std::endl;
    }
  }

  //......................................................................
  //
  // Make an analysis of the variables corresponding to a single slice.
  //
  // Return "true" if this is a slice we want to keep in the analysis
  // "false" if we thinkg this slice is ~100% likely to be background.
  //
  // This selection is reduction pass 1 for the analysis
  //
  //
  bool AnaSlice(const rb::Cluster* s) {
    unsigned int i;
    //
    // Define the cell boundaries at the edges of the detector
    //
    int icl  =  5; // 5 cells on the left side
    int icr  = 90; // 5 cells on the right side
    int icb  =  5; // 5 cells on the bottom face
    int ict  = 90; // 5 cells on the top face
    int ipf  =  5; // 5 planes on the front face
    int nc5l =  0; // number hits in the boundary on the left
    int nc5r =  0; // "" on the right
    int nc5b =  0; // "" on the bottom
    int nc5t =  0; // "" on the top
    int npl  =  0; // Total number of planes in cluster
    int np5  =  0; // Number of hits in the first 5 planes
    std::vector<int> hitp(1000);
    for (i=0; i<s->NCell(); ++i) {
      int p = s->Cell(i)->Plane();
      int c = s->Cell(i)->Cell();
      //
      // Count the number of hit planes
      //
      if (hitp[s->Cell(i)->Plane()]==0) ++npl;
      hitp[s->Cell(i)->Plane()] = 1;
      //
      // Count cell hits near walls. Don't double count hits in the
      // corners (hence the "else")
      //
      if (p<ipf) { 
        ++np5;
      }
      else {
        if (s->Cell(i)->View()==geo::kX) {
          if (c<icl) ++nc5r;
          if (c>icr) ++nc5l;
        }
        else if (s->Cell(i)->View()==geo::kY) {
          if (c<icb) ++nc5b;
          if (c>ict) ++nc5t;
        }
      }
    }
    float x[16];
    x[0]  = fRun;
    x[1]  = fEvent;
    x[2]  = fSlice;
    x[3]  = npl;
    x[4]  = s->NXCell();
    x[5]  = s->NYCell();
    x[6]  = np5;
    x[7]  = nc5l;
    x[8]  = nc5r;
    x[9]  = nc5b;
    x[10] = nc5t;
    fPass1Nt->Fill(x);

    return true;
  }
  //......................................................................

  double FindVertexZ(const rb::Cluster* c) {
    if (fVerbosity>20) {
      std::cout << __FILE__ << ":" << __LINE__
                << " Start 'FindVertexZ'. "
                << " Cluster has " << c->NCell() << " cells" << std::endl;
    }
    unsigned int i,  j;
    unsigned int ip, ic;
    //
    // The way we are going to find a candidate vertex z location is
    // to identify the first plane which is followed by a reasonably
    // continuous set of activity. So, we will loop over all the hits,
    // record their plane locations. Then we will look for the first
    // plane in the list where if I step n positions forward in the
    // list the plane number doesn't change more than m. While we are
    // making this list, keep track of the lowest z position we've
    // seen. In cases where we don't have enough hits to apply the
    // algorithm above, we'll just use that lowest z value.
    //
    std::vector<unsigned int> p(c->NCell());
    unsigned int              pvtx = 999999;
    for (i=0; i<c->NCell(); ++i) {
      ip = c->Cell(i)->Plane();
      if (ip<pvtx) pvtx = ip;
      p[i] = ip;
    }
    if (fVerbosity>20) {
      std::cout << __FILE__ << ":" << __LINE__ 
                << " Set default vertex plane to pvtx=" 
                << pvtx << std::endl;
    }
    std::sort(p.begin(), p.end());

    const unsigned int N = 2;   // How far to skip forward in list
    const unsigned int M = N+2; // Plane numbers should be this close or closer
    for (i=0, j=N; j+1<c->NCell(); ++i,++j) {
      if ( (p[j]-p[i])<=M ) {
        pvtx = p[i];
        if (fVerbosity>20) {
          std::cout << __FILE__ ":" << __LINE__ 
                    << " Plane " << pvtx << " satisfies "
                    << " N= " << N
                    << " M= " << M
                    << " pi = " << p[i] << " condition."
                    << " pj = " << p[j] << std::endl;
        }
        break;
      }
    }
    //
    // Now we have to match the plane number to a z location.
    //
    art::ServiceHandle<geo::Geometry> geom;
    for (i=0; i<c->NCell(); ++i) {
      ip = c->Cell(i)->Plane();
      if (ip==pvtx) {
        ic = c->Cell(i)->Cell();
        double xyz[3];
        geom->Plane(ip)->Cell(ic)->GetCenter(xyz);
        if (fVerbosity>10) {
          std::cout << __FILE__ << ":" << __LINE__ 
                    << " Found vertex z to z=" 
                    << xyz[2] << " cm." << std::endl;
        }
        return xyz[2];
      }
    }
    //
    // Shouldn't ever get here.
    //
    abort();
    return 0.;
  }
  
  //......................................................................

  void FitView(const rb::Cluster& clust, geo::View_t view, double zvtx)
  {
    if (fVerbosity>20) {
      std::cout << __FILE__ << ":" << __LINE__ 
                << " Starting FitView for view " << view
                << std::endl;
    }
    unsigned int i;
    const double zero_weight = 1e-9;
    //
    // My initial weighting scheme depends on the location of the hits
    // inside the slice so first I need to work out the bounding box
    // that contains the cluster
    //
    unsigned int plo = 999999;
    unsigned int phi = 0;
    unsigned int clo = 999999;
    unsigned int chi = 0;
    for (i=0; i<clust.NCell(); ++i) {
      const rb::CellHit* h = clust.Cell(i).get();
      unsigned int ip = h->Plane();
      unsigned int ic = h->Cell();
      if (ip<plo) plo = ip;
      if (ip>phi) phi = ip;
      if (ic<clo) clo = ic;
      if (ic>chi) chi = ic;
    }
    if (fVerbosity>20) {
      std::cout << __FILE__ << ":" << __LINE__ 
                << " Found slice extent " 
                << " [" << plo << "," << phi << "] "
                << " [" << clo << "," << chi << "]"
                << std::endl;
    }

    //
    // Initialize the data for the fit. Need x,z locations, error bars
    // on x and weights for each of the hitss
    //
    std::vector<double> x1; // Transverse hit locations (cm)
    std::vector<double> z1; // Longitudinal hit locations (cm)
    std::vector<double> s1; // Transverse uncertainty hit locations (cm)
    std::vector<double> w1; // Weight of the hit to be associated with track 1
    
    std::vector<double> x2; // Transverse hit locations (cm)
    std::vector<double> z2; // Longitudinal hit locations (cm)
    std::vector<double> s2; // Transverse uncertainty hit locations (cm)
    std::vector<double> w2; // Weight of the hit to be associated with track 1

    double zlo  =  999999;
    double zhi  = -999999;
    double rt12 = sqrt(12.0);
    for (i=0; i<clust.NCell(); ++i) {
      const rb::CellHit* h = clust.Cell(i).get();
      unsigned int ip = h->Plane();
      unsigned int ic = h->Cell();

      art::ServiceHandle<geo::Geometry> geom;
      const geo::PlaneGeo* p = geom->Plane(ip);
      const geo::CellGeo*  c = p->Cell(ic);
      geo::View_t          v = p->View();
      //
      // If the hit doesn't match the view we've been asked to fit,
      // skip
      //
      if (v != view) continue;
      
      double xyz[3];
      c->GetCenter(xyz, clust.W(h));
      if (xyz[2]<zlo) zlo = xyz[2];
      if (xyz[2]>zhi) zhi = xyz[2];
      double sigxy = c->HalfW()/rt12;

      
      if      (v==geo::kX) { x1.push_back(xyz[0]); x2.push_back(xyz[0]); }
      else if (v==geo::kY) { x1.push_back(xyz[1]); x2.push_back(xyz[1]); }
      else                 { abort(); }
      z1.push_back(xyz[2]);
      z2.push_back(xyz[2]);
      s1.push_back(sigxy);
      s2.push_back(sigxy);
      //
      // We have to break the symmetry in the weights somehow or the
      // two tracks will converge to the same solution. I'm splitting
      // according to cell number to that the "top" track generally
      // will be track one and the "bottom" track generally will be
      // track 2
      //
      double wp  = (float)(phi-ip)/(float)(phi-plo);
      double wc1 = (float)(chi-ic)/(float)(chi-clo);
      double wc2 = (float)(ic-clo)/(float)(chi-clo);
      wc1 = (wp*wp)*(wc1*wc1);
      wc2 = (wp*wp)*(wc2*wc2);
      wc1 = wc1/((wc1+wc2)+zero_weight);
      wc2 = 1.0-wc1;
      wc1 += zero_weight;
      wc2 += zero_weight;
      w1.push_back(wc1);
      w2.push_back(wc2);
    }
    if (fVerbosity>30) { 
      for (i=0; i<z1.size(); ++i) {
        std::cout << __FILE__ << ":" << __LINE__ << "\t"
                  << i << "\t"
                  << z1[i] << "\t" << z2[i] << "\t" 
                  << x1[i] << "\t" << x2[i] << "\t"
                  << s1[i] << "\t" << s2[i] << "\t" 
                  << w1[i] << "\t" << w2[i] << std::endl;
      }
    }
    if (fVerbosity>20) {
      std::cout << __FILE__ << ":" << __LINE__ 
                << " Initialized fits." 
                << std::endl;
    }

    //
    // Construct the scattering surfaces. Do something simple to get
    // started: A straight line down the middle of the
    // detector. Should be good enough to estimate the total radiation
    // lengths the tracks encounter.
    //
    int    pdg_muon = 13;
    double ds       = 50.0;
    double dx       = 1.0;
    double dt       = 0.3;
    double p1[3]    = {3,3,zlo};
    double p2[3]    = {3,3,zhi};
    std::cout << __FILE__ << ":" << __LINE__ << std::endl;
    std::vector<double>               s(z1);
    art::ServiceHandle<geo::Geometry> geom;
    bpfit::ScatteringSurfaces         scatter(pdg_muon, ds, dx, dt);
    sort(s.begin(), s.end());
    scatter.MakeSurfaces(*geom, p1, p2, s);
    if (fVerbosity>30) {
      std::cout << __FILE__ << ":" << __LINE__ << " Made "
                << scatter.N() << " scattering surfaces between z=" 
                << zlo << ", " << zhi << std::endl;
    }
    //
    // Now perform the fits
    //
    double lambda = 1/25.;
    for (unsigned int itr=0; itr<5; ++itr) {
      if (fVerbosity>30) {
        std::cout << __FILE__ << ":" << __LINE__ 
                  << " Fitting iteration " << itr 
                  << std::endl;
      }
      bpfit::Lutz lutz1(z1.size(), scatter.N());
      bpfit::Lutz lutz2(z2.size(), scatter.N());
      for (i=0; i<z1.size(); ++i) {
        lutz1.SetMeasurement(i, z1[i], x1[i], s1[i]/w1[i]);
      }
      for (i=0; i<z2.size(); ++i) {
        lutz2.SetMeasurement(i, z2[i], x2[i], s2[i]/w2[i]);
      }
      for (i=0; i<scatter.N(); ++i) {
        lutz1.SetScatteringPlane(i, scatter.S(i), scatter.SigAlpha(i));
        lutz2.SetScatteringPlane(i, scatter.S(i), scatter.SigAlpha(i));
      }

      double a1, b1, chi1;
      double a2, b2, chi2;
      lutz1.Fit(&a1, &b1, 0, &chi1);
      lutz2.Fit(&a2, &b2, 0, &chi2);
      if (fVerbosity>20) {
        std::cout << __FILE__ << ":" << __LINE__ 
                  << " Finished fits in view " << view << "."
                  << " chi1 = " << chi1 
                  << " chi2 = " << chi2
                  << std::endl;
      }
      //
      // Update the weights for the next iteration
      //
      lambda = (itr+1)*(itr+1)/25.0;
      for (i=0; i<z1.size(); ++i) {
        double ex1 = exp(-lutz1.Chi2XIi(i));
        double ex2 = exp(-lutz2.Chi2XIi(i));
        double ex3 = exp(-25.0*lambda);
        w1[i] = ex1/(ex3+ex1+ex2);
        w2[i] = ex2/(ex3+ex1+ex2);
        if (w1[i]<zero_weight) w1[i] = zero_weight;
        if (w2[i]<zero_weight) w2[i] = zero_weight;
      }
    }
  }

  //......................................................................
  
  void produce(art::Event& evt) 
  {
    //
    // Allocations for the track products and their associations
    //
    std::unique_ptr< std::vector<rb::Vertex> > 
      vtx(new std::vector<rb::Vertex>);

    fRun    = evt.run();
    fSubrun = evt.subRun();
    fEvent  = evt.id().event();
    fSlice  = 0;
    //
    // Extract the slices, vertices, and tracks we need as inputs
    //
    art::Handle< std::vector<rb::Cluster> > sliceh;
    art::PtrVector<rb::Cluster>             slicep;
    this->FindSlices(evt, sliceh, slicep);

    unsigned int islice;
    for (islice=0; islice<slicep.size(); ++islice) {
      //
      // Get the index of this slices and a pointer to it for
      // convenient access. Skip the noice slice, no physics there for
      // us.
      //
      fSlice                            = islice;
      const art::Ptr<rb::Cluster> slice = slicep[islice];
      if (slice->IsNoise()) continue;
      //
      // Analyze the slice and see if it passes selections for analysis
      //
      this->AnaSlice(slice.get());

      //
      // The track fits require a candidate z vertex location
      //
      double zvtx = this->FindVertexZ(slice.get());
      vtx->push_back(rb::Vertex(0,0,zvtx,0));

      // Then we want to fit the x hits to two tracks and the y hits to
      // two tracks. Likely this means creating a Lutz object, feeding in
      // the hits from the x or y view. Will need to settle on some
      // locations for scattering planes. Probably the hit associations
      // can be done with some sort of annealing...
      //
      this->FitView(*slice.get(), geo::kX, zvtx);
      this->FitView(*slice.get(), geo::kY, zvtx);
      
      // Now we'll need to match the four tracks across the two views.
      // this->ViewMatch(...);
      
      //
      // Annouce that we ran to the log file...
      //
      std::cout << __FILE__ << ":" << __LINE__ 
                << " Slice #: " << islice << std::endl;
    } // Loop on slices
    
    //
    // To keep ART quiet for now. Eventually should produce real track
    // fits
    //
    std::unique_ptr<std::vector<rb::Track> >
      tracks(new std::vector<rb::Track>);
    evt.put(std::move(vtx));
    evt.put(std::move(tracks));
  } // produce() method
};

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

DEFINE_ART_MODULE(bpfit::DimuonFitter)

////////////////////////////////////////////////////////////////////////