NOvA Offline: trk::CirceFit Class Reference

00017                                       : 
00018   jobc::Module("trk::CirceFit"),
00019   fFinalVtxX(0),
00020   fFinalVtxY(0),
00021   fFinalVtxZ(0)
00022 {
00023   this->SetCfgVersion(version); 
00024 }

00048 {
00049   // The Monte Carlo true 4-vectors for the event
00050   std::vector<const simb::MCTruth*> mctruth;
00051   try {
00052     evt.MC().Get("./",mctruth);
00053   }
00054   catch (...) {
00055     return jobc::kFailed;
00056   }
00057   
00058   // We're only set up to handle one MCTruth object at a time...
00059   if (mctruth.size()==0) {
00060     return jobc::kFailed;
00061   }
00062   else if (mctruth.size()>1) {
00063     std::cout << __FILE__ << ":" << __LINE__ 
00064               << " Not ready for overlap events yet..." << std::endl;
00065     return jobc::kFailed;
00066   }
00067   
00068   const simb::MCNeutrino* nu = 0;
00069   nu = mctruth[0]->GetNeutrino();
00070   if (nu==0) {
00071     std::cout << __FILE__ << ":" << __LINE__ 
00072               << " No neutrino??" << std::endl;
00073     return jobc::kFailed;
00074   }
00075   
00076   // Plot where the true vertex is compared to where we are seeding
00077   double x0mc = nu->Lepton().Vx();
00078   double y0mc = nu->Lepton().Vy();
00079   double z0mc = nu->Lepton().Vz();
00080   
00081   // Get seed from the circe fitter...
00082   // double x0seed, y0seed, z0seed;
00083   // Circe fitter;
00084   // fitter.Seed(....)
00085 
00086   // Compare seed and true locations.
00087   // fVtxXHisto->Fill(x0seed-x0mc);
00088   // fVtxYHisto->Fill(y0seed-y0mc);
00089   // fVtxZHisto->Fill(z0seed-z0mc);
00090 
00091   // Pull the prongs we fit out of the event
00092   std::vector<const recobase::Prong*> prong;
00093   try {
00094     evt.Reco().Get("./circe",prong);
00095   }
00096   catch (...) { }
00097   
00098   if (prong.size()>1 && nu->CCNC()==simb::kCC) {
00099     if (fFinalVtxX==0) {
00100       fFinalVtxX = new TH1F("fFinalVtxX",
00101                             "Circe Vertex X; Xfit-Xmc (cm); Events",
00102                             200, -50.0,50.0);
00103       fFinalVtxY = new TH1F("fFinalVtxY",
00104                             "Circe Vertex Y; Yfit-Ymc (cm); Events",
00105                             200, -50.0,50.0);
00106       fFinalVtxZ = new TH1F("fFinalVtxZ",
00107                             "Circe Vertex Z; Zfit-Zmc (cm); Events",
00108                             200, -50.0,50.0);
00109     }
00110     fFinalVtxX->Fill(prong[0]->StartPos()[0]-x0mc);
00111     fFinalVtxY->Fill(prong[0]->StartPos()[1]-y0mc);
00112     fFinalVtxZ->Fill(prong[0]->StartPos()[2]-z0mc);
00113   }
00114   
00115 
00116   return jobc::kPassed;
00117 }

00122 {
00123   // Pull the calibrated hits out of the event
00124   std::vector<const recobase::CellHit*> hits;
00125   try {
00126     evt.Cal().Get(fInputCalHits.c_str(), hits);
00127   }
00128   catch (...) {
00129     // Failed to find any CalHits. Fail the event and keep going...
00130     if (fDebug>0) {
00131       std::cout << "Failed to find any CalHits" << std::endl;
00132     }
00133     return jobc::kFailed;
00134   }
00135   if (fDebug>0) {
00136     std::cout << __FILE__ << ":" << __LINE__ 
00137               << " Found " << hits.size() << " hits in event." 
00138               << std::endl;
00139   }
00140 
00141   // Select hits to try to reconstruct and assign them weights
00142   unsigned int nxview = 0;
00143   unsigned int nyview = 0;
00144   trk::Measurement m;
00145   for (unsigned int i=0; i<hits.size(); ++i) {
00146     // Assign weight to hit.
00147     double w = 0.0;
00148     bool   ok;
00149     float  pe;
00150     double svw = fOtherViewW0;
00151     switch (fWeightMethod) {
00152     case 1:
00153       ok = hits[i]->PEAC(pe);
00154       if (ok) w = pe;
00155       else    w = 1.0;
00156       break;
00157     case 2:
00158       // Weight based on nearest neighbor count
00159       w = 1.0;
00160       for (unsigned int j=0; j<hits.size(); ++j) {
00161         if (i==j) continue;
00162         double dT = fabs(hits[i]->TCPos()-hits[j]->TCPos());
00163         double dZ = fabs(hits[i]->ZCPos()-hits[j]->ZCPos());
00164         double wT = 1.0 - 1.0/(1.0+exp((fNNTrange-dT)/fNNTedge));
00165         double wZ = 1.0 - 1.0/(1.0+exp((fNNZrange-dZ)/fNNZedge));
00166         if (hits[i]->View() == hits[j]->View()) w += wT+wZ;
00167         else                                    w += svw*wZ;
00168       }
00169       break;
00170     default:
00171       w = 1.0;
00172       break;
00173     }
00174     if (fDebug>0) std::cout << "Weight=" << w << std::endl;
00175     m.fHits.push_back(hits[i]);
00176     m.fW.push_back(w);
00177     
00178     if (hits[i]->View()==geo::kX) ++nxview;
00179     if (hits[i]->View()==geo::kY) ++nyview;
00180   }
00181 
00182   // Only bother if we have some minimum number of hits in each view
00183   if (nxview<fMinXhits || nyview<fMinYhits) return jobc::kFailed;
00184   
00185   // Build the fitter and pass it the data to fit
00186   double stdev;
00187   Circe fitter;
00188   for (int i=1; i<=fNmax; ++i) {
00189     stdev = fitter.Fit(i,&m);
00190     if (fDebug>0) {
00191       std::cout << __FILE__ << ":" << __LINE__ << " " 
00192                 << i << " stdev=" << stdev << std::endl;
00193     }
00194     if (stdev<fStop) break;
00195   }
00196 
00197   std::vector<recobase::Prong> prong;
00198   fitter.MakeProngs(prong);
00199   
00200   // Store results into the event
00201   if (evt.Reco().GetFolder("./circe")==0) {
00202     evt.Reco().MakeFolder("./circe");
00203   }
00204   for (unsigned int i=0; i<prong.size(); ++i) {
00205     evt.Reco().Put(prong[i],"./circe");
00206   }
00207   return jobc::kPassed;
00208 }

00029 {
00030   c("Debug").       Get(fDebug);
00031   c("InputCalHits").Get(fInputCalHits);
00032   c("WeightMethod").Get(fWeightMethod);
00033   c("MinXhits").    Get(fMinXhits);
00034   c("MinYhits").    Get(fMinYhits);
00035   c("Nmax").        Get(fNmax);
00036   c("Stop").        Get(fStop);
00037   c("NNZrange").    Get(fNNZrange);
00038   c("NNZedge").     Get(fNNZedge);
00039   c("NNTrange").    Get(fNNTrange);
00040   c("NNTedge").     Get(fNNTedge);
00041   c("OtherViewW0"). Get(fOtherViewW0);
00042 
00043 }


Public Member Functions
	CirceFit (const char *version)
void	Update (const cfg::Config &c)
jobc::Result	Reco (edm::EventHandle &evt)
jobc::Result	Ana (const edm::EventHandle &evt)
Private Attributes
int	fDebug
	Controls debug output 0=none, 99=scream.
std::string	fInputCalHits
	Where to find CalHits to reconstruct.
unsigned int	fMinXhits
	Required # of hits in x view.
unsigned int	fMinYhits
	Required # of hits in x view.
int	fNmax
	Hard stop after this many prongs.
double	fStop
	Stopping criteria.
int	fWeightMethod
	What method to use to weight hits.
double	fOtherViewW0
	Weight for NN hits in other view.
double	fNNZrange
	Range in z for NN weight.
double	fNNZedge
	Width of edge for NN weight.
double	fNNTrange
	Rnage in x or y for NN weight.
double	fNNTedge
	Width of edge for NN weight.
TH1F *	fFinalVtxX
	X vertex resolution.
TH1F *	fFinalVtxY
	Y vertex resolution.
TH1F *	fFinalVtxZ
	Z vertex resolution.

trk::CirceFit Class Reference

Public Member Functions

Private Attributes

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation