$\nu_e$ PID More...

Namespaces
	dec
	Decision tree PID.

Classes
class	AlignLibToVtx
	Rewrite a library file with events aligned according to vertex. More...

class	ApplyFilter
	Filter according to FilterList status of largest slice. More...

class	BackgroundEnricher
	Filter to remove non-pizero interactions. More...

class	DistanceMap
	Calculate and cache electrostatic potential between cells. More...

class	EventSummary
	Simple representation of event for LEM use. More...

class	FindLEMMatches
	Module to find best LEM matches. More...

class	FindMatchesAlg

class	FindMatchesTranspose

struct	FlippedPotentials
	Collection of Potential objects with odd and/or even view flipped. More...

class	GenFromLib
	Reconstitute ART events from a LEM library. More...

struct	GetWorkArgs

class	Heads
	Mapping from a subset of the library to their best matches. More...

class	HeadsTranspose

class	LEM
	Module to find best LEM matches. More...

class	LEMAssociator

class	LEMClient
	Fetch LEM PID result from remote server. More...

class	LEMInput

class	LEMLibrary
	Provide access to the LEM library to all modules. More...

class	LEMMixer

class	LEMSummarizer
	Transform slices into the minimal description LEM needs to PID them. More...

struct	LEMWebSettings

class	LEMWorker
	This runs on a high-mem node and reconstitutes LEMSummary events. More...

class	LEMWorkerOutput
	Runs after LEM on the high-mem node, returns results to the server. More...

class	Library
	Collection of events for matching. More...

class	LibrarySummary
	Details of the library LEM matches were made against. More...

class	LiteHit
	Compressed hit info, basic component of LEM events. More...

class	MakeLibrary
	Summarize ART events into LEM libary format. More...

class	MakeMatches
	Calculate final PID variables from match details. More...

class	MakePID
	Calculate final PID variables from match details. More...

class	Match
	Information about a LEM match. More...

class	MatchableEvent
	Attach some information used in matching to an EventSummary. More...

class	MatchIndices
	List of indices of matched events. More...

class	MatchList
	Collection of MatchSummary objects. More...

struct	MatchSummary
	Simplified Match information, suitable for serialization. More...

class	MergeMatches
	Merge LEM match lists from split running. More...

class	PIDDetails
	Attach LEM-specific info to the base PID object. More...

class	PIDExtraVars
	Outputs of the MakePID module. More...

class	PIDInput
	Attach MatchList and other useful stuff. More...

class	Potential
	Map of electrostatic potential at each cell. More...

class	Preselection

struct	ReturnResultsArgs

class	SortedMatchList

class	UnTranspose

Typedefs
typedef art::Source< GenFromLib >	GenFromLibSource

Functions
double	SimpleSurvivalProb (const EventSummary &evt)

float	CalcEnergy (const Potential &Va, const MatchableEvent &a, const MatchableEvent &b)

double	FracChargeMatched (const EventSummary &a, const EventSummary &b, bool flipEven, bool flipOdd)

void	FillPotential (const EventSummary &trial, Potential &V, bool flipEven, bool flipOdd)

void	FillPotential (const EventSummary &trial, FlippedPotentials &Vs)

std::vector< Match >	FindMatches (const MatchableEvent &trial, unsigned int libSize, const MatchableEvent lib, const MatchableEvent libDownsample, int factor, const std::multiset< float > &alreadyInput, unsigned int numMatches, int enrich, const Heads *heads, int headIdx, bool flipHeadEven, bool flipHeadOdd, bool useDownsample)

std::ostream &	operator<< (std::ostream &os, const lem::Heads &hs)

std::ostream &	operator<< (std::ostream &os, const lem::HeadsTranspose &hs)

std::ostream &	operator<< (std::ostream &os, const lem::LibrarySummary &ls)

bool	compareByCell (const LiteHit &a, const LiteHit &b)

bool	compareByPlane (const LiteHit &a, const LiteHit &b)

std::ostream &	operator<< (std::ostream &os, const lem::LiteHit &h)

void	SetOnDemandMemoryLimit (long limit_mb)

void	CreateOnDemandMapping (const std::string &fname, const char *base)

long	MMapFileAtAddress (const std::string &fname, void *addr, bool touchAll, const std::string &lockName)

void *	GetWorkThreadFunc (void *x)

void *	ReturnResultsThreadFunc (void *x)

int	Clamp (int x)

int	MeanCell (geo::View_t view, const rb::Cluster &slice)

void	DefaultVertex (const rb::Cluster &slice, int &plane, int &cell, int &cellOtherView, bool reverse)

geo::CellUniqueId	FindCellRobust (double x, double y, double z)
	Steps around a bit. Can still return zero for total failure. More...

void	VertexToPlaneAndCell (const TVector3 vtx, const rb::Cluster &slice, int &plane, int &cell, int &cellOtherView, bool reverse)

size_t	writeMemoryCallback (void contents, size_t size, size_t nmemb, void userp)

std::string	get_http_response_internal (const std::string &host, const std::string &query, const std::string &postdata, int minPort, int maxPort, int queryTimeout, int retryTimeout, bool &ok)

std::string	get_query (const std::string &host, const std::string &query, int minPort, int maxPort, int queryTimeout, int retryTimeout, bool &ok)
	Make an HTTP GET query. More...

std::string	post_query (const std::string &host, const std::string &query, const std::string &postdata, int minPort, int maxPort, int queryTimeout, int retryTimeout, bool &ok)
	Make an HTTP POST query. More...

std::string	get_release ()
	Figure out the current release. TODO: does this belong somewhere else? More...

Variables
const int	kNumPlanes = 256

const int	kNumCells = 256

const int	kVertexPlane = 64

const int	kVertexCell = 128

const unsigned int	kMaxNumMatches = 1001

Detailed Description

$\nu_e$ PID

Provide access to the LEM library to all modules.

Author: Christopher Backhouse - bckho.nosp@m.use@.nosp@m.calte.nosp@m.ch.e.nosp@m.du

Typedef Documentation

typedef art::Source<GenFromLib> lem::GenFromLibSource

Definition at line 216 of file GenFromLib_source.cc.

Function Documentation

float lem::CalcEnergy	(	const Potential &	Va,
		const MatchableEvent &	a,
		const MatchableEvent &	b
	)

Definition at line 21 of file FindMatches.cxx.

References lem::LiteHit::cellIdx, energy, lem::EventSummary::hits, MECModelEnuComparisons::i, Nb, lem::EventSummary::nhits, lem::LiteHit::pecorr, lem::MatchableEvent::selfEnergy, and lem::Potential::V.

Referenced by lem::FindMatchesAlg::BestHeadFor(), lem::FindMatchesTranspose::endRun(), FindMatches(), lem::FindMatchesAlg::FindMatchesHeads(), lem::FindMatchesAlg::FindMatchesNoHeads(), and lem::LEM::produce().

   {
     float energy = a.selfEnergy + b.selfEnergy;
 
     const LiteHit* hb = b.hits;
     const unsigned int Nb = b.nhits;
 
     for(unsigned int i = 0; i < Nb; ++i){
       energy -= hb->pecorr*Va.V[hb->cellIdx];
       ++hb;
     }
 
     return energy;
   }

int lem::Clamp ( int x )

Definition at line 154 of file MakeLibrary_module.cc.

References cet::sqlite::max(), and min().

   {
     return std::min(std::max(0, x), 255);
   }

bool lem::compareByCell	(	const LiteHit &	a,
		const LiteHit &	b
	)

Definition at line 11 of file LiteHit.cxx.

References lem::LiteHit::Cell().

Referenced by lem::LiteHit::operator<().

   {
     return a.Cell() < b.Cell();
   }

bool lem::compareByPlane	(	const LiteHit &	a,
		const LiteHit &	b
	)

Definition at line 16 of file LiteHit.cxx.

References lem::LiteHit::Plane().

Referenced by lem::LiteHit::operator<().

   {
     return a.Plane() < b.Plane();
   }

void lem::CreateOnDemandMapping	(	const std::string &	fname,
		const char *	base
	)

Definition at line 123 of file OnDemand.cxx.

References ana::assert(), fclose(), and fileno.

Referenced by lem::Heads::FromMMapOnDemand(), lem::Library::FromMMapOnDemand(), and lem::HeadsTranspose::FromMMapOnDemand().

   {
     // Must be aligned on OS pages
     assert(long(base)%(4*1024) == 0);
 
     // Open the file to check it exists and find out its length
     FILE* f = fopen(fname.c_str(), "r");
     assert(f);
 
     struct stat st;
     fstat(fileno(f), &st);
     const long fsize = st.st_size;
     fclose(f);
 
 
     gMappings.push_back(Mapping(fname, base, fsize));
 
     // Set up signal handler to fault this file (and others) in on demand. It
     // should be OK to do this multiple times.
     struct sigaction sa;
     memset(&sa, 0, sizeof(sa));
     sigemptyset(&sa.sa_mask);
 
     sa.sa_sigaction = segfault_sigaction;
     sa.sa_flags = SA_SIGINFO;
 
     sigaction(SIGSEGV, &sa, NULL);
   }

void lem::DefaultVertex	(	const rb::Cluster &	slice,
		int &	plane,
		int &	cell,
		int &	cellOtherView,
		bool	reverse
	)

Definition at line 20 of file Util.cxx.

References geom(), rb::Cluster::MaxPlane(), MeanCell(), rb::Cluster::MinPlane(), geo::GeometryBase::NPlanes(), geo::GeometryBase::Plane(), POTSpillRate::view, and geo::PlaneGeo::View().

Referenced by lem::LEMSummarizer::produce(), lem::AlignLibToVtx::produce(), and VertexToPlaneAndCell().

   {
     art::ServiceHandle<geo::Geometry> geom;
 
     // Fallback when everything goes wrong.
 
     // Just put the vertex at the top of the slice in the middle
 
     plane = reverse ? slice.MaxPlane() : slice.MinPlane();
     if(plane < 0 || plane >= int(geom->NPlanes())) plane = 0; // arghhh
     const geo::View_t view = geom->Plane(plane)->View();
     cell = MeanCell(view, slice);
     cellOtherView = MeanCell(geo::View_t(1-view), slice);
   }

void lem::FillPotential	(	const EventSummary &	trial,
		Potential &	V,
		bool	flipEven,
		bool	flipOdd
	)

Definition at line 76 of file FindMatches.cxx.

References abs(), getBrightness::cell, lem::LiteHit::Cell(), make_syst_table_plots::h, lem::EventSummary::hits, MECModelEnuComparisons::i, lem::DistanceMap::Instance(), lem::DistanceMap::InvDist(), kVertexCell, lem::EventSummary::nhits, lem::LiteHit::pecorr, lem::LiteHit::Plane(), NDAPDHVSetting::plane, lem::Potential::V, and POTSpillRate::view.

Referenced by lem::FindLEMMatches::AddHists(), lem::LEM::AddHists(), FindMatches(), lem::FindMatchesAlg::FindMatchesHeads(), lem::FindMatchesAlg::FindMatchesNoHeads(), lem::FindMatchesTranspose::produce(), lem::LEM::produce(), and show_potential().

   {
     const DistanceMap& dm = DistanceMap::Instance();
 
     const unsigned int N = trial.nhits;
 
     for(unsigned int i = 0; i < N; ++i){
       const LiteHit& h = trial.hits[i];
 
       // Doesn't matter? Has to match in qfrac definition
       const int view = abs(h.Plane())%2;
 
       int hc = h.Cell();
       if((flipEven && view == 0) || (flipOdd && view == 1))
         hc = 2*lem::kVertexCell-hc;
 
       for(int plane = 0; plane < 256; ++plane){
         // Only in the same view
         if(abs(plane-h.Plane())%2) continue;
         for(int cell = 0; cell < 256; ++cell){
 
           const double contrib = h.pecorr*dm.InvDist(plane, h.Plane(),
                                                      cell, hc);
 
           V.V[256*plane+cell] += contrib;
         }
       }
     }
   }

void lem::FillPotential	(	const EventSummary &	trial,
		FlippedPotentials &	Vs
	)

Definition at line 108 of file FindMatches.cxx.

References abs(), getBrightness::cell, lem::LiteHit::Cell(), make_associated_cosmic_defs::fc, make_syst_table_plots::h, lem::EventSummary::hits, MECModelEnuComparisons::i, lem::DistanceMap::Instance(), lem::DistanceMap::InvDist(), kVertexCell, lem::EventSummary::nhits, lem::LiteHit::pecorr, lem::LiteHit::Plane(), NDAPDHVSetting::plane, lem::Potential::V, lem::FlippedPotentials::V, and POTSpillRate::view.

   {
     // We could just call FillPotential 4 times with each flip variant like so
     /*
     for(bool flipEven: {false, true}){
       for(bool flipOdd: {false, true}){
         FillPotential(trial, Vs.V[flipEven][flipOdd], flipEven, flipOdd);
       }
     }
     */
 
     // But this function actually takes noticeable runtime, so instead just
     // fill it once, and then copy into the over three flips.
     const DistanceMap& dm = DistanceMap::Instance();
 
     const unsigned int N = trial.nhits;
 
     // There's a complication because 256 is odd there isn't a perfect
     // centreline to flip round. The way things are defined, we end up needing
     // to know what cell 256 in the original map would have been, which we
     // don't have. Make a place to store it here.
     float overflow[256] = {0,};
 
     for(unsigned int i = 0; i < N; ++i){
       const LiteHit& h = trial.hits[i];
 
       for(int plane = 0; plane < 256; ++plane){
         // Only in the same view
         if(abs(plane-h.Plane())%2) continue;
         for(int cell = 0; cell < 256; ++cell){
 
           const double contrib = h.pecorr*dm.InvDist(plane, h.Plane(),
                                                      cell,  h.Cell());
 
           Vs.V[false][false].V[256*plane+cell] += contrib;
         }
         // And one for cell 256, stored seperately
         overflow[plane] += h.pecorr*dm.InvDist(plane, h.Plane(),
                                                256,   h.Cell());
       }
     }
 
     // Now copy over into the flipped variants
 
     for(bool flipEven: {false, true}){
       for(bool flipOdd: {false, true}){
         if(!flipEven && !flipOdd) continue;
 
         for(int plane = 0; plane < 256; ++plane){
           const int view = plane%2;
 
           for(int cell = 0; cell < 256; ++cell){
 
             int fc = cell; // flipped cell
             if((flipEven && view == 0) || (flipOdd && view == 1))
               fc = 2*lem::kVertexCell-fc;
 
             if(fc >= 0 && fc < 256)
               Vs.V[flipEven][flipOdd].V[plane*256+cell] = Vs.V[false][false].V[plane*256+fc];
             if(fc == 256) Vs.V[flipEven][flipOdd].V[plane*256+cell] = overflow[plane];
           }
         }
       }
     }
   }

geo::CellUniqueId lem::FindCellRobust	(	double	x,
		double	y,
		double	z
	)

Steps around a bit. Can still return zero for total failure.

Definition at line 38 of file Util.cxx.

References geo::GeometryBase::CellId(), and geom().

Referenced by VertexToPlaneAndCell().

   {
     art::ServiceHandle<geo::Geometry> geom;
 
     for(int trial = 0; trial < 4; ++trial){
       // Figure out what cell the vertex was in. Step 3cm in z in case of
       // failure. If that didn't work, it must be in plastic between cells,
       // move laterally in x or y.
       double offsetX = 0, offsetY = 0;
       if(trial == 1) offsetX = -2;
       if(trial == 2) offsetY = -2;
       if(trial == 3) offsetX = offsetY = -2;
 
       const geo::CellUniqueId id = geom->CellId(x+offsetX, y+offsetY, z,
                                                 0, 0, 1, 3);
       // We found one that works
       if(id) return id;
     }
     return 0;
   }

std::vector< Match > lem::FindMatches	(	const MatchableEvent &	trial,
		unsigned int	libSize,
		const MatchableEvent *	lib,
		const MatchableEvent *	libDownsample,
		int	factor,
		const std::multiset< float > &	alreadyInput,
		unsigned int	numMatches,
		int	enrich,
		const Heads *	heads,
		int	headIdx,
		bool	flipHeadEven,
		bool	flipHeadOdd,
		bool	useDownsample
	)

Definition at line 175 of file FindMatches.cxx.

References a, ana::assert(), CalcEnergy(), lem::Heads::ChildIdxFor(), lem::Potential::Downsampled(), E, FillPotential(), lem::Heads::FlipEvenFor(), lem::Heads::FlipOddFor(), FracChargeMatched(), lem::Heads::HeadSeqsLen(), MECModelEnuComparisons::i, calib::j, m, runNovaSAM::ret, lem::MatchableEvent::selfEnergy, lem::FlippedPotentials::V, lem::MatchableEvent::weight, and weight.

Referenced by lem::LEM::produce().

   {
     //    assert(libsize >= kMaxNumMatches);
 
     assert(!heads || headIdx >= 0);
 
     float target = FLT_MAX; // infinity
 
     // "De-normalize" values back into the range most of this function deals
     // with. Need to copy to a new container, since set iterators are const
     std::multiset<float> already;
     for(auto a: alreadyInput) already.insert(a*trial.selfEnergy);
 
     // Truncate to the right length. Update target if possible.
     while(already.size() > numMatches){
       auto itLargest = already.end();
       --itLargest;
       already.erase(itLargest);
       target = *already.rbegin();
     }
 
     // Precalculate potential
     // Try it with either/both of x&y flipped. Effectively gain factor of 4 library stats.
     // TODO: This isn't entirely right because of attenuation. It /would/
     // actually be possible to introduce this properly during library building
     // when truth information is available.
 
     FlippedPotentials Vs;
     FillPotential(trial, Vs);
 
     // Make downsampled versions of the four potential cases
     FlippedPotentials Vs_down;
     if(useDownsample){
       for(int flipEven = 0; flipEven <= 1; ++flipEven){
         for(int flipOdd = 0; flipOdd <= 1; ++flipOdd){
           Vs_down.V[flipEven][flipOdd] = Vs.V[flipEven][flipOdd].Downsampled(factor);
         }
       }
     }
 
     // Candidates for return, remains sorted
     std::multiset<Match> Elist;
 
     // This is the quickest ordering of the loops (locality into lib?)
     for(unsigned int j = 0; j < libSize; ++j){
       int i = j;
       if(heads){
         if(int(j) >= heads->HeadSeqsLen()) break;
         i = heads->ChildIdxFor(headIdx, j);
       }
 
       if(enrich != 2){
         if(lib[i].enrich != bool(enrich)) continue;
       }
 
       const MatchableEvent& downEvt = libDownsample[i];
 
       // There's beam nue in the library, but we don't want to match to it
       if( useDownsample && downEvt.weight == 0) continue;
       if(!useDownsample && lib[i].weight  == 0) continue;
 
       for(int flipEven = 0; flipEven <= 1; ++flipEven){
         if(heads && (heads->FlipEvenFor(headIdx, j) ^ flipHeadEven) != flipEven) continue;
         for(int flipOdd = 0; flipOdd <= 1; ++flipOdd){
           if(heads && (heads->FlipOddFor(headIdx, j) ^ flipHeadOdd) != flipOdd) continue;
 
           // Compare the downsampled event to the downsampled potential. The
           // resulting energy is a lower-bound on what the real match would
           // have done.
           float E = useDownsample ? CalcEnergy(Vs_down.V[flipEven][flipOdd], trial, downEvt) : 0;
 
           if(E < target){
             // Calculate the energy fully
             E = CalcEnergy(Vs.V[flipEven][flipOdd], trial, lib[i]);
 
             // If it's still good, use it
             if(E < target){
               Elist.insert(Match(E, &lib[i], flipEven, flipOdd, j));
               already.insert(E);
 
               if(Elist.size() > numMatches){
                 // Always drop the last one, to keep list correct size
                 auto itLargest = Elist.end();
                 --itLargest;
                 Elist.erase(itLargest);
                 // In case "already" isn't long enough
                 target = Elist.rbegin()->energy;
               }
               if(already.size() > numMatches){
                 // Always drop the last one, to keep list correct size
                 auto itLargest = already.end();
                 --itLargest;
                 already.erase(itLargest);
                 target = *already.rbegin();
               }
             }
           }
         } // end for flipOdd
       } // end for flipEven
     } // end for i
 
     // Convert back to a vector, and calculate the qfrac variable for surviving
     // matches.
     std::vector<Match> ret;
     ret.reserve(numMatches);
     for(const auto& m: Elist){
       ret.push_back(m);
       ret.back().qfrac = FracChargeMatched(trial, *m.evt, m.flipEven, m.flipOdd);
       ret.back().energy /= trial.selfEnergy; // Normalize to approx 0-1
     }
     return ret;
   }

double lem::FracChargeMatched	(	const EventSummary &	a,
		const EventSummary &	b,
		bool	flipEven,
		bool	flipOdd
	)

Definition at line 37 of file FindMatches.cxx.

References abs(), lem::LiteHit::Cell(), lem::LiteHit::cellIdx, lem::EventSummary::hits, kVertexCell, min(), Na, Nb, lem::EventSummary::nhits, lem::LiteHit::pecorr, lem::LiteHit::Plane(), lem::EventSummary::totalPE, and registry_explorer::v.

Referenced by lem::FindMatchesTranspose::endRun(), FindMatches(), and lem::FindMatchesAlg::PackageMatches().

   {
     // TODO: Duplicate hits definitely screw this up
 
     float qmatched = 0;
 
     // Quicker to spend this (not much) memory than the more obvious
     // nested-loops implementation.
     float acharge[256*256] = {0,};
 
     const unsigned int Na = a.nhits;
     const unsigned int Nb = b.nhits;
     for(unsigned int ia = 0; ia < Na; ++ia){
       const LiteHit& ha = a.hits[ia];
       acharge[ha.cellIdx] += ha.pecorr;
     } // end for ia
 
     for(unsigned int ib = 0; ib < Nb; ++ib){
       const LiteHit& hb = b.hits[ib];
 
       int hbc = hb.Cell();
 
       const int v = abs(hb.Plane())%2; // Matches definition in FillPotential
       if((v == 0 && flipEven) ||
          (v == 1 && flipOdd)) hbc = 2*lem::kVertexCell-hbc;
 
       const int bidx = 256*hb.Plane()+hbc;
 
       // This is the MINOS definition
       //      if(acharge[bidx]) qmatched += acharge[bidx]+hb.pecorr;
 
       qmatched += 2*std::min(acharge[bidx], hb.pecorr);
     } // end for ib
 
     return qmatched/(a.totalPE+b.totalPE);
   }

std::string lem::get_http_response_internal	(	const std::string &	host,
		const std::string &	query,
		const std::string &	postdata,
		int	minPort,
		int	maxPort,
		int	queryTimeout,
		int	retryTimeout,
		bool &	ok
	)

Definition at line 31 of file WebUtils.cxx.

References ana::assert(), om::cerr, om::cout, allTimeWatchdog::endl, genie::utils::style::Format(), central_limit::rand, sw, getGoodRuns4SAM::t0, getGoodRuns4SAM::t1, recentWatchdog::url, and writeMemoryCallback().

Referenced by get_query(), and post_query().

   {
     ok = false;
 
     const time_t t0 = time(0);
 
     // No, this is explicitly not thread-safe. Caller should call this themselves
     //    curl_global_init(CURL_GLOBAL_ALL);
     CURL* curl_handle = curl_easy_init();
     assert(curl_handle);
 
     // Don't send "Expect: 100-continue" for long queries, just forge
     // ahead. The chances of being rejected are pretty low and the server
     // doesn't send a continue anyway.
     curl_slist* slist = 0;
     slist = curl_slist_append(slist, "Expect:");
     curl_easy_setopt(curl_handle, CURLOPT_HTTPHEADER, slist);
 
     // Configure callback function
     std::string response;
     curl_easy_setopt(curl_handle, CURLOPT_WRITEFUNCTION, writeMemoryCallback);
     curl_easy_setopt(curl_handle, CURLOPT_WRITEDATA, &response);
 
     // Timeout for the HTTP connection. Will try again with a new connection up
     // to retryTimeout.
     curl_easy_setopt(curl_handle, CURLOPT_TIMEOUT, queryTimeout);
 
     // Work around this crash
     // http://stackoverflow.com/questions/9191668/error-longjmp-causes-uninitialized-stack-frame
     curl_easy_setopt(curl_handle, CURLOPT_NOSIGNAL, 1);
 
     // Allow redirection
     curl_easy_setopt(curl_handle, CURLOPT_FOLLOWLOCATION, 1);
 
     if(!postdata.empty()){
       curl_easy_setopt(curl_handle, CURLOPT_POST, 1);
       curl_easy_setopt(curl_handle, CURLOPT_POSTFIELDS, postdata.c_str());
     }
 
     double realTime = -1;
 
     int k = 0;
     while(true){
       const TString url = TString::Format("%s:%d/%s",
                                           host.c_str(),
                                           minPort+rand()%(maxPort-minPort+1),
                                           query.c_str());
       curl_easy_setopt(curl_handle, CURLOPT_URL, url.Data());
 
       std::cout << "Making query " << url;
       if(!postdata.empty()) std::cout << " with " << postdata.size() << " bytes of postdata";
       std::cout << std::endl;
 
       response.clear();
       TStopwatch sw;
       CURLcode code = curl_easy_perform(curl_handle);
       sw.Stop();
       realTime = sw.RealTime();
 
       if(code != CURLE_OK){
         std::cerr << "Connection failed: " << curl_easy_strerror(code)
                   << std::endl;
       }
       else {
         long http_code;
         curl_easy_getinfo(curl_handle, CURLINFO_RESPONSE_CODE, &http_code);
         if(http_code == 200 && code != CURLE_ABORTED_BY_CALLBACK){
           //Succeeded
           ok = true;
           break;
         }
         else{
           std::cerr << "Connection failed with HTTP " << http_code
                     << ": " << response << std::endl;
         }
       }
 
       // For subsequent iterations round the loop let's have a lot of debugging
       // output.
       curl_easy_setopt(curl_handle, CURLOPT_VERBOSE, 1);
 
       const int delay = 1 + ((double)random()/(double)RAND_MAX) * (1 << k++);
       const time_t t1 = time(0);
       if(t1+delay > t0+retryTimeout) break;
 
       std::cerr << "Iteration " << k
                 << ": will wait for " << delay << " sec and try again. "
                 << t1-t0 << " sec since we started. Will retry until "
                 << retryTimeout << " sec." << std::endl;
       sleep(delay);
     } // end while
 
     curl_easy_cleanup(curl_handle);
     // No, this is explicitly not thread-safe. Caller should call this themselves
     //    curl_global_cleanup();
 
     curl_slist_free_all(slist);
 
     if(!response.empty() && response.size() < 20){
       std::cout << "Response: " << response << " in " << realTime << " sec" << std::endl;
     }
     else{
       std::cout << "Response: " << response.size() << " bytes in " << realTime << " sec" << std::endl;
     }
 
     return response;
   }

std::string lem::get_query	(	const std::string &	host,
		const std::string &	query,
		int	minPort,
		int	maxPort,
		int	queryTimeout,
		int	retryTimeout,
		bool &	ok
	)

Make an HTTP GET query.

Definition at line 147 of file WebUtils.cxx.

References get_http_response_internal().

Referenced by GetWorkThreadFunc().

   {
     return get_http_response_internal(host, query, "",
                                       minPort, maxPort,
                                       queryTimeout, retryTimeout,
                                       ok);
   }

std::string lem::get_release ( )

Figure out the current release. TODO: does this belong somewhere else?

Definition at line 178 of file WebUtils.cxx.

References om::cerr, e, allTimeWatchdog::endl, and cet::getenv().

Referenced by lem::LEMWorker::GetWork(), and lem::LEMClient::produce().

   {
     char* e = getenv("NOVA_RELEASE");
     if(e) return e;
     e = getenv("SRT_BASE_RELEASE");
     if(e) return e;
     std::cerr << "Unable to determine current release. Aborting." << std::endl;
     abort();
   }

void* lem::GetWorkThreadFunc ( void * x )

Definition at line 101 of file LEMWorker_module.cc.

References runNovaSAM::args, om::cerr, allTimeWatchdog::endl, get_query(), lem::LEMWebSettings::host, lem::LEMWebSettings::maxPort, lem::LEMWebSettings::minPort, lem::GetWorkArgs::params, lem::LEMWebSettings::queryTimeout, lem::LEMWebSettings::retryTimeout, and lem::GetWorkArgs::web.

Referenced by lem::LEMWorker::GetWork().

   {
     const GetWorkArgs* args = (GetWorkArgs*)x;
 
     bool ok;
     const std::string work = get_query(args->web.host(),
                                        args->params,
                                        args->web.minPort(),
                                        args->web.maxPort(),
                                        args->web.queryTimeout(),
                                        args->web.retryTimeout(),
                                        ok);
     if(!ok){
       std::cerr << "No response from server after many retries. "
                 << "Maybe there's no work for us. Aborting"
                 << std::endl;
       abort();
     }
 
     delete args;
     return new std::string(work);
   }

int lem::MeanCell	(	geo::View_t	view,
		const rb::Cluster &	slice
	)

Definition at line 14 of file Util.cxx.

References rb::Cluster::MaxCell(), and rb::Cluster::MinCell().

Referenced by DefaultVertex().

   {
     return (slice.MinCell(view)+slice.MaxCell(view))/2;
   }

long lem::MMapFileAtAddress	(	const std::string &	fname,
		void *	addr,
		bool	touchAll,
		const std::string &	lockName = `""`
	)

Parameters

fname	File to map in
addr	Where to map it
touchAll	Read from every page to make sure it's cached
lockName	If touchAll, coordinate using lock here

Returns: Size of mapping in bytes

Definition at line 18 of file Util.cxx.

References ana::assert(), om::cout, allTimeWatchdog::endl, MakeMiniprodValidationCuts::f, fclose(), fileno, MECModelEnuComparisons::i, makeTrainCVSamples::int, runNovaSAM::ret, and PandAna.Demos.pi0_spectra::tot.

Referenced by lem::Heads::FromMMap(), lem::Library::FromMMap(), and lem::HeadsTranspose::FromMMap().

   {
     // Touching all the memory ("preloading") means we can be sure that it's
     // all quickly available to us later. The first process to do this has to
     // fetch the whole file from disk. We don't want more than one process
     // doing that, it's counter-productive. But once that's happened once, it's
     // in disk cache, and all other processes should feel free to preload in
     // parallel. This locking scheme ensures that.
 
     FILE* lockMain = 0;
     FILE* lockDone = 0;
     if(touchAll && !lockName.empty()){
       // No two processes should be preloading from disk at once
 
       // Make sure not to block out other users from this lock
       chmod(lockName.c_str(), S_IRUSR | S_IWUSR | S_IRGRP | S_IWGRP | S_IROTH);
       chmod((lockName+".done").c_str(), S_IRUSR | S_IWUSR | S_IRGRP | S_IWGRP | S_IROTH);
 
       lockMain = fopen(lockName.c_str(), "w");
       int ret = flock(fileno(lockMain), LOCK_EX);
       assert(ret == 0); // Blocking flock call, must have got it
 
       // We're now holding the loading lock. Either it's our responsibility to
       // do the preload, or someone else beat us to it but is now done, so we
       // can go ahead freely. Check the ".done" lock to see which.
 
       lockDone = fopen((lockName+".done").c_str(), "w");
       ret = flock(fileno(lockDone), LOCK_EX | LOCK_NB);
       // If we succesfully get this lock, that means that no running processes
       // have preloaded the file for us.
       if(ret == 0){
         // Drop it again, we were just testing
         fclose(lockDone);
         lockDone = 0;
       }
       else{
         // Otherwise, someone else has already done the preload. Don't hold the
         // rest of the processes up.
         fclose(lockMain);
         lockMain = 0;
       }
     }
 
 
     FILE* f = fopen(fname.c_str(), "r");
     assert(f);
 
     struct stat st;
     fstat(fileno(f), &st);
     const long fsize = st.st_size;
 
     void* ret = mmap(addr, fsize, PROT_READ, MAP_SHARED | MAP_FIXED, fileno(f), 0);
     fclose(f);
     assert(ret != (void*)-1);
 
 
     if(touchAll){
       // Touch all the memory, one byte on each page
       char tot = 0;
       for(long i = 0; i < fsize; i += 4096){
         if((i/4096)%(fsize/409600) == 0) std::cout << (100*i)/fsize << "%" << std::endl;
         tot += ((char*)addr)[i];
       }
       // Printing forces it not to get optimized out
       std::cout << "Fully loaded: " << int(tot) << std::endl;
 
       if(!lockName.empty()){
         lockDone = fopen((lockName+".done").c_str(), "w");
         // Mark that we succesfully preloaded the file. This indication lasts
         // until the end of this process.
         flock(fileno(lockDone), LOCK_SH);
 
         // Release the main lock if we still have it
         if(lockMain) fclose(lockMain);
       }
     }
 
     return fsize;
   }

std::ostream & lem::operator<<	(	std::ostream &	os,
		const lem::LiteHit &	h
	)

Definition at line 21 of file LiteHit.cxx.

References lem::LiteHit::Cell(), lem::LiteHit::pecorr, and lem::LiteHit::Plane().

   {
     os << h.pecorr << " PECorr on plane " << h.Plane() << ", cell " << h.Cell();
     return os;
   }

std::ostream & lem::operator<<	(	std::ostream &	os,
		const lem::LibrarySummary &	ls
	)

Definition at line 32 of file LibrarySummary.cxx.

References lem::LibrarySummary::N, lem::LibrarySummary::totBkg, lem::LibrarySummary::totEnrich, and lem::LibrarySummary::totSig.

   {
     os << ls.N << " events, weighted to "
        << ls.totSig << " signal, "
        << ls.totBkg << " background and "
        << ls.totEnrich << " enriched background";
 
     return os;
   }

std::ostream & lem::operator<<	(	std::ostream &	os,
		const lem::Heads &	hs
	)

Definition at line 140 of file Heads.cxx.

References lem::Heads::fHeadSeqsLen, and lem::Heads::fNHeads.

Referenced by lem::LiteHit::operator<().

   {
     os << hs.fNHeads << " heads, with " << hs.fHeadSeqsLen << " events each.";
     return os;
   }

std::ostream & lem::operator<<	(	std::ostream &	os,
		const lem::HeadsTranspose &	hs
	)

Definition at line 147 of file Heads.cxx.

References lem::HeadsTranspose::fLibSize, and lem::HeadsTranspose::fNHeads.

   {
     os << hs.fNHeads << " heads, for a librar of " << hs.fLibSize << " events.";
     return os;
   }

std::string lem::post_query	(	const std::string &	host,
		const std::string &	query,
		const std::string &	postdata,
		int	minPort,
		int	maxPort,
		int	queryTimeout,
		int	retryTimeout,
		bool &	ok
	)

Make an HTTP POST query.

Definition at line 162 of file WebUtils.cxx.

References get_http_response_internal().

Referenced by lem::LEMClient::produce(), and ReturnResultsThreadFunc().

   {
     return get_http_response_internal(host, query, postdata,
                                       minPort, maxPort,
                                       queryTimeout, retryTimeout,
                                       ok);
   }

void* lem::ReturnResultsThreadFunc ( void * x )

Definition at line 100 of file LEMWorkerOutput_module.cc.

References runNovaSAM::args, om::cerr, allTimeWatchdog::endl, lem::LEMWebSettings::host, lem::LEMWebSettings::maxPort, lem::LEMWebSettings::minPort, lem::ReturnResultsArgs::params, post_query(), lem::ReturnResultsArgs::postdata, lem::LEMWebSettings::queryTimeout, lem::LEMWebSettings::retryTimeout, and lem::ReturnResultsArgs::web.

Referenced by lem::LEMWorkerOutput::ReturnResults().

   {
     const ReturnResultsArgs* args = (ReturnResultsArgs*)x;
 
     bool ok;
     post_query(args->web.host(),
                args->params, args->postdata,
                args->web.minPort(), args->web.maxPort(),
                args->web.queryTimeout(), args->web.retryTimeout(),
                ok);
 
     if(!ok){
       std::cerr << "No response from server after many retries. Aborting"
                 << std::endl;
       abort();
     }
 
     delete args;
 
     return 0;
   }

void lem::SetOnDemandMemoryLimit ( long limit_mb )

Definition at line 117 of file OnDemand.cxx.

Referenced by lem::LEMLibrary::reconfigure(), lem::FindMatchesTranspose::reconfigure(), and lem::LEM::reconfigure().

   {
     // Convert to bytes
     gMemoryLimit = limit_mb<<20;
   }

double lem::SimpleSurvivalProb ( const EventSummary & evt )

Definition at line 159 of file FindMatchesTranspose_module.cc.

References abs(), lem::EventSummary::ccnc, e, lem::EventSummary::IsSig(), lem::EventSummary::origPdg, lem::EventSummary::pdg, std::sin(), util::sqr(), and lem::EventSummary::trueEVis.

   {
     if(evt.IsSig()){
       // Scale up to keep as much signal as possible in sample, weighted
       // correctly internally. The missing prefactor is (ss2th13*ssth23)
       return util::sqr(sin(1.267*2.35e-3*810/evt.trueEVis));
     }
     if(evt.ccnc == 0 && abs(evt.pdg) == 14 && abs(evt.origPdg) == 14){
       // CC mu
       return 1-util::sqr(sin(1.267*2.35e-3*810/evt.trueEVis));
     }
     if(evt.ccnc == 0 && abs(evt.pdg) == 12 && abs(evt.origPdg) == 12){
       // Beam nue
       //      return 1-.1*util::sqr(sin(1.267*2.35e-3*810/evt.trueEVis));
       // Want library matching to ignore nues
       return 0;
     }
     if(evt.ccnc == 1 && evt.pdg == evt.origPdg){
       // NC
       return 1;
     }
     if(evt.ccnc == 1 && evt.pdg != evt.origPdg){
       // Swapped NC. We play clever reweighting games so this is OK
       return 1;
     }
     if(evt.ccnc == 0 && abs(evt.pdg) == 14 && abs(evt.origPdg) == 12){
       // CC mu appearance, take ss2th13=0.1, maximal th23
       return .1*.5*util::sqr(sin(1.267*2.35e-3*810/evt.trueEVis));
     }
     // Don't know what this is
     return 0;
   }

void lem::VertexToPlaneAndCell	(	const TVector3	vtx,
		const rb::Cluster &	slice,
		int &	plane,
		int &	cell,
		int &	cellOtherView,
		bool	reverse
	)

Definition at line 60 of file Util.cxx.

References geo::PlaneGeo::Cell(), DefaultVertex(), FindCellRobust(), geom(), geo::CellGeo::GetCenter(), geo::GeometryBase::IdToCell(), geo::kPLANE_NOT_FOUND, geo::PlaneGeo::Ncells(), ncells, geo::GeometryBase::NextPlaneOtherView(), and geo::GeometryBase::Plane().

Referenced by lem::LEMSummarizer::produce(), lem::AlignLibToVtx::produce(), and lem::MakeLibrary::produce().

   {
     // Picked up by failure cases below
     DefaultVertex(slice, plane, cell, cellOtherView, reverse);
 
     art::ServiceHandle<geo::Geometry> geom;
 
     const geo::CellUniqueId id = FindCellRobust(vtx.X(),
                                                 vtx.Y(),
                                                 vtx.Z());
 
     // Take defaults
     if(!id) return;
 
     geom->IdToCell(id, &plane, &cell);
 
     const int nextPlane = geom->NextPlaneOtherView(plane);
     // Take default in other view
 
     if(nextPlane == geo::kPLANE_NOT_FOUND) return;
 
     double xyz[3];
     int ncells = geom->Plane(nextPlane)->Ncells();
     if(cellOtherView >= ncells)
       cellOtherView = ncells-1;
 
     geom->Plane(nextPlane)->Cell(cellOtherView)->GetCenter(xyz);
 
     // z position of plane in other view
     double zOther = xyz[2];
 
     // Sometimes this goes wrong and we have to iterate a bit
     while(true){
       const geo::CellUniqueId idOther = FindCellRobust(vtx.X(),
                                                        vtx.Y(),
                                                        zOther);
       // Take default in other view
       if(!idOther) return;
 
       // The cell we found must be in the plane we wanted
       int checkPlane;
       geom->IdToCell(idOther, &checkPlane, &cellOtherView);
       // If so, we're good and return the answer
       if(checkPlane == nextPlane) return;
       // Otherwise, walk 1cm at a time in the right direction
       if(checkPlane < nextPlane) zOther += 1;
       if(checkPlane > nextPlane) zOther -= 1;
     }
   }

size_t lem::writeMemoryCallback	(	void *	contents,
		size_t	size,
		size_t	nmemb,
		void *	userp
	)

Definition at line 17 of file WebUtils.cxx.

References ana::assert(), and runNovaSAM::str.

Referenced by get_http_response_internal().

   {
     assert(userp);
     std::string* str = (std::string*)userp;
 
     const size_t recvdSize = size*nmemb;
     const size_t prevSize = str->size();
     str->resize(prevSize+recvdSize);
     memcpy((char*)str->c_str()+prevSize, contents, recvdSize);
     return recvdSize;
   }