ana::SBNOsc Namespace Reference

Classes
class	Chi2Sensitivity
class	Covariance
class	NueSelection
	Electron neutrino event selection. More...
class	NumuSelection
	Electron neutrino event selection. More...
struct	VisibleEnergyCalculator
class	AuxDetSimChannelPrinter
	Electron neutrino event selection. More...
class	CaloPrinter
	Electron neutrino event selection. More...
class	CosmicArrivalTimes
	Electron neutrino event selection. More...
class	CRTHitTiming
	Electron neutrino event selection. More...
class	CRTSimHitViewer
	Electron neutrino event selection. More...
class	GetFV
	Electron neutrino event selection. More...
class	MCParticleTreePrinter
	Electron neutrino event selection. More...
class	OpDetBackTrackerPrinter
	Electron neutrino event selection. More...
class	OpDetWaveformMaker
	Electron neutrino event selection. More...
class	OpSimHitPrinter
	Electron neutrino event selection. More...
class	PandoraIDPrinter
	Electron neutrino event selection. More...
class	PandoraMetadataPrinter
	Electron neutrino event selection. More...
class	SimPhotonPrinter
	Electron neutrino event selection. More...
struct	CosmicHistos
struct	CRTHistos
class	HistoList
struct	InteractionHistos
struct	TrackProfiles
struct	TrackHistos
class	Cuts
class	Flatten
struct	Histograms
class	Normalize
struct	ROC
class	Selection
class	NumuReco
	Electron neutrino event selection. More...

Functions
double	osc_factor_L_integrated (double energy, double l_min, double l_max, double dm2)
double	numu_to_numu (double x, double sin, double dm2)
double	numu_to_nue (double x, double sin, double dm2)
double	NumuOscillate (ROOT::Math::IntegratorOneDim &integrator, double l_min, double l_max, double e_min, double e_max, double dm2, double sinth)
std::vector< float >	GetUniWeights (const std::map< std::string, std::vector< float > > &weights, const std::vector< std::string > &keys, int n_unis, int uni_offset)
double	aaBoxesMin (const std::vector< geoalgo::AABox > &boxes, unsigned dim)
double	aaBoxesMax (const std::vector< geoalgo::AABox > &boxes, unsigned dim)
geoalgo::AABox	shaveVolume (const geoalgo::AABox &select_volume, double delta)
void	hello ()
event::Interaction	TruthReco (const simb::MCTruth &mctruth)
double	ECCQE (const TVector3 &l_momentum, double l_energy, double energy_distortion, double angle_distortion)
double	NuMuOscillation (double numu_energy, double numu_dist, double osc_dm2, double osc_angle)
double	containedLength (const TVector3 &v0, const TVector3 &v1, const std::vector< geoalgo::AABox > &boxes)
double	visibleEnergyProposalMCParticles (TRandom &rand, const simb::MCTruth &mctruth, const std::vector< sim::MCTrack > mctrack_list, const VisibleEnergyCalculator &calculator)
double	visibleEnergyProposal (TRandom &rand, const simb::MCTruth &mctruth, const std::vector< sim::MCTrack > &mctrack_list, const VisibleEnergyCalculator &calculator)
double	visibleEnergy (TRandom &rand, const simb::MCTruth &mctruth, const std::vector< sim::MCTrack > &mctrack_list, const std::vector< sim::MCShower > &mcshower_list, const VisibleEnergyCalculator &calculator, bool include_showers)
double	smearLeptonEnergy (TRandom &rand, const sim::MCTrack &mct, const VisibleEnergyCalculator &calculator)
static const TDatabasePDG *	PDGTable (new TDatabasePDG)
double	PDGMass (int pdg)
double	PDGCharge (int pdg)
bool	isFromNuVertex (const simb::MCTruth &mc, const simb::MCParticle &mcp, float distance)
bool	isFromNuVertex (const simb::MCTruth &mc, const sim::MCShower &show, float distance)
bool	isFromNuVertex (const simb::MCTruth &mc, const sim::MCTrack &track, float distance)
double	closestDistance (const TVector3 &line0, const TVector3 &line1, const TVector3 &p)
double	closestDistanceDim (const TVector3 &line0, const TVector3 &line1, const TVector3 &p, int dim)
template<typename T , std::size_t N1, std::size_t N2>
constexpr std::array< T, N1+N2 >	concat (std::array< T, N1 > lhs, std::array< T, N2 > rhs)
void	SetEvent (numu::RecoEvent &event, const event::Event &core, const ana::SBNOsc::Cuts &cuts, numu::MCType file_type, bool use_calorimetry=true)
void	DumpTrueStart (const gallery::Event &ev, int mcparticle_id)
numu::Wall	GetWallCross (const geo::BoxBoundedGeo &volume, const TVector3 p0, const TVector3 p1)
sbn::crt::CRTHit	ICARUS2SBNDCrtHit (const sbn::crt::CRTHit &inp)
double	RecoTrackLength (const art::Ptr< recob::Track > &track)

Variables
static const TVector3	InvalidTVector3 = TVector3(-999, -999, -999)

Function Documentation

double ana::SBNOsc::aaBoxesMax	(	const std::vector< geoalgo::AABox > &	boxes,
		unsigned	dim
	)

Definition at line 65 of file NumuSelection.cxx.

                                                                      {
  return std::max_element(boxes.begin(), boxes.end(), [dim](auto &lhs, auto &rhs) { return lhs.Max()[dim] < rhs.Max()[dim]; })->Max()[dim];
}

double ana::SBNOsc::aaBoxesMin	(	const std::vector< geoalgo::AABox > &	boxes,
		unsigned	dim
	)

Definition at line 61 of file NumuSelection.cxx.

                                                                      {
  return std::min_element(boxes.begin(), boxes.end(), [dim](auto &lhs, auto &rhs) { return lhs.Min()[dim] < rhs.Min()[dim]; })->Min()[dim];
}

double ana::SBNOsc::closestDistance	(	const TVector3 &	line0,
		const TVector3 &	line1,
		const TVector3 &	p
	)

Definition at line 490 of file sbnana/sbnanalysis/ana/SBNOsc/Utilities.cxx.

                                                                                        {
  // Return minimum distance between line segment vw and point p
  //   const float l2 = length_squared(v, w);  // i.e. |w-v|^2 -  avoid a sqrt
  double length2 = (line0 - line1).Mag2();
  if (length2 == 0.0) return (line0 - p).Mag();
  // Consider the line extending the segment, parameterized as v + t (w - v).
  // We find projection of point p onto the line. 
  // It falls where t = [(p-v) . (w-v)] / |w-v|^2
  // We clamp t from [0,1] to handle points outside the segment vw.
  double t = std::max(0., std::min(1., (p - line0).Dot(line1 - line0) / length2));
  TVector3 projection = line0 + t * (line1 - line0);
  return (projection - p).Mag();
}

double ana::SBNOsc::closestDistanceDim	(	const TVector3 &	line0,
		const TVector3 &	line1,
		const TVector3 &	p,
		int	dim
	)

Definition at line 504 of file sbnana/sbnanalysis/ana/SBNOsc/Utilities.cxx.

                                                                                                    {
  double l0 = line0[dim];
  double l1 = line1[dim];
  double pp = p[dim];

  if ((l0 < pp && pp < l1) || (l1 < pp && pp < l0)) return 0.;
  return std::min(abs(l0 - pp) , abs(l1 - pp)); 
}

template<typename T , std::size_t N1, std::size_t N2>

constexpr std::array<T, N1+N2> ana::SBNOsc::concat	(	std::array< T, N1 >	lhs,
		std::array< T, N2 >	rhs
	)

Definition at line 28 of file Histograms.h.

                                                                                {
    std::array<T,N1+N2> result{};
    size_t index = 0;
    for (auto &el : lhs)  {
      result[index] = std::move(el);
      index ++;
    }
    for (auto &el : rhs) {
      result[index] = std::move(el);
      index ++;
    }
    return result;
  }

double ana::SBNOsc::containedLength	(	const TVector3 &	v0,
		const TVector3 &	v1,
		const std::vector< geoalgo::AABox > &	boxes
	)

Finds length of line segment contained inside AABox. Make sure that AABox and TVector's use the same units.

Parameters

v0	the first point of the line segment
v1	the second point of the line segment
boxes	a list of fiducial volumes instantiated as AABoxes

Returns

Length of line segment contained in the list of AABox's.

Definition at line 92 of file sbnana/sbnanalysis/ana/SBNOsc/Utilities.cxx.

                                                             {
  static const geoalgo::GeoAlgo algo;

  // if points are the same, return 0
  if ((v0 - v1).Mag() < 1e-6) return 0;

  // construct individual points
  geoalgo::Point_t p0(v0);
  geoalgo::Point_t p1(v1);

  // construct line segment
  geoalgo::LineSegment line(p0, p1);

  double length = 0;

  // total contained length is sum of lengths in all boxes
  // assuming they are non-overlapping
  for (auto const &box: boxes) {
    int n_contained = box.Contain(p0) + box.Contain(p1);
    // both points contained -- length is total length (also can break out of loop)
    if (n_contained == 2) {
      length = (v1 - v0).Mag();
      break;
    }
    // one contained -- have to find intersection point (which must exist)
    if (n_contained == 1) {
      auto intersections = algo.Intersection(line, box);
      // Because of floating point errors, it can sometimes happen
      // that there is 1 contained point but no "Intersections"
      // if one of the points is right on the edge
      if (intersections.size() == 0) {
        // determine which point is on the edge
        double tol = 1e-5;
        bool p0_edge = algo.SqDist(p0, box) < tol;
        bool p1_edge = algo.SqDist(p1, box) < tol;
        assert(p0_edge || p1_edge);
        // contained one is on edge -- can treat both as not contained
        //
        // In this case, no length
        if ((p0_edge && box.Contain(p0)) || (box.Contain(p1) && p1_edge))
          continue;
        // un-contaned one is on edge -- treat both as contained
        else if ((p0_edge && box.Contain(p1)) || (box.Contain(p0) && p1_edge)) {
          length = (v1 - v0).Mag();
          break;
        }
        else {
          assert(false); // bad
        }
      }
      // floating point errors can also falsely cause 2 intersection points
      //
      // in this case, one of the intersections must be very close to the 
      // "contained" point, so the total contained length will be about
      // the same as the distance between the two intersection points
      else if (intersections.size() == 2) {
        length += (intersections.at(0).ToTLorentzVector().Vect() - intersections.at(1).ToTLorentzVector().Vect()).Mag();
        continue;
      }
      // "Correct"/ideal case -- 1 intersection point
      else if (intersections.size() == 1) {
        // get TVector at intersection point
        TVector3 int_tv(intersections.at(0).ToTLorentzVector().Vect());
        length += ( box.Contain(p0) ? (v0 - int_tv).Mag() : (v1 - int_tv).Mag() ); 
      }
      else assert(false); // bad
    }
    // none contained -- either must have zero or two intersections
    if (n_contained == 0) {
      auto intersections = algo.Intersection(line, box);
      if (!(intersections.size() == 0 || intersections.size() == 2)) {
        // more floating point error fixes...
        //
        // figure out which points are near the edge
        double tol = 1e-5;
        bool p0_edge = algo.SqDist(p0, box) < tol;
        bool p1_edge = algo.SqDist(p1, box) < tol;
        // and which points are near the intersection
        TVector3 vint = intersections.at(0).ToTLorentzVector().Vect();
        bool p0_int = (v0 - vint).Mag() < tol;
        bool p1_int = (v1 - vint).Mag() < tol;
        // exactly one of them should produce the intersection
        assert((p0_int && p0_edge) != (p1_int && p1_edge));
        // both close to edge -- full length is contained
        if (p0_edge && p1_edge) {
          length += (v0 - v1).Mag();
        }
        // otherwise -- one of them is not on an edge, no length is contained
        else {}
      }
      // assert(intersections.size() == 0 || intersections.size() == 2);
      else if (intersections.size() == 2) {
        TVector3 start(intersections.at(0).ToTLorentzVector().Vect());
        TVector3 end(intersections.at(1).ToTLorentzVector().Vect());
        length += (start - end).Mag();
      }
    }
  }

  return length;
}

void ana::SBNOsc::DumpTrueStart	(	const gallery::Event &	ev,
		int	mcparticle_id
	)

Definition at line 67 of file NumuReco.cxx.

                                                              {
  // track.match.mcparticle_id);
  const std::vector<simb::MCParticle> &mcparticle_list = *ev.getValidHandle<std::vector<simb::MCParticle>>("largeant");
      std::cout << "Track: " << mcparticle_id;
  for (const simb::MCParticle &part: mcparticle_list) {
    if (mcparticle_id == part.TrackId()) {
      std::cout << " start T: " << part.Position().T() << " to: " << part.EndPosition().T() << std::endl;
    }
  }
  return;
}

double ana::SBNOsc::ECCQE	(	const TVector3 &	l_momentum,
		double	l_energy,
		double	energy_distortion = 0.,
		double	angle_distortion = 0.
	)

Calculate CCQE energy from associated lepton information (and optional distortion). Energy in GeV.

Parameters

l_momentum	Lepton momentum (in any units – used only to get angle info)
l_energy	Lepton energy in GeV
energy_distortion	Optional energy distortion in GeV
angle_distortion	Optiona langle distortion

Returns

CCQE energy in GeV.

Definition at line 59 of file sbnana/sbnanalysis/ana/SBNOsc/Utilities.cxx.

                                                                {
  // Based on D. Kaleko, LowEnergyExcess LArLite module ECCQECalculator
  double M_n = 0.939565; // GeV/c^2
  double M_p = 0.938272; // GeV/c^2
  double M_e = 0.000511; // GeV/c^2
  double bindingE = 0.0300; // GeV

  double mp2 = M_p * M_p;
  double me2 = M_e * M_e;
  double mnb = M_n - bindingE;

  // mess around with lorentz vector
  TVector3 v(l_momentum);
  v.SetTheta( v.Theta() + angle_distortion);
  l_energy = l_energy + energy_distortion;
  double l_mom = sqrt(l_energy*l_energy - me2);
  double l_theta = \
    acos(v.Pz() / sqrt(v.Px()*v.Px() + v.Py()*v.Py() + v.Pz()*v.Pz()));
  double enu_top = mp2 - mnb*mnb - me2 + 2.0 * mnb * l_energy;
  double enu_bot = 2.0 * (mnb - l_energy + l_mom * cos(l_theta));
  return enu_top / enu_bot;
}

std::vector<float> ana::SBNOsc::GetUniWeights	(	const std::map< std::string, std::vector< float > > &	weights,
		const std::vector< std::string > &	keys,
		int	n_unis,
		int	uni_offset
	)

Definition at line 61 of file Covariance.cxx.

                                                                                                                                                     {
    std::vector <float> uweights;
    for (int u = 0; u < n_unis; u++) {
        float weight = 1.;
        for (auto const &key: keys) {
            const std::vector<float>& this_weights = weights.at(key);
            int wind = u + uni_offset;
            weight *= this_weights.at(wind);
        }
        uweights.push_back(weight);
    }
    return uweights;
}

numu::Wall ana::SBNOsc::GetWallCross	(	const geo::BoxBoundedGeo &	volume,
		const TVector3	p0,
		const TVector3	p1
	)

Definition at line 79 of file NumuReco.cxx.

                                                                                            {
  TVector3 direction = (p1 - p0) * ( 1. / (p1 - p0).Mag());
  std::vector<TVector3> intersections = volume.GetIntersections(p0, direction);
  
  /*
  std::cout << "p0: " << p0.X() << " " << p0.Y() << " " << p0.Z() << std::endl;
  std::cout << "p1: " << p1.X() << " " << p1.Y() << " " << p1.Z() << std::endl;
  std::cout << "i0: " << intersections[0].X() << " " << intersections[0].Y() << " " << intersections[0].Z() << std::endl; 
  std::cout << "i1: " << intersections[1].X() << " " << intersections[1].Y() << " " << intersections[1].Z() << std::endl; 
  */

  assert(intersections.size() == 2);

  // get the intersection point closer to p0
  int intersection_i = ((intersections[0] - p0).Mag() < (intersections[1] - p0).Mag()) ? 0 : 1;
  

  double eps = 1e-3;
  if (abs(intersections[intersection_i].X() - volume.MinX()) < eps) {
    //std::cout << "Left\n";
    return numu::wLeft;
  }
  else if (abs(intersections[intersection_i].X() - volume.MaxX()) < eps) {
    //std::cout << "Right\n";
    return numu::wRight;
  }
  else if (abs(intersections[intersection_i].Y() - volume.MinY()) < eps) {
    //std::cout << "Bottom\n";
    return numu::wBottom;
  }
  else if (abs(intersections[intersection_i].Y() - volume.MaxY()) < eps) {
    //std::cout << "Top\n";
    return numu::wTop;
  }
  else if (abs(intersections[intersection_i].Z() - volume.MinZ()) < eps) {
    //std::cout << "Front\n";
    return numu::wFront;
  }
  else if (abs(intersections[intersection_i].Z() - volume.MaxZ()) < eps) {
    //std::cout << "Back\n";
    return numu::wBack;
  }
  else assert(false);
  //std::cout << "None\n";  

  return numu::wNone;
}

void ana::SBNOsc::hello

(

)

A function that says hello.

Definition at line 21 of file sbnana/sbnanalysis/ana/SBNOsc/Utilities.cxx.

             {
  std::cout << "Hello SBNOsc!" << std::endl;
}

sbn::crt::CRTHit ana::SBNOsc::ICARUS2SBNDCrtHit

(

const sbn::crt::CRTHit &

inp

)

Definition at line 138 of file NumuReco.cxx.

                                                        {
  sbn::crt::CRTHit ret;
  ret.feb_id = inp.feb_id;
  ret.pesmap = inp.pesmap;
  // convert ADC -> PE's
  // TODO: fix -- hardcoded for now as temporary hack
  unsigned n_strip = 2;
  double baseline = 63.6; // ADC
  double gain = 70; // ADC / PE
  ret.peshit = (inp.peshit - n_strip*baseline) / (gain * n_strip);
  ret.ts0_s = inp.ts0_s;
  ret.ts0_s_corr = inp.ts0_s_corr;
  ret.ts0_ns = inp.ts0_ns;
  ret.ts0_ns_corr = inp.ts0_ns_corr;
  ret.ts1_ns = inp.ts1_ns;
  ret.plane = inp.plane;
  ret.x_pos = inp.x_pos;
  ret.x_err = inp.x_err;
  ret.y_pos = inp.y_pos;
  ret.y_err = inp.y_err;
  ret.z_pos = inp.z_pos;
  ret.z_err = inp.z_err;
  ret.tagger = inp.tagger;
  return ret;
}

bool ana::SBNOsc::isFromNuVertex	(	const simb::MCTruth &	mc,
		const simb::MCParticle &	mcp,
		float	distance
	)

Definition at line 468 of file sbnana/sbnanalysis/ana/SBNOsc/Utilities.cxx.

                                                                                      {
  TVector3 nuVtx = mc.GetNeutrino().Nu().Position().Vect();
  TVector3 partStart = mcp.Position().Vect();
  return (partStart - nuVtx).Mag() < distance;
}

bool ana::SBNOsc::isFromNuVertex	(	const simb::MCTruth &	mc,
		const sim::MCShower &	show,
		float	distance = 5.0
	)

Returns whether track/shower object is from the neutrino vertex

Parameters

mc	MCTruth corresponding to neutrino interaction
show	The object to be matched
distance	between shower start and interaction vertex

Returns

Whether track/shower object is from neutrino vertex

Definition at line 474 of file sbnana/sbnanalysis/ana/SBNOsc/Utilities.cxx.

                                     {
  TVector3 nuVtx = mc.GetNeutrino().Nu().Trajectory().Position(0).Vect();
  TVector3 showStart = show.Start().Position().Vect();
  return (showStart - nuVtx).Mag() < distance;
}

bool ana::SBNOsc::isFromNuVertex	(	const simb::MCTruth &	mc,
		const sim::MCTrack &	track,
		float	distance = 5.0
	)

Returns whether track/shower object is from the neutrino vertex

Parameters

mc	MCTruth corresponding to neutrino interaction
track	The object to be matched
distance	between track start and interaction vertex

Returns

Whether track/shower object is from neutrino vertex

Definition at line 482 of file sbnana/sbnanalysis/ana/SBNOsc/Utilities.cxx.

                                    {
  TVector3 nuVtx = mc.GetNeutrino().Nu().Trajectory().Position(0).Vect();
  TVector3 trkStart = track.Start().Position().Vect();
  return (trkStart - nuVtx).Mag() < distance;
}

double ana::SBNOsc::numu_to_nue	(	double	x,
		double	sin,
		double	dm2
	)

Definition at line 42 of file Chi2Sensitivity.cxx.

                                                     {
    
    // x is L/E, sin is sin^2(2theta) and dm2 is delta m^2
    return sin * TMath::Power(TMath::Sin(1.27 * dm2 * x), 2);
}

double ana::SBNOsc::numu_to_numu	(	double	x,
		double	sin,
		double	dm2
	)

Definition at line 36 of file Chi2Sensitivity.cxx.

                                                      {
    
    // x is L/E, sin is sin^2(2theta) and dm2 is delta m^2
    return 1 - sin * TMath::Power(TMath::Sin(1.27 * dm2 * x), 2);
}

double ana::SBNOsc::NumuOscillate	(	ROOT::Math::IntegratorOneDim &	integrator,
		double	l_min,
		double	l_max,
		double	e_min,
		double	e_max,
		double	dm2,
		double	sinth
	)

Definition at line 48 of file Chi2Sensitivity.cxx.

                                                                                                                                             {
    auto integrand = [l_min, l_max, dm2](double energy) { return osc_factor_L_integrated(energy, l_min, l_max, dm2); };
    double integral = integrator.Integral(integrand, e_min, e_max);
    double factor = (1. / 2. ) + integral / ( (l_max - l_min)*(e_max - e_min));
    return 1 - sinth * factor;
}

double ana::SBNOsc::NuMuOscillation	(	double	numu_energy,
		double	numu_dis,
		double	osc_dm2,
		double	osc_angle
	)

Get oscillation probability of muon neutrino in a 3+1 model. I.e. probability that the numu will stay a numu.

Parameters

numu_energy	Energy of incident muon neutrino in GeV
numu_dist	Distance travelled by muon neutrino in km
osc_dm2	dm^2 of sterile netrino in eV^2
osc_angle	Sterile neutrino mixing angle

Returns

Probability of muon neutrino not oscillating in 3+1 model.

Definition at line 84 of file sbnana/sbnanalysis/ana/SBNOsc/Utilities.cxx.

                                                         {
  double overlap = sin(2*osc_angle);
  double energy_factor = sin(1.27 * osc_dm2 * numu_dist / numu_energy);
  return 1 - overlap * overlap * energy_factor * energy_factor;
}

double ana::SBNOsc::osc_factor_L_integrated	(	double	energy,
		double	l_min,
		double	l_max,
		double	dm2
	)

Definition at line 30 of file Chi2Sensitivity.cxx.

                                                                                      {
    double f = 1.27 * dm2;
    return (energy / (4 * f)) * (TMath::Sin(2 * f * l_min / energy) - TMath::Sin(2 * f * l_max / energy));
}

double ana::SBNOsc::PDGCharge

(

int

pdg

)

Get charge from PDGID of particle in |e|/3.

Parameters

pdg	The Particle Data Group ID of the particle (as returned by i.e. an MCTruth object)

Returns

Charge of particle in |e|/3.

Definition at line 455 of file sbnana/sbnanalysis/ana/SBNOsc/Utilities.cxx.

                          {
  // regular particle
  if (pdg < 1000000000) {
    TParticlePDG* ple = PDGTable->GetParticle(pdg);
    return ple->Charge();
  }
  // ion
  else {
    int p = (pdg % 10000000) / 10000;
    return p * 3;
  }
}

double ana::SBNOsc::PDGMass

(

int

pdg

)

Get mass from PDGID of particle in MeV/c^2.

Parameters

pdg	The Particle Data Group ID of the particle (as returned by i.e. an MCTruth object)

Returns

Mass of particle in MeV/c^2

Definition at line 438 of file sbnana/sbnanalysis/ana/SBNOsc/Utilities.cxx.

                        {
  // regular particle
  if (pdg < 1000000000) {
    TParticlePDG* ple = PDGTable->GetParticle(pdg);
    if (ple == NULL) return -1;
    return ple->Mass() * 1000.0;
  }
  // ion
  else {
    int p = (pdg % 10000000) / 10000;
    int n = (pdg % 10000) / 10 - p;
    return (PDGTable->GetParticle(2212)->Mass() * p +
            PDGTable->GetParticle(2112)->Mass() * n) * 1000.0;
  }
}

static const TDatabasePDG* ana::SBNOsc::PDGTable ( new  TDatabasePDG )

static

double ana::SBNOsc::RecoTrackLength

(

const art::Ptr< recob::Track > &

track

)

Definition at line 537 of file NumuReco.cxx.

                                                        {
  if (track->CountValidPoints() == 0) return 0.;
  double dist = 0.;
  geo::Point_t first = track->Start();
  for (size_t i = 1; i < track->CountValidPoints(); i++) {
    geo::Point_t second = track->LocationAtPoint(i);
    dist += sqrt((second - first).Mag2());
    first = second;
  }
  return dist;
}

void ana::SBNOsc::SetEvent	(	numu::RecoEvent &	event,
		const event::Event &	core,
		const ana::SBNOsc::Cuts &	cuts,
		numu::MCType	file_type,
		bool	use_calorimetry = true
	)

Definition at line 5 of file SetEvent.cc.

                                                                                                                                            {
  // update each reco Interaction to a smarter primary track selector
  unsigned i = 0;
  while (i < event.reco.size()) {
    for (size_t particle_ind: event.reco[i].slice.tracks) {
      numu::RecoTrack &track = event.tracks.at(particle_ind);

      // set containment
      if (cuts.InCalorimetricContainment(track.start) && cuts.InCalorimetricContainment(track.end)) {
        track.is_contained = true;
      }
      else {
        track.is_contained = false;
      }
    }

    int primary_track;
    if (use_calorimetry) {
      primary_track = numu::SelectLongestIDdMuon(event.tracks, event.reco[i].slice);
    }
    else {
      primary_track = numu::SelectLongestTrack(event.tracks, event.reco[i].slice);
    }

    // remove vertices without a good primary track
    if (primary_track < 0) {
      event.reco.erase(event.reco.begin() + i); 
      continue;
    }

    event.reco[i].slice.primary_track_index = primary_track;

    // re-do truth matching
    event.reco[i].match = numu::InteractionTruthMatch(core.truth, event.tracks, event.reco[i]);

    // if this is an in-time cosmic file, update the cosmic mode
    if (file_type == numu::fIntimeCosmic) {
      assert(event.reco[i].match.mode == numu::mCosmic || event.reco[i].match.mode == numu::mOther);
      event.reco[i].match.mode = numu::mIntimeCosmic;
    }

    i += 1;
  }
  // make sure no two vertices match to the same true neutrino interaction
  numu::CorrectMultiMatches(event, event.reco);

  // re-do the truth matching one more time with the multi-matches fixed
  for (unsigned i = 0; i < event.reco.size(); i++) {
    event.reco[i].match = numu::InteractionTruthMatch(core.truth, event.tracks, event.reco[i]);
  }
}

geoalgo::AABox ana::SBNOsc::shaveVolume	(	const geoalgo::AABox &	select_volume,
		double	delta
	)

Definition at line 69 of file NumuSelection.cxx.

                                                                          {

  double xmin = select_volume.Min()[0] + 5.;
  double ymin = select_volume.Min()[1] + 5.;
  double zmin = select_volume.Min()[2] + 5.;

  double xmax = select_volume.Max()[0] - 5.;
  double ymax = select_volume.Max()[1] - 5.;
  double zmax = select_volume.Max()[2] - 5.;

  return geoalgo::AABox(xmin, ymin, zmin, xmax, ymax, zmax);
}

double ana::SBNOsc::smearLeptonEnergy	(	TRandom &	rand,
		const sim::MCTrack &	mct,
		const VisibleEnergyCalculator &	calculator = VisibleEnergyCalculator()
	)

Get the smeared energy from a lepton.

Parameters

mctrack	The MCTrack object corresponding to the lepton
calculator	Struct containing values to be used in energy calculation

Definition at line 407 of file sbnana/sbnanalysis/ana/SBNOsc/Utilities.cxx.

                                                                                                          {
  // if not contained and zero length, return
  if (!calculator.lepton_contained && calculator.lepton_contained_length < 1e-4) {
    // std::cout << "Out of detector lepton\n";
    return 0;
  }

  double smearing_percentage;
  if (calculator.lepton_contained) {
    smearing_percentage = calculator.lepton_energy_distortion_contained;
  } else {
    double A = calculator.lepton_energy_distortion_leaving_A;
    double B = calculator.lepton_energy_distortion_leaving_B;
    smearing_percentage = -A * TMath::Log(B * calculator.lepton_contained_length);
  }
  // smear visible energy
  double lepton_visible_energy = (mct.Start().E()) / 1000.; /* MeV to GeV */
  double smeared_lepton_visible_energy = rand.Gaus(lepton_visible_energy, smearing_percentage * lepton_visible_energy);
  // clamp to 0
  smeared_lepton_visible_energy = std::max(smeared_lepton_visible_energy, 0.);

  // std::cout << "Lepton -- is_contained: " << calculator.lepton_contained << " length: " << calculator.lepton_contained_length << " smearing: " << smearing_percentage << " true E: " << lepton_visible_energy << " smeared E: " << smeared_lepton_visible_energy << std::endl; 

  return smeared_lepton_visible_energy;
}

event::Interaction ana::SBNOsc::TruthReco

(

const simb::MCTruth &

mctruth

)

Extract truth information to approximate reconstruction.

Definition at line 26 of file sbnana/sbnanalysis/ana/SBNOsc/Utilities.cxx.

                                                       {
  event::Interaction interaction;

  // Neutrino
  const simb::MCNeutrino& nu = mctruth.GetNeutrino();
  interaction.neutrino.energy = nu.Nu().EndMomentum().Energy();

  // Primary lepton
  const simb::MCParticle& lepton = nu.Lepton();
  interaction.lepton.pdg = lepton.PdgCode();
  interaction.lepton.energy = lepton.Momentum(0).Energy();
  interaction.lepton.momentum = lepton.Momentum(0).Vect();

  // Hadronic system
  for (int iparticle=0; iparticle<interaction.finalstate.size(); iparticle++) {
    event::FinalStateParticle fsp;
    const simb::MCParticle& particle = mctruth.GetParticle(iparticle);

    if (particle.Process() != "primary") {
      continue;
    }

    fsp.pdg = particle.PdgCode();
    fsp.energy = particle.Momentum(0).Energy();
    fsp.momentum = particle.Momentum(0).Vect();

    interaction.finalstate.push_back(fsp);
  }

  return interaction;
}

double ana::SBNOsc::visibleEnergy	(	TRandom &	rand,
		const simb::MCTruth &	mctruth,
		const std::vector< sim::MCTrack > &	mctrack_list,
		const std::vector< sim::MCShower > &	mcshower_list,
		const VisibleEnergyCalculator &	calculator = VisibleEnergyCalculator(),
		bool	include_showers = true
	)

Get the "visible" energy from a neutrino interaction. Is equal to sum of non-neutral hadronic kinetic energies and lepton total energies.

Parameters

ev	The gallery event.
mctruth	The MCTruth object corresponding to the interaction.
mctrack_list	Vector of MCTrack objects in the gallery event.
mcshower_list	Vector of MCShower objects in the gallery event.
calculator	Struct containing values to be used in energy calculation
smeared_lepton_energy	lepton energy to be used in calculation – will default to smearLeptonEnergy(mctruth, calculator) if not set

Returns

Visble energy in GeV.

Definition at line 330 of file sbnana/sbnanalysis/ana/SBNOsc/Utilities.cxx.

                                                                     {
  double visible_E = 0;

  // primary leptron track
  const sim::MCTrack *lepton_track = NULL;
  bool lepton_track_exists = false;

  // total up visible energy from tracks...
  double track_visible_energy = 0.;
  for (unsigned ind = 0; ind < mctrack_list.size(); ind++) {
    auto const &mct = mctrack_list[ind];
    // ignore particles not from nu vertex, non primary particles, and uncharged particles
    if (!isFromNuVertex(mctruth, mct) || abs(PDGCharge(mct.PdgCode())) < 1e-4 || mct.Process() != "primary")
       continue;
    // account for primary lepton later
    if ((abs(mct.PdgCode()) == 13 || abs(mct.PdgCode()) == 11) && calculator.lepton_index == ind) {
      continue;
    }
    // ignore non-primary pion??
    //if ((abs(mct.PdgCode()) == 211) && calculator.lepton_index != ind) {
    //  continue;
    //}

    double mass = PDGMass(mct.PdgCode());
    double this_visible_energy = (mct.Start().E() - mass) / 1000. /* MeV to GeV */;
    if (this_visible_energy > calculator.track_threshold) {
      track_visible_energy += this_visible_energy;
    }
  }

  // do energy smearing
  if (calculator.track_energy_distortion > 1e-4) {
    track_visible_energy = rand.Gaus(track_visible_energy, track_visible_energy*calculator.track_energy_distortion);
    // clamp to 0
    track_visible_energy = std::max(track_visible_energy, 0.);
  }

  // ...and showers
  double shower_visible_energy = 0.;
  if (include_showers) {
    for (auto const &mcs: mcshower_list) {
      // ignore particles not from nu vertex, non primary particles, and uncharged particles
      if (!isFromNuVertex(mctruth, mcs) || abs(PDGCharge(mcs.PdgCode())) < 1e-4 || mcs.Process() != "primary")
        continue; 
      // account for primary lepton later
      if ((abs(mcs.PdgCode()) == 13 || abs(mcs.PdgCode()) == 11) && isFromNuVertex(mctruth, mcs))
        continue; 

      double mass = PDGMass(mcs.PdgCode());
      double this_visible_energy = (mcs.Start().E() - mass) / 1000. /* MeV to GeV */;
      
      if (this_visible_energy > calculator.shower_threshold) {
        shower_visible_energy += this_visible_energy;
      }
    }
  }

  // do energy smearing
  if (calculator.shower_energy_distortion > 1e-4) {
    shower_visible_energy = rand.Gaus(shower_visible_energy, shower_visible_energy*calculator.shower_energy_distortion);
    // clamp to 0
    shower_visible_energy = std::max(shower_visible_energy, 0.);
  }

  visible_E = track_visible_energy + shower_visible_energy;

  // ...and primary lepton energy (for CC events)
  // only add in extra here if identified "lepton" is actually a lepton
  if (calculator.lepton_index >= 0 && (abs(mctrack_list[calculator.lepton_index].PdgCode()) == 13 || abs(mctrack_list[calculator.lepton_index].PdgCode()) == 11)) {
    visible_E += smearLeptonEnergy(rand, mctrack_list[calculator.lepton_index], calculator);
  }

  return visible_E;
}

double ana::SBNOsc::visibleEnergyProposal	(	TRandom &	rand,
		const simb::MCTruth &	mctruth,
		const std::vector< sim::MCTrack > &	mctrack_list,
		const VisibleEnergyCalculator &	calculator
	)

Definition at line 262 of file sbnana/sbnanalysis/ana/SBNOsc/Utilities.cxx.

                                                                                                                                                              {
  double total = 0.;

  // std::cout << "\n\nNew Interaction" << std::endl;
  // std::cout << "True E: " << mctruth.GetNeutrino().Nu().E() << std::endl;
  for (unsigned ind = 0; ind < mctrack_list.size(); ind++) {
    auto const &mct = mctrack_list[ind];
    int pdg = mct.PdgCode();

    double track_energy = (mct.Start().E() - PDGMass(pdg)) / 1000. /* MeV -> GeV */;
    double smear_energy = track_energy;

    // In the proposal, different PDG codes got different energies

    // kaons did not have the mass subtracted
    if (abs(pdg) == 321) {
      track_energy = mct.Start().E() / 1000.;
      smear_energy = track_energy;
    }
    // pions did not have mass subtracted in the smearing
    if (abs(pdg) == 211) {
      smear_energy = mct.Start().E() / 1000.;
    }

    bool skip = false;

    // threshold -- only for protons
    if (abs(pdg) == 2212 && track_energy < calculator.track_threshold) {
      skip = true;
    } 

    // add in smeared energy
    double this_smeared_energy = rand.Gaus(track_energy, smear_energy * calculator.track_energy_distortion);
    // clamp to zero
    this_smeared_energy = std::max(this_smeared_energy, 0.);

    // Ignore particles not from nu vertex
    if (!isFromNuVertex(mctruth, mct) || mct.Process() != "primary") {
      skip = true;
    }
    // ignore everything except for protons and kaons and pions
    if (!(abs(pdg) == 2212 || abs(pdg) == 211 || abs(pdg) == 321)) {
      skip = true;
    }
    // account for primary track later
    if ((abs(pdg) == 13 || abs(pdg) == 11) && calculator.lepton_index == ind) {
      skip = true;
    }
    // double distance = (mctruth.GetNeutrino().Nu().Trajectory().Position(0).Vect() - mct.Start().Position().Vect()).Mag();
    // double l_distance = (mctruth.GetNeutrino().Nu().Trajectory().Position(0) - mct.Start().Position()).Mag();
    // std::cout << "New particle -- true E: " << track_energy << " smeared E: " << this_smeared_energy << " PDG: " << pdg << " is primary: " << (calculator.lepton_index == ind) << " is skipped: " << skip << " distance: " << distance << " lorentz distance: " << l_distance << " process: " << mct.Process() << std::endl; 
    if (skip) continue;

    total += this_smeared_energy;
  }

  // ...and primary lepton energy (for CC events)
  // only add in extra here if identified "lepton" is actually a lepton
  if (calculator.lepton_index >= 0 && (abs(mctrack_list[calculator.lepton_index].PdgCode()) == 13 || abs(mctrack_list[calculator.lepton_index].PdgCode()) == 11)) {
    double lepton_energy = smearLeptonEnergy(rand, mctrack_list[calculator.lepton_index], calculator);
    total += lepton_energy;
  }
  // std::cout << "Reco Energy: " << total << std::endl;

  return total;
}

double ana::SBNOsc::visibleEnergyProposalMCParticles	(	TRandom &	rand,
		const simb::MCTruth &	mctruth,
		const std::vector< sim::MCTrack >	mctrack_list,
		const VisibleEnergyCalculator &	calculator
	)

Definition at line 196 of file sbnana/sbnanalysis/ana/SBNOsc/Utilities.cxx.

                                                                                                                                                                        {
  double total = 0;
  for (int iparticle=0; iparticle<mctruth.NParticles(); iparticle++) {
    const simb::MCParticle& particle = mctruth.GetParticle(iparticle);

    int pdg = particle.PdgCode();

    // if there's no mass entry, we probably don't care about it
    if (PDGMass(pdg) < 0) continue;

    double track_energy = (particle.Momentum().E() - PDGMass(pdg) / 1000. /* MeV -> GeV */);
    double smear_energy = track_energy;

    // In the proposal, different PDG codes got different energies

    // kaons did not have the mass subtracted
    if (abs(pdg) == 321) {
      track_energy = particle.Momentum().E();
      smear_energy = track_energy;
    }
    // pions did not have mass subtracted in the smearing
    if (abs(pdg) == 211) {
      smear_energy = particle.Momentum().E();
    }

    bool skip = false;

    // threshold -- only for protons
    if (abs(pdg) == 2212 && track_energy < calculator.track_threshold) {
      skip = true;
    } 

    // add in smeared energy
    double this_smeared_energy = rand.Gaus(track_energy, smear_energy * calculator.track_energy_distortion);
    // clamp to zero
    this_smeared_energy = std::max(this_smeared_energy, 0.);

    // Ignore particles not from nu vertex
    if (!isFromNuVertex(mctruth, particle) || particle.Process() != "primary") {
      skip = true;
    }
    // ignore everything except for protons and kaons and pions
    if (!(abs(pdg) == 2212 || abs(pdg) == 211 || abs(pdg) == 321)) {
      skip = true;
    }
    // double distance = (mctruth.GetNeutrino().Nu().Trajectory().Position(0).Vect() - mct.Start().Position().Vect()).Mag();
    // double l_distance = (mctruth.GetNeutrino().Nu().Trajectory().Position(0) - mct.Start().Position()).Mag();
    // std::cout << "New particle -- true E: " << track_energy << " smeared E: " << this_smeared_energy << " PDG: " << pdg << " is primary: " << (calculator.lepton_index == ind) << " is skipped: " << skip << " distance: " << distance << " lorentz distance: " << l_distance << " process: " << mct.Process() << std::endl; 
    if (skip) continue;

    total += this_smeared_energy;
  }

  // ...and primary lepton energy (for CC events)
  // only add in extra here if identified "lepton" is actually a lepton
  if (calculator.lepton_index >= 0 && (abs(mctrack_list[calculator.lepton_index].PdgCode()) == 13 || abs(mctrack_list[calculator.lepton_index].PdgCode()) == 11)) {
    double lepton_energy = smearLeptonEnergy(rand, mctrack_list[calculator.lepton_index], calculator);
    total += lepton_energy;
  }
  // std::cout << "Reco Energy: " << total << std::endl;

  return total;
  
  }

Variable Documentation

const TVector3 ana::SBNOsc::InvalidTVector3 = TVector3(-999, -999, -999)

static

Definition at line 65 of file NumuReco.cxx.