CirceFit.cxx

Go to the documentation of this file.

#include "TrackFit/CirceFit.h"
#include <vector>
#include "TH1F.h"
#include "RecoBase/CellHit.h"
#include "RecoBase/Prong.h"
#include "SimulationBase/MCTruth.h"
#include "TrackFit/Circe.h"
using namespace trk;
MODULE_DECL(trk::CirceFit);

CirceFit::CirceFit(const char* version) : 
  jobc::Module("trk::CirceFit"),
  fFinalVtxX(0),
  fFinalVtxY(0),
  fFinalVtxZ(0)
{
  this->SetCfgVersion(version); 
}

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

void CirceFit::Update(const cfg::Config& c)
{
  c("Debug").       Get(fDebug);
  c("InputCalHits").Get(fInputCalHits);
  c("WeightMethod").Get(fWeightMethod);
  c("MinXhits").    Get(fMinXhits);
  c("MinYhits").    Get(fMinYhits);
  c("Nmax").        Get(fNmax);
  c("Stop").        Get(fStop);
  c("NNZrange").    Get(fNNZrange);
  c("NNZedge").     Get(fNNZedge);
  c("NNTrange").    Get(fNNTrange);
  c("NNTedge").     Get(fNNTedge);
  c("OtherViewW0"). Get(fOtherViewW0);

}

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

jobc::Result CirceFit::Ana (const edm::EventHandle& evt) 
{
  // The Monte Carlo true 4-vectors for the event
  std::vector<const simb::MCTruth*> mctruth;
  try {
    evt.MC().Get("./",mctruth);
  }
  catch (...) {
    return jobc::kFailed;
  }
  
  // We're only set up to handle one MCTruth object at a time...
  if (mctruth.size()==0) {
    return jobc::kFailed;
  }
  else if (mctruth.size()>1) {
    std::cout << __FILE__ << ":" << __LINE__ 
              << " Not ready for overlap events yet..." << std::endl;
    return jobc::kFailed;
  }
  
  const simb::MCNeutrino* nu = 0;
  nu = mctruth[0]->GetNeutrino();
  if (nu==0) {
    std::cout << __FILE__ << ":" << __LINE__ 
              << " No neutrino??" << std::endl;
    return jobc::kFailed;
  }
  
  // Plot where the true vertex is compared to where we are seeding
  double x0mc = nu->Lepton().Vx();
  double y0mc = nu->Lepton().Vy();
  double z0mc = nu->Lepton().Vz();
  
  // Get seed from the circe fitter...
  // double x0seed, y0seed, z0seed;
  // Circe fitter;
  // fitter.Seed(....)

  // Compare seed and true locations.
  // fVtxXHisto->Fill(x0seed-x0mc);
  // fVtxYHisto->Fill(y0seed-y0mc);
  // fVtxZHisto->Fill(z0seed-z0mc);

  // Pull the prongs we fit out of the event
  std::vector<const recobase::Prong*> prong;
  try {
    evt.Reco().Get("./circe",prong);
  }
  catch (...) { }
  
  if (prong.size()>1 && nu->CCNC()==simb::kCC) {
    if (fFinalVtxX==0) {
      fFinalVtxX = new TH1F("fFinalVtxX",
                            "Circe Vertex X; Xfit-Xmc (cm); Events",
                            200, -50.0,50.0);
      fFinalVtxY = new TH1F("fFinalVtxY",
                            "Circe Vertex Y; Yfit-Ymc (cm); Events",
                            200, -50.0,50.0);
      fFinalVtxZ = new TH1F("fFinalVtxZ",
                            "Circe Vertex Z; Zfit-Zmc (cm); Events",
                            200, -50.0,50.0);
    }
    fFinalVtxX->Fill(prong[0]->StartPos()[0]-x0mc);
    fFinalVtxY->Fill(prong[0]->StartPos()[1]-y0mc);
    fFinalVtxZ->Fill(prong[0]->StartPos()[2]-z0mc);
  }
  

  return jobc::kPassed;
}

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

jobc::Result CirceFit::Reco(edm::EventHandle& evt) 
{
  // Pull the calibrated hits out of the event
  std::vector<const recobase::CellHit*> hits;
  try {
    evt.Cal().Get(fInputCalHits.c_str(), hits);
  }
  catch (...) {
    // Failed to find any CalHits. Fail the event and keep going...
    if (fDebug>0) {
      std::cout << "Failed to find any CalHits" << std::endl;
    }
    return jobc::kFailed;
  }
  if (fDebug>0) {
    std::cout << __FILE__ << ":" << __LINE__ 
              << " Found " << hits.size() << " hits in event." 
              << std::endl;
  }

  // Select hits to try to reconstruct and assign them weights
  unsigned int nxview = 0;
  unsigned int nyview = 0;
  trk::Measurement m;
  for (unsigned int i=0; i<hits.size(); ++i) {
    // Assign weight to hit.
    double w = 0.0;
    bool   ok;
    float  pe;
    double svw = fOtherViewW0;
    switch (fWeightMethod) {
    case 1:
      ok = hits[i]->PEAC(pe);
      if (ok) w = pe;
      else    w = 1.0;
      break;
    case 2:
      // Weight based on nearest neighbor count
      w = 1.0;
      for (unsigned int j=0; j<hits.size(); ++j) {
        if (i==j) continue;
        double dT = fabs(hits[i]->TCPos()-hits[j]->TCPos());
        double dZ = fabs(hits[i]->ZCPos()-hits[j]->ZCPos());
        double wT = 1.0 - 1.0/(1.0+exp((fNNTrange-dT)/fNNTedge));
        double wZ = 1.0 - 1.0/(1.0+exp((fNNZrange-dZ)/fNNZedge));
        if (hits[i]->View() == hits[j]->View()) w += wT+wZ;
        else                                    w += svw*wZ;
      }
      break;
    default:
      w = 1.0;
      break;
    }
    if (fDebug>0) std::cout << "Weight=" << w << std::endl;
    m.fHits.push_back(hits[i]);
    m.fW.push_back(w);
    
    if (hits[i]->View()==geo::kX) ++nxview;
    if (hits[i]->View()==geo::kY) ++nyview;
  }

  // Only bother if we have some minimum number of hits in each view
  if (nxview<fMinXhits || nyview<fMinYhits) return jobc::kFailed;
  
  // Build the fitter and pass it the data to fit
  double stdev;
  Circe fitter;
  for (int i=1; i<=fNmax; ++i) {
    stdev = fitter.Fit(i,&m);
    if (fDebug>0) {
      std::cout << __FILE__ << ":" << __LINE__ << " " 
                << i << " stdev=" << stdev << std::endl;
    }
    if (stdev<fStop) break;
  }

  std::vector<recobase::Prong> prong;
  fitter.MakeProngs(prong);
  
  // Store results into the event
  if (evt.Reco().GetFolder("./circe")==0) {
    evt.Reco().MakeFolder("./circe");
  }
  for (unsigned int i=0; i<prong.size(); ++i) {
    evt.Reco().Put(prong[i],"./circe");
  }
  return jobc::kPassed;
}

---