/w/halld-scifs17exp/halld2/home/sdobbs/Software/jana/jana_0.8.2/Linux_CentOS7.7-x86

Bug Summary

File:	alld2/home/sdobbs/Software/jana/jana_0.8.2/Linux_CentOS7.7-x86_64-gcc4.8.5/include/JANA/JEventLoop.h
Warning:	line 398, column 7 Null pointer passed to 2nd parameter expecting 'nonnull'

Annotated Source Code

Press '?' to see keyboard shortcuts

Show analyzer invocation

clang -cc1 -cc1 -triple x86_64-unknown-linux-gnu -analyze -disable-free -main-file-name DAnalysisUtilities.cc -analyzer-store=region -analyzer-opt-analyze-nested-blocks -analyzer-checker=core -analyzer-checker=apiModeling -analyzer-checker=unix -analyzer-checker=deadcode -analyzer-checker=cplusplus -analyzer-checker=security.insecureAPI.UncheckedReturn -analyzer-checker=security.insecureAPI.getpw -analyzer-checker=security.insecureAPI.gets -analyzer-checker=security.insecureAPI.mktemp -analyzer-checker=security.insecureAPI.mkstemp -analyzer-checker=security.insecureAPI.vfork -analyzer-checker=nullability.NullPassedToNonnull -analyzer-checker=nullability.NullReturnedFromNonnull -analyzer-output plist -w -setup-static-analyzer -mrelocation-model pic -pic-level 2 -fhalf-no-semantic-interposition -mframe-pointer=none -fmath-errno -fno-rounding-math -mconstructor-aliases -munwind-tables -target-cpu x86-64 -tune-cpu generic -fno-split-dwarf-inlining -debugger-tuning=gdb -resource-dir /w/halld-scifs17exp/home/sdobbs/clang/llvm-project/install/lib/clang/12.0.0 -D HAVE_CCDB -D HAVE_RCDB -D HAVE_EVIO -D HAVE_TMVA=1 -D RCDB_MYSQL=1 -D RCDB_SQLITE=1 -D SQLITE_USE_LEGACY_STRUCT=ON -I .Linux_CentOS7.7-x86_64-gcc4.8.5/libraries/ANALYSIS -I libraries/ANALYSIS -I . -I libraries -I libraries/include -I /w/halld-scifs17exp/home/sdobbs/clang/halld_recon/Linux_CentOS7.7-x86_64-gcc4.8.5/include -I external/xstream/include -I /usr/include/tirpc -I /group/halld/Software/builds/Linux_CentOS7.7-x86_64-gcc4.8.5/root/root-6.08.06/include -I /w/halld-scifs17exp/halld2/home/sdobbs/Software/jana/jana_0.8.2/Linux_CentOS7.7-x86_64-gcc4.8.5/include -I /group/halld/Software/builds/Linux_CentOS7.7-x86_64-gcc4.8.5/ccdb/ccdb_1.06.06/include -I /group/halld/Software/builds/Linux_CentOS7.7-x86_64-gcc4.8.5/rcdb/rcdb_0.06.00/cpp/include -I /usr/include/mysql -I /group/halld/Software/builds/Linux_CentOS7.7-x86_64-gcc4.8.5/sqlitecpp/SQLiteCpp-2.2.0^bs130/include -I /group/halld/Software/builds/Linux_CentOS7.7-x86_64-gcc4.8.5/sqlite/sqlite-3.13.0^bs130/include -I /group/halld/Software/builds/Linux_CentOS7.7-x86_64-gcc4.8.5/hdds/hdds-4.9.0/Linux_CentOS7.7-x86_64-gcc4.8.5/src -I /group/halld/Software/builds/Linux_CentOS7.7-x86_64-gcc4.8.5/xerces-c/xerces-c-3.1.4/include -I /group/halld/Software/builds/Linux_CentOS7.7-x86_64-gcc4.8.5/evio/evio-4.4.6/Linux-x86_64/include -internal-isystem /usr/lib/gcc/x86_64-redhat-linux/4.8.5/../../../../include/c++/4.8.5 -internal-isystem /usr/lib/gcc/x86_64-redhat-linux/4.8.5/../../../../include/c++/4.8.5/x86_64-redhat-linux -internal-isystem /usr/lib/gcc/x86_64-redhat-linux/4.8.5/../../../../include/c++/4.8.5/backward -internal-isystem /usr/local/include -internal-isystem /w/halld-scifs17exp/home/sdobbs/clang/llvm-project/install/lib/clang/12.0.0/include -internal-externc-isystem /include -internal-externc-isystem /usr/include -O2 -std=c++11 -fdeprecated-macro -fdebug-compilation-dir /home/sdobbs/work/clang/halld_recon/src -ferror-limit 19 -fgnuc-version=4.2.1 -fcxx-exceptions -fexceptions -vectorize-loops -vectorize-slp -analyzer-output=html -faddrsig -o /tmp/scan-build-2021-01-21-110224-160369-1 -x c++ libraries/ANALYSIS/DAnalysisUtilities.cc

libraries/ANALYSIS/DAnalysisUtilities.cc

→

1#ifdef VTRACE
2#include "vt_user.h"
3#endif

5#include "DAnalysisUtilities.h"
6#include "ANALYSIS/DParticleComboCreator.h"

8DAnalysisUtilities::DAnalysisUtilities(JEventLoop* locEventLoop)
9{
DEBUG_LEVEL=0;
gPARMS->SetDefaultParameter("DAnalysisUtilities:DEBUG_LEVEL",DEBUG_LEVEL);

locEventLoop->GetSingle(dPIDAlgorithm);

dTargetZCenter = 65.0;
// Get Target parameters from XML
DApplication *locApplication = dynamic_cast<DApplication*> (locEventLoop->GetJApplication());
DGeometry *locGeometry = locApplication ? locApplication->GetDGeometry(locEventLoop->GetJEvent().GetRunNumber()) : NULL__null;
locGeometry->GetTargetZ(dTargetZCenter);

//Get magnetic field map
dMagneticFieldMap = locApplication->GetBfield(locEventLoop->GetJEvent().GetRunNumber());
dIsNoFieldFlag = (dynamic_cast<const DMagneticFieldMapNoField*>(dMagneticFieldMap) != NULL__null);

//For "Unused" tracks/showers
//BEWARE: IF THIS IS CHANGED, CHANGE IN THE BLUEPRINT FACTORY AND THE EVENT WRITER ALSO!!
dTrackSelectionTag = "PreSelect";
dShowerSelectionTag = "PreSelect";
gPARMS->SetDefaultParameter("COMBO:TRACK_SELECT_TAG", dTrackSelectionTag);
gPARMS->SetDefaultParameter("COMBO:SHOWER_SELECT_TAG", dShowerSelectionTag);
31}

33bool DAnalysisUtilities::Check_IsBDTSignalEvent(JEventLoop* locEventLoop, const DReaction* locReaction, bool locExclusiveMatchFlag, bool locIncludeDecayingToReactionFlag) const
34{
35#ifdef VTRACE
VT_TRACER("DAnalysisUtilities::Check_IsBDTSignalEvent()");
37#endif

//IF DREACTION HAS A MISSING UNKNOWN PARTICLE, MUST USE locExclusiveMatchFlag = false

//if locIncludeDecayingToReactionFlag = true, will test whether the thrown reaction could decay to the DReaction
//Note that resonances, phi's, and omega's are automatically decayed
	//e.g. if DReaction or thrown is g, p -> pi+, pi-, omega, p; will instead treat it as g, p -> 2pi+, 2pi-, pi0, p (or whatever the omega decay products are)
//e.g. g, p -> pi+, pi0, K0, Lambda can decay to g, p -> 2pi+, 2pi-, pi0, p
//if locIncludeDecayingToReactionFlag = true, then it would be included as "Signal," if false, then background
//locIncludeDecayingToReactionFlag should be true UNLESS you are explicitly checking all possible reactions that could decay to your channel in your BDT
	//e.g. could kinfit to g, p -> pi+, pi0, K0, Lambda and include it as a BDT variable

if(dParticleComboCreator == nullptr) //Can't create in constructor: infinite recursion
dParticleComboCreator = new DParticleComboCreator(locEventLoop, nullptr, nullptr, nullptr);
DReaction_factory_Thrown* dThrownReactionFactory = static_cast<DReaction_factory_Thrown*>(locEventLoop->GetFactory("DReaction", "Thrown"));
      if(!dThrownReactionFactory->brun_was_called())
{
          dThrownReactionFactory->brun(locEventLoop, locEventLoop->GetJEvent().GetRunNumber());
          dThrownReactionFactory->Set_brun_called();
}

vector<const DReaction*> locThrownReactions;
locEventLoop->Get(locThrownReactions, "Thrown");
if(locThrownReactions.empty())
return false;
auto locActualThrownReaction = locThrownReactions[0];

//Replace omega & phi in DReaction with their decay products
size_t locStepIndex = 0;
vector<DReactionStep> locNewReactionSteps;
const DReaction* locCurrentReaction = locReaction;
DReaction locNewReaction("Interim"); //if needed

do
{
const DReactionStep* locReactionStep = locCurrentReaction->Get_ReactionStep(locStepIndex);
Particle_t locInitialPID = locReactionStep->Get_InitialPID();
if((locInitialPID != omega) && (locInitialPID != phiMeson))
{
	++locStepIndex;
	continue;
}

//is decaying phi or omega, replace it with its decay products
auto locDecayProducts = locReactionStep->Get_FinalPIDs();

//find the production step
for(size_t loc_i = 0; loc_i < locStepIndex; ++loc_i)
{
	const DReactionStep* locProductionStep = locCurrentReaction->Get_ReactionStep(loc_i);
	auto locPIDs = locProductionStep->Get_FinalPIDs();

	//search for the decaying particle. when found, replace it with its decay products
	bool locFoundFlag = false;
	for(size_t loc_j = 0; loc_j < locPIDs.size(); ++loc_j)
	{
		Particle_t locFinalStatePID = locPIDs[loc_j];
		if(locFinalStatePID != locInitialPID)
			continue;
		locPIDs.erase(locPIDs.begin() + loc_j);
		locPIDs.insert(locPIDs.begin(), locDecayProducts.begin(), locDecayProducts.end());
		locFoundFlag = true;
		break;
	}
	if(!locFoundFlag)
		continue;

	//make a new reaction step
	DReactionStep locNewReactionStep;
	locNewReactionStep.Set_InitialParticleID(locProductionStep->Get_InitialPID());
	locNewReactionStep.Set_TargetParticleID(locProductionStep->Get_TargetPID());
	int locMissingParticleIndex = locProductionStep->Get_MissingParticleIndex();
	for(int loc_j = 0; loc_j < int(locPIDs.size()); ++loc_j)
		locNewReactionStep.Add_FinalParticleID(locPIDs[loc_j], (loc_j == locMissingParticleIndex));
	locNewReactionSteps.push_back(locNewReactionStep);

	//make a new reaction
	DReaction locReactionBuild("Build");
	locReactionBuild.Clear_ReactionSteps();
	for(size_t loc_j = 0; loc_j < locCurrentReaction->Get_NumReactionSteps(); ++loc_j)
	{
		if(loc_j == loc_i)
			locReactionBuild.Add_ReactionStep(&locNewReactionSteps.back());
		else if(loc_j == locStepIndex)
			continue;
		else
			locReactionBuild.Add_ReactionStep(locCurrentReaction->Get_ReactionStep(loc_j));
	}
	locNewReaction = locReactionBuild;
	locCurrentReaction = &locNewReaction;
	break;
}
}
while(locStepIndex < locCurrentReaction->Get_NumReactionSteps());

//Get thrown steps
deque<pair<const DMCThrown*, deque<const DMCThrown*> > > locThrownSteps;
Get_ThrownParticleSteps(locEventLoop, locThrownSteps);

if(locThrownSteps.size() == 1) //nothing to replace: it either works or it doesn't
return Check_ThrownsMatchReaction(locEventLoop, locCurrentReaction, locExclusiveMatchFlag);

//build maps of counts of the thrown & DReaction decaying particles
map<Particle_t, size_t> locNumDecayingParticles_Thrown;
for(size_t loc_i = 0; loc_i < locThrownSteps.size(); ++loc_i)
{
if(locThrownSteps[loc_i].first == NULL__null)
	continue; //beam particle
Particle_t locPID = locThrownSteps[loc_i].first->PID();
if(locNumDecayingParticles_Thrown.find(locPID) == locNumDecayingParticles_Thrown.end())
	locNumDecayingParticles_Thrown[locPID] = 1;
else
	++locNumDecayingParticles_Thrown[locPID];
}

map<Particle_t, size_t> locNumDecayingParticles_Reaction;
for(size_t loc_i = 0; loc_i < locCurrentReaction->Get_NumReactionSteps(); ++loc_i)
{
Particle_t locPID = locCurrentReaction->Get_ReactionStep(loc_i)->Get_InitialPID();
if(locPID == Gamma)
	continue; //beam particle
if(locNumDecayingParticles_Reaction.find(locPID) == locNumDecayingParticles_Reaction.end())
	locNumDecayingParticles_Reaction[locPID] = 1;
else
	++locNumDecayingParticles_Reaction[locPID];
}

//if there is a particle type in the DReaction that is not present in the thrown, it's gonna fail no matter what: check now
map<Particle_t, size_t>::iterator locIterator = locNumDecayingParticles_Reaction.begin();
for(; locIterator != locNumDecayingParticles_Reaction.end(); ++locIterator)
{
Particle_t locPID = locIterator->first;
if(locNumDecayingParticles_Thrown.find(locPID) == locNumDecayingParticles_Thrown.end())
	return false;
}

//loop through thrown steps: replace phi's, omega's with their decay products
//also replace particles not in the DReaction with their decay products IF locIncludeDecayingToReactionFlag = true
locStepIndex = 1;
do
{
pair<const DMCThrown*, deque<const DMCThrown*> > locStepPair = locThrownSteps[locStepIndex];

//check to see if the thrown decaying particle is present in the dreaction
const DMCThrown* locThrownParent = locStepPair.first;
Particle_t locInitialPID = locThrownParent->PID();
bool locInitialPIDFoundFlag = (locNumDecayingParticles_Reaction.find(locInitialPID) != locNumDecayingParticles_Reaction.end());

//if it was not found, the only way it can be a match is if the decay products ARE in the reaction step: decay it in place
	//also: omega and phi have non-negligible width: cannot constrain their peaks in kinfit or BDT: replace with their decay products
if(((!locInitialPIDFoundFlag) && locIncludeDecayingToReactionFlag) || (locInitialPID == omega) || (locInitialPID == phiMeson))
{
	Replace_DecayingParticleWithProducts(locThrownSteps, locStepIndex);
	if(locNumDecayingParticles_Thrown[locInitialPID] == 1)
		locNumDecayingParticles_Thrown.erase(locInitialPID);
	else
		--locNumDecayingParticles_Thrown[locInitialPID];
}
else
	++locStepIndex;
}
while(locStepIndex < locThrownSteps.size());

if(!locIncludeDecayingToReactionFlag)
{
//don't try decaying thrown particles: compare as-is
DReaction* locThrownReaction = dThrownReactionFactory->Build_ThrownReaction(locEventLoop, locThrownSteps);
auto locThrownCombo = dParticleComboCreator->Build_ThrownCombo(locEventLoop, locThrownReaction, locThrownSteps);
bool locCheckResult = Check_ThrownsMatchReaction(locActualThrownReaction, locThrownCombo, locCurrentReaction, locExclusiveMatchFlag);

dParticleComboCreator->Reset();
dThrownReactionFactory->Recycle_Reaction(locThrownReaction);
return locCheckResult;
}

//at this point, #-steps are final. If exclusive, bail if they don't match
if(locExclusiveMatchFlag)
{
if(locReaction->Get_NumReactionSteps() != locThrownSteps.size())
	return false;
}

//ok, if there are still an unequal # of parents for a given PID between thrown & DReaction (i.e. both #'s are non-zero):
//it's too confusing to figure out which should be replaced by their decay products and which shouldn't
//so, try all possibilities, and see if any of them match. 

//build PIDs-to-replace vectors //easiest to manipulate a 1D vector
vector<Particle_t> locPIDVector;
vector<int> locResumeAtIndex;
for(locIterator = locNumDecayingParticles_Reaction.begin(); locIterator != locNumDecayingParticles_Reaction.end(); ++locIterator)
{
Particle_t locPID = locIterator->first;
size_t locNumThrown = locNumDecayingParticles_Thrown[locPID];
if(locNumThrown < locIterator->second)
	return false; //thrown doesn't match
if(locNumThrown == locIterator->second)
	continue; //all is well
for(size_t loc_i = 0; loc_i < locNumThrown - locIterator->second; ++loc_i)
{
	locPIDVector.push_back(locPID);
	locResumeAtIndex.push_back(0);
}
}

//if no additional replacements to make: check it
if(locPIDVector.empty())
{
DReaction* locThrownReaction = dThrownReactionFactory->Build_ThrownReaction(locEventLoop, locThrownSteps);
auto locThrownCombo = dParticleComboCreator->Build_ThrownCombo(locEventLoop, locThrownReaction, locThrownSteps);
bool locCheckResult = Check_ThrownsMatchReaction(locActualThrownReaction, locThrownCombo, locCurrentReaction, locExclusiveMatchFlag);

dParticleComboCreator->Reset();
dThrownReactionFactory->Recycle_Reaction(locThrownReaction);
return locCheckResult;
}

//loop through, making all combos of particles, (almost) just like in DParticleComboBlueprint_factory
//unlike that factory, the below may try a given combo twice: too hard to code against it though 
	//for a given PID, hard to compare current locResumeAtIndex to previous ones, since the thrown steps are continually replaced
deque<int> locDecayReplacementIndices;
int locParticleIndex = 0;
//below: contains "locThrownSteps" after each replacement
//1st deque index is replacement index, 
deque<deque<pair<const DMCThrown*, deque<const DMCThrown*> > > > locReplacementThrownSteps(1, locThrownSteps);
do
{
if(locParticleIndex == -1)
	break; //no success

deque<pair<const DMCThrown*, deque<const DMCThrown*> > > locCurrentThrownSteps = locReplacementThrownSteps.back();
if(locParticleIndex == int(locPIDVector.size()))
{
	//combo defined: try it
	DReaction* locThrownReaction = dThrownReactionFactory->Build_ThrownReaction(locEventLoop, locCurrentThrownSteps);
	auto locThrownCombo = dParticleComboCreator->Build_ThrownCombo(locEventLoop, locThrownReaction, locCurrentThrownSteps);
	bool locCheckResult = Check_ThrownsMatchReaction(locActualThrownReaction, locThrownCombo, locCurrentReaction, locExclusiveMatchFlag);

	dParticleComboCreator->Reset();
	dThrownReactionFactory->Recycle_Reaction(locThrownReaction);

	if(locCheckResult)
		return true; //it worked!
	locReplacementThrownSteps.pop_back();
	--locParticleIndex;
	continue;
}

Particle_t locPID = locPIDVector[locParticleIndex];
//find the next instance of this step & replace it
bool locFoundFlag = false;
for(size_t loc_i = locResumeAtIndex[locParticleIndex]; loc_i < locCurrentThrownSteps.size(); ++loc_i)
{
	if(locCurrentThrownSteps[loc_i].first == NULL__null)
		continue;
	Particle_t locInitialPID = locCurrentThrownSteps[loc_i].first->PID();
	if(locInitialPID != locPID)
		continue;
	locFoundFlag = true;

	Replace_DecayingParticleWithProducts(locCurrentThrownSteps, loc_i);
	locReplacementThrownSteps.push_back(locCurrentThrownSteps);
	locResumeAtIndex[locParticleIndex] = loc_i + 1;

	break;
}

if(locFoundFlag)
	++locParticleIndex;
else
{
	locResumeAtIndex[locParticleIndex] = 0;
	locReplacementThrownSteps.pop_back();
	--locParticleIndex;
}
}
while(true);

return false;
315}

317void DAnalysisUtilities::Replace_DecayingParticleWithProducts(deque<pair<const DMCThrown*, deque<const DMCThrown*> > >& locThrownSteps, size_t locStepIndex) const
318{
if(locStepIndex >= locThrownSteps.size())
return;

//find the step where this particle is produced at
int locInitialParticleDecayFromStepIndex = -1;
const DMCThrown* locThrownParent = locThrownSteps[locStepIndex].first;
if(locThrownParent == NULL__null)
return;

for(size_t loc_j = 0; loc_j < locStepIndex; ++loc_j)
{
for(size_t loc_k = 0; loc_k < locThrownSteps[loc_j].second.size(); ++loc_k)
{
	if(locThrownParent != locThrownSteps[loc_j].second[loc_k])
		continue;
	locInitialParticleDecayFromStepIndex = loc_j;
	break;
}
if(locInitialParticleDecayFromStepIndex != -1)
	break;
}

if(locInitialParticleDecayFromStepIndex == -1)
{
cout << "ERROR: SOMETHING IS WRONG WITH DAnalysisUtilities::Replace_DecayingParticleWithProducts(). ABORTING." << endl;
abort();
}

//insert the decay products where the production occured
deque<const DMCThrown*>& locProductionStepFinalParticles = locThrownSteps[locInitialParticleDecayFromStepIndex].second;
deque<const DMCThrown*>& locDecayStepFinalParticles = locThrownSteps[locStepIndex].second;
for(size_t loc_j = 0; loc_j < locProductionStepFinalParticles.size(); ++loc_j)
{
if(locProductionStepFinalParticles[loc_j] != locThrownParent)
	continue;
locProductionStepFinalParticles.erase(locProductionStepFinalParticles.begin() + loc_j);
locProductionStepFinalParticles.insert(locProductionStepFinalParticles.end(), locDecayStepFinalParticles.begin(), locDecayStepFinalParticles.end());
break;
}

locThrownSteps.erase(locThrownSteps.begin() + locStepIndex);
360}

362bool DAnalysisUtilities::Check_ThrownsMatchReaction(JEventLoop* locEventLoop, const DReaction* locReaction, bool locExclusiveMatchFlag) const
363{
//IF DREACTION HAS A MISSING UNKNOWN PARTICLE, MUST USE locExclusiveMatchFlag = false

//note, if you decay a final state particle (e.g. k+, pi+) in your input DReaction*, a match will NOT be found: the thrown reaction/combo is truncated
//if locExclusiveMatchFlag = false, then allow the input DReaction to be a subset of the thrown
if(dParticleComboCreator == nullptr) //Can't create in constructor: infinite recursion
dParticleComboCreator = new DParticleComboCreator(locEventLoop, nullptr, nullptr, nullptr);
auto locThrownCombo = dParticleComboCreator->Build_ThrownCombo(locEventLoop);

vector<const DReaction*> locReactions;
locEventLoop->Get(locReactions, "Thrown");
if(locReactions.empty())
return false;

auto locResult = Check_ThrownsMatchReaction(locReactions[0], locThrownCombo, locReaction, locExclusiveMatchFlag);
dParticleComboCreator->Reset();
return locResult;
380}

382bool DAnalysisUtilities::Check_ThrownsMatchReaction(const DReaction* locThrownReaction, const DParticleCombo* locThrownCombo, const DReaction* locReaction, bool locExclusiveMatchFlag) const
383{
384#ifdef VTRACE
VT_TRACER("DAnalysisUtilities::Check_ThrownsMatchReaction()");
386#endif

//IF DREACTION HAS A MISSING UNKNOWN PARTICLE, MUST USE locExclusiveMatchFlag = false

//note, if you decay a final state particle (e.g. k+, pi+) in your input DReaction*, a match will NOT be found: the thrown reaction/combo is truncated
//if locExclusiveMatchFlag = false, then allow the input DReaction to be a subset of the thrown

if(locThrownCombo == NULL__null)
return false;
if(locExclusiveMatchFlag)
{
if(locReaction->Get_NumReactionSteps() != locThrownCombo->Get_NumParticleComboSteps())
	return false;
}

//build map of InitialParticleDecayFromStepIndex's for input reaction: assume that if it matters, the user wanted them in the input order
map<size_t, int> locReactionInitialParticleDecayFromStepIndexMap; //first is step index (of locReaction) where particle is initial, second is where is final state particle
for(size_t loc_i = 0; loc_i < locReaction->Get_NumReactionSteps(); ++loc_i)
{
const DReactionStep* locReactionStep = locReaction->Get_ReactionStep(loc_i);
if(loc_i == 0)
{
	if(locReactionStep->Get_InitialPID() == Gamma)
		locReactionInitialParticleDecayFromStepIndexMap[0] = -1;
}
//loop over final state particles, and if decaying, find the step they are a parent in
for(size_t loc_j = 0; loc_j < locReactionStep->Get_NumFinalPIDs(); ++loc_j)
{
	Particle_t locFinalStatePID = locReactionStep->Get_FinalPID(loc_j);
	//see if this pid is a parent in a future step
	for(size_t loc_k = loc_i; loc_k < locReaction->Get_NumReactionSteps(); ++loc_k)
	{
		if(locReaction->Get_ReactionStep(loc_k)->Get_InitialPID() != locFinalStatePID)
			continue;
		if(locReactionInitialParticleDecayFromStepIndexMap.find(loc_k) != locReactionInitialParticleDecayFromStepIndexMap.end())
			continue; //this step already accounted for
		locReactionInitialParticleDecayFromStepIndexMap[loc_k] = loc_i;
		break;
	}
}
}

//since throwns are more organized, loop through those and back-check to dreaction
set<size_t> locMatchedInputStepIndices; //step indices in input locReaction that are already matched-to
map<int, int> locStepMatching; //locReaction map to thrown step map
map<int, int> locReverseStepMatching; //thrown map to locReaction step map
for(size_t loc_i = 0; loc_i < locThrownCombo->Get_NumParticleComboSteps(); ++loc_i)
{
int locInitialParticleDecayFromStepIndex = DAnalysis::Get_InitialParticleDecayFromIndices(locThrownReaction, loc_i).first;

//find where to start the search for this step in locReaction
size_t locStartSearchIndex = 0; //locReaction could be a subset of the total; start at the beginning unless ...:
if(locReverseStepMatching.find(locInitialParticleDecayFromStepIndex) != locReverseStepMatching.end())
	locStartSearchIndex = locReverseStepMatching[locInitialParticleDecayFromStepIndex] + 1; //parent step was matched, start search for this one after it

//loop through locReaction and try to find this thrown step
bool locMatchFoundFlag = false;
int locPossibleMatchIndex = -1;
for(size_t loc_j = locStartSearchIndex; loc_j < locReaction->Get_NumReactionSteps(); ++loc_j)
{
	const DReactionStep* locReactionStep = locReaction->Get_ReactionStep(loc_j);
	if(locMatchedInputStepIndices.find(loc_j) != locMatchedInputStepIndices.end())
		continue; //this step was already accounted for

	//when not exact match: allow user step to have a missing unknown particle
	if(!DAnalysis::Check_ChannelEquality(locThrownReaction->Get_ReactionStep(loc_i), locReactionStep, false, !locExclusiveMatchFlag))
		continue; //particles aren't the same

	//ok, now check to make sure that the parent particle in this step was produced the same way in both thrown & locReaction
	if(locReactionInitialParticleDecayFromStepIndexMap.find(loc_j) == locReactionInitialParticleDecayFromStepIndexMap.end())
	{
		//this (loc_j) parent particle's production step wasn't listed in the locReaction: locReaction is probably a subset of the total
		//a match is possible but not certain: (e.g. locReaction is pi0, eta -> pi0, pi0 and this step (loc_i) is a pi0)
		locPossibleMatchIndex = loc_j;
		continue; //keep searching in case a different one should be used instead //will use this if another isn't found
	}

	int locReactionInitialParticleDecayFromStepIndex = locReactionInitialParticleDecayFromStepIndexMap[loc_j];
	if(locReactionInitialParticleDecayFromStepIndex != -1) //locReaction is not beam particle
	{
		if(locStepMatching.find(locReactionInitialParticleDecayFromStepIndex) == locStepMatching.end())
			continue; //the step in locReaction where this (loc_j) parent particle was produced was not mapped to the thrown steps yet: this is not the step we want
		int locReactionInitialParticleDecayFromStepIndexMappedBackToThrown = locStepMatching[locReactionInitialParticleDecayFromStepIndex];
		if(locInitialParticleDecayFromStepIndex != locReactionInitialParticleDecayFromStepIndexMappedBackToThrown)
			continue; //the decaying parent particle in this step (loc_j) comes from a different step in thrown (locInitialParticleDecayFromStepIndex)/reaction: continue
	}
	else if(locInitialParticleDecayFromStepIndex != -1)
		continue; //reaction is beam but thrown is not beam

	//finally, a match! register it
	locMatchedInputStepIndices.insert(loc_j);
	locStepMatching[loc_j] = loc_i;
	locReverseStepMatching[loc_i] = loc_j;
	locMatchFoundFlag = true;
	break;
}

if((!locMatchFoundFlag) && (locPossibleMatchIndex != -1))
{
	//need to use the possible match
	locMatchedInputStepIndices.insert(locPossibleMatchIndex);
	locStepMatching[locPossibleMatchIndex] = loc_i;
	locReverseStepMatching[loc_i] = locPossibleMatchIndex;
	locMatchFoundFlag = true;
}
if(locExclusiveMatchFlag && (!locMatchFoundFlag))
	return false; //needed an exact match and it wasn't found: bail
}
if(locExclusiveMatchFlag)
return true;

//locReaction could be a subset of thrown: check if all locReaction steps found
for(size_t loc_i = 0; loc_i < locReaction->Get_NumReactionSteps(); ++loc_i)
{
if(locMatchedInputStepIndices.find(loc_i) == locMatchedInputStepIndices.end())
	return false; //one of the input steps wasn't matched: abort!
}

return true;
505}

507void DAnalysisUtilities::Get_UnusedChargedTracks(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo, vector<const DChargedTrack*>& locUnusedChargedTracks) const
508{
locUnusedChargedTracks.clear();
locEventLoop->Get(locUnusedChargedTracks, "Combo");
std::sort(locUnusedChargedTracks.begin(), locUnusedChargedTracks.end());

auto locChargedSourceObjects = locParticleCombo->Get_FinalParticle_SourceObjects(d_Charged);
for(size_t loc_i = 0; loc_i < locChargedSourceObjects.size(); ++loc_i)
{
for(auto locIterator = locUnusedChargedTracks.begin(); locIterator != locUnusedChargedTracks.end();)
{
	if(locChargedSourceObjects[loc_i] != *locIterator)
		++locIterator; //not-used (yet)
	else
		locIterator = locUnusedChargedTracks.erase(locIterator); //used
}
}
524}

526void DAnalysisUtilities::Get_UnusedTimeBasedTracks(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo, vector<const DTrackTimeBased*>& locUnusedTimeBasedTracks) const
527{
locUnusedTimeBasedTracks.clear();
locEventLoop->Get(locUnusedTimeBasedTracks);

vector<const DTrackTimeBased*> locComboTimeBasedTracks;
locEventLoop->Get(locComboTimeBasedTracks, "Combo");
locUnusedTimeBasedTracks.insert(locUnusedTimeBasedTracks.end(), locComboTimeBasedTracks.begin(), locComboTimeBasedTracks.end());

auto locChargedSourceObjects = locParticleCombo->Get_FinalParticle_SourceObjects(d_Charged);
for(size_t loc_i = 0; loc_i < locChargedSourceObjects.size(); ++loc_i)
{
//only need the candidate id: same for all hypotheses for a given track
auto locChargedTrack = static_cast<const DChargedTrack*>(locChargedSourceObjects[loc_i]);
for(auto locIterator = locUnusedTimeBasedTracks.begin(); locIterator != locUnusedTimeBasedTracks.end();)
{
	if(locChargedTrack->candidateid != (*locIterator)->candidateid)
		++locIterator; //not-used (yet)
	else
		locIterator = locUnusedTimeBasedTracks.erase(locIterator); //used
}
}
548}

550void DAnalysisUtilities::Get_UnusedWireBasedTracks(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo, vector<const DTrackWireBased*>& locUnusedWireBasedTracks) const
551{
locUnusedWireBasedTracks.clear();
locEventLoop->Get(locUnusedWireBasedTracks);

auto locChargedSourceObjects = locParticleCombo->Get_FinalParticle_SourceObjects(d_Charged);
for(size_t loc_i = 0; loc_i < locChargedSourceObjects.size(); ++loc_i)
{
//only need the candidate id: same for all hypotheses for a given track
auto locChargedTrack = static_cast<const DChargedTrack*>(locChargedSourceObjects[loc_i]);
for(auto locIterator = locUnusedWireBasedTracks.begin(); locIterator != locUnusedWireBasedTracks.end();)
{
	if(locChargedTrack->candidateid != (*locIterator)->candidateid)
		++locIterator; //not-used (yet)
	else
		locIterator = locUnusedWireBasedTracks.erase(locIterator); //used
}
}
568}

570void DAnalysisUtilities::Get_UnusedTrackCandidates(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo, vector<const DTrackCandidate*>& locUnusedTrackCandidates) const
571{
locUnusedTrackCandidates.clear();
locEventLoop->Get(locUnusedTrackCandidates);
1
Calling 'JEventLoop::Get'→

auto locChargedSourceObjects = locParticleCombo->Get_FinalParticle_SourceObjects(d_Charged);
set<unsigned int> locUsedCandidateIndices;
for(size_t loc_i = 0; loc_i < locChargedSourceObjects.size(); ++loc_i)
{
//only need the candidate id: same for all hypotheses for a given track
auto locChargedTrack = static_cast<const DChargedTrack*>(locChargedSourceObjects[loc_i]);
locUsedCandidateIndices.insert(locChargedTrack->candidateid - 1); //id = index + 1
}

for(int loc_i = locUnusedTrackCandidates.size() - 1; loc_i >= 0; --loc_i)
{
if(locUsedCandidateIndices.find(loc_i) != locUsedCandidateIndices.end())
	locUnusedTrackCandidates.erase(locUnusedTrackCandidates.begin() + loc_i);
}
589}

591void DAnalysisUtilities::Get_UnusedNeutralShowers(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo, vector<const DNeutralShower*>& locUnusedNeutralShowers) const
592{
locUnusedNeutralShowers.clear();
locEventLoop->Get(locUnusedNeutralShowers, dShowerSelectionTag.c_str());

auto locNeutralSourceObjects = locParticleCombo->Get_FinalParticle_SourceObjects(d_Neutral);
for(size_t loc_i = 0; loc_i < locNeutralSourceObjects.size(); ++loc_i)
{
auto locNeutralShower = static_cast<const DNeutralShower*>(locNeutralSourceObjects[loc_i]);
for(auto locIterator = locUnusedNeutralShowers.begin(); locIterator != locUnusedNeutralShowers.end();)
{
	if(locNeutralShower != *locIterator)
		++locIterator;
	else
		locIterator = locUnusedNeutralShowers.erase(locIterator);
}
}
608}

610void DAnalysisUtilities::Get_UnusedNeutralParticles(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo, vector<const DNeutralParticle*>& locUnusedNeutralParticles) const
611{
locUnusedNeutralParticles.clear();
locEventLoop->Get(locUnusedNeutralParticles, dShowerSelectionTag.c_str());

auto locNeutralSourceObjects = locParticleCombo->Get_FinalParticle_SourceObjects(d_Neutral);
for(size_t loc_i = 0; loc_i < locNeutralSourceObjects.size(); ++loc_i)
{
auto locNeutralShower = static_cast<const DNeutralShower*>(locNeutralSourceObjects[loc_i]);
for(auto locIterator = locUnusedNeutralParticles.begin(); locIterator != locUnusedNeutralParticles.end();)
{
	if(locNeutralShower != (*locIterator)->dNeutralShower)
		++locIterator;
	else
		locIterator = locUnusedNeutralParticles.erase(locIterator);
}
}
627}

629void DAnalysisUtilities::Get_ThrownParticleSteps(JEventLoop* locEventLoop, deque<pair<const DMCThrown*, deque<const DMCThrown*> > >& locThrownSteps) const
630{
vector<const DMCThrown*> locMCThrowns;
locEventLoop->Get(locMCThrowns);

map<size_t, const DMCThrown*> locIDMap; //size_t is the myid
for(size_t loc_i = 0; loc_i < locMCThrowns.size(); ++loc_i)
locIDMap[locMCThrowns[loc_i]->myid] = locMCThrowns[loc_i];

locThrownSteps.clear();
if(locMCThrowns.empty())
return;
locThrownSteps.push_back(pair<const DMCThrown*, deque<const DMCThrown*> >(NULL__null, deque<const DMCThrown*>() ) );

for(size_t loc_i = 0; loc_i < locMCThrowns.size(); ++loc_i)
{
if(IsResonance(locMCThrowns[loc_i]->PID()))
	continue; //don't include resonances in DReaction!!

if(locMCThrowns[loc_i]->PID() == Unknown)
	continue; //could be some weird pythia "resonance" like a diquark: just ignore them all

//initial checks of parent id
int locParentID = locMCThrowns[loc_i]->parentid;
if(locParentID == 0) //photoproduced
{
	locThrownSteps[0].second.push_back(locMCThrowns[loc_i]);
	continue;
}
if(locIDMap.find(locParentID) == locIDMap.end()) //produced from a particle that was not saved: spurious, don't save (e.g. product of BCAL shower)
	continue;

//initial checks of parent pid
Particle_t locParentPID = locIDMap[locParentID]->PID();
bool locDoneFlag = false;
while(((locParentPID == Unknown) || IsResonance(locParentPID)) && (!locDoneFlag))
{
	//intermediate particle, continue towards the source
	locParentID = locIDMap[locParentID]->parentid; //parent's parent
	if(locParentID == 0) //photoproduced
	{
		locThrownSteps[0].second.push_back(locMCThrowns[loc_i]);
		locDoneFlag = true;
	}
	else if(locIDMap.find(locParentID) == locIDMap.end()) //produced from a particle that was not saved: spurious, don't save (e.g. product of BCAL shower)
		locDoneFlag = true;
	else
		locParentPID = locIDMap[locParentID]->PID();
}

if(Is_FinalStateParticle(locParentPID) == 1)
	continue; //e.g. parent is a final state particle (e.g. this is a neutrino from a pion decay)

//check to see if the parent is already listed as a decaying particle //if so, add it to that step
bool locListedAsDecayingFlag = false;
for(size_t loc_j = 1; loc_j < locThrownSteps.size(); ++loc_j)
{
	if(locThrownSteps[loc_j].first->myid != locParentID)
		continue;
	locThrownSteps[loc_j].second.push_back(locMCThrowns[loc_i]);
	locListedAsDecayingFlag = true;
	break;
}
if(locListedAsDecayingFlag)
	continue;

//would add a new decay step, but first make sure that its parent is a decay product of a previous step
	//if the parent was not saved as a product, it may have been a decay product of a final state particle: don't save
const DMCThrown* locThrownParent = locIDMap[locParentID];
bool locFoundFlag = false;
for(size_t loc_j = 0; loc_j < locThrownSteps.size(); ++loc_j)
{
	for(size_t loc_k = 0; loc_k < locThrownSteps[loc_j].second.size(); ++loc_k)
	{
		if(locThrownSteps[loc_j].second[loc_k] != locThrownParent)
			continue;
		locFoundFlag = true;
		break;
	}
	if(locFoundFlag)
		break;
}
if(!locFoundFlag)
	continue;

//else add a new decay step and add this particle to it
locThrownSteps.push_back(pair<const DMCThrown*, deque<const DMCThrown*> >(locThrownParent, deque<const DMCThrown*>(1, locMCThrowns[loc_i]) ));
}

718/*
719cout << "THROWN STEPS: " << endl;
720for(size_t loc_i = 0; loc_i < locThrownSteps.size(); ++loc_i)
721{
cout << ((loc_i == 0) ? 0 : locThrownSteps[loc_i].first->myid) << ": ";
for(size_t loc_j = 0; loc_j < locThrownSteps[loc_i].second.size(); ++loc_j)
cout << locThrownSteps[loc_i].second[loc_j]->myid << ", ";
cout << endl;
726}
727*/
728}

730bool DAnalysisUtilities::Are_ThrownPIDsSameAsDesired(JEventLoop* locEventLoop, const deque<Particle_t>& locDesiredPIDs, Particle_t locMissingPID) const
731{
vector<const DMCThrown*> locMCThrowns;
locEventLoop->Get(locMCThrowns, "FinalState");
deque<Particle_t> locDesiredPIDs_Copy = locDesiredPIDs;

bool locMissingPIDMatchedFlag = false;
for(size_t loc_i = 0; loc_i < locMCThrowns.size(); ++loc_i)
{
Particle_t locPID = (Particle_t)(locMCThrowns[loc_i]->type);

if((!locMissingPIDMatchedFlag) && (locMissingPID == locPID))
{
	//matched missing
	locMissingPIDMatchedFlag = true;
	continue;
}

bool locPIDFoundFlag = false;
for(deque<Particle_t>::iterator locIterator = locDesiredPIDs_Copy.begin(); locIterator != locDesiredPIDs_Copy.end(); ++locIterator)
{
	if(*locIterator != locPID)
		continue;
	locDesiredPIDs_Copy.erase(locIterator);
	locPIDFoundFlag = true;
	break;
}
if(!locPIDFoundFlag)
	return false;
}

return (locDesiredPIDs_Copy.empty());
762}

764DLorentzVector DAnalysisUtilities::Calc_MissingP4(const DReaction* locReaction, const DParticleCombo* locParticleCombo, bool locUseKinFitDataFlag) const
765{
set<pair<const JObject*, unsigned int> > locSourceObjects;
return Calc_MissingP4(locReaction, locParticleCombo, 0, -1, set<size_t>(), locSourceObjects, locUseKinFitDataFlag);
768}

770DLorentzVector DAnalysisUtilities::Calc_MissingP4(const DReaction* locReaction, const DParticleCombo* locParticleCombo, set<pair<const JObject*, unsigned int> >& locSourceObjects, bool locUseKinFitDataFlag) const
771{
return Calc_MissingP4(locReaction, locParticleCombo, 0, -1, set<size_t>(), locSourceObjects, locUseKinFitDataFlag);
773}

775DLorentzVector DAnalysisUtilities::Calc_MissingP4(const DReaction* locReaction, const DParticleCombo* locParticleCombo, size_t locStepIndex, int locUpToStepIndex, set<size_t> locUpThroughIndices, bool locUseKinFitDataFlag) const
776{
set<pair<const JObject*, unsigned int> > locSourceObjects;
return Calc_MissingP4(locReaction, locParticleCombo, locStepIndex, locUpToStepIndex, locUpThroughIndices, locSourceObjects, locUseKinFitDataFlag);
779}

781DLorentzVector DAnalysisUtilities::Calc_MissingP4(const DReaction* locReaction, const DParticleCombo* locParticleCombo, size_t locStepIndex, int locUpToStepIndex, set<size_t> locUpThroughIndices, set<pair<const JObject*, unsigned int> >& locSourceObjects, bool locUseKinFitDataFlag) const
782{
//NOTE: this routine assumes that the p4 of a charged decaying particle with a detached vertex is the same at both vertices!
//assumes missing particle is not the beam particle
if(locUseKinFitDataFlag && (locParticleCombo->Get_KinFitResults() == NULL__null))
return Calc_MissingP4(locReaction, locParticleCombo, locStepIndex, locUpToStepIndex, locUpThroughIndices, locSourceObjects, false); //kinematic fit failed

DLorentzVector locMissingP4;
const DParticleComboStep* locParticleComboStep = locParticleCombo->Get_ParticleComboStep(locStepIndex);
auto locReactionStep = locReaction->Get_ReactionStep(locStepIndex);

const DKinematicData* locKinematicData = NULL__null;
if(locStepIndex == 0)
{
//initial particle
locKinematicData = locParticleComboStep->Get_InitialParticle_Measured();
locSourceObjects.insert(pair<const JObject*, unsigned int>(locKinematicData, abs(PDGtype(locKinematicData->PID())))); //want to use source objects for comparing
if(locUseKinFitDataFlag) //kinfit
	locKinematicData = locParticleComboStep->Get_InitialParticle();
locMissingP4 += locKinematicData->lorentzMomentum();
}

//target particle
Particle_t locPID = locReactionStep->Get_TargetPID();
if(locPID != Unknown)
{
double locMass = ParticleMass(locPID);
locMissingP4 += DLorentzVector(DVector3(0.0, 0.0, 0.0), locMass);
}

auto locParticles = locUseKinFitDataFlag ? locParticleComboStep->Get_FinalParticles() : locParticleComboStep->Get_FinalParticles_Measured();
for(size_t loc_j = 0; loc_j < locParticles.size(); ++loc_j)
{
if(int(loc_j) == locReactionStep->Get_MissingParticleIndex())
	continue; //exclude missing particle
if((int(locStepIndex) == locUpToStepIndex) && (locUpThroughIndices.find(loc_j) == locUpThroughIndices.end()))
	continue; //skip it: don't want to include it

int locDecayStepIndex = DAnalysis::Get_DecayStepIndex(locReaction, locStepIndex, loc_j);
if(locDecayStepIndex > 0) //decaying-particle
{
	//why plus? because the minus-signs are already applied during the call below
	locMissingP4 += Calc_MissingP4(locReaction, locParticleCombo, locDecayStepIndex, locUpToStepIndex, locUpThroughIndices, locSourceObjects, locUseKinFitDataFlag); //p4 returned is already < 0
}
else //detected
{
	Particle_t locPID = locReactionStep->Get_FinalPID(loc_j);
	auto locDetectedP4 = locParticles[loc_j]->lorentzMomentum();
	locMissingP4 -= locDetectedP4;
	locSourceObjects.insert(pair<const JObject*, unsigned int>(locParticleComboStep->Get_FinalParticle_SourceObject(loc_j), abs(PDGtype(locPID))));
}
}

return locMissingP4;
835}

837TMatrixFSym DAnalysisUtilities::Calc_MissingP3Covariance(const DReaction* locReaction, const DParticleCombo* locParticleCombo) const
838{
//uses measured data!
return Calc_MissingP3Covariance(locReaction, locParticleCombo, 0, -1, set<size_t>());
841}

843TMatrixFSym DAnalysisUtilities::Calc_MissingP3Covariance(const DReaction* locReaction, const DParticleCombo* locParticleCombo, size_t locStepIndex, int locUpToStepIndex, set<size_t> locUpThroughIndices) const
844{
//uses measured data!

//NOTE: this routine assumes that the p4 of a charged decaying particle with a detached vertex is the same at both vertices!
//assumes missing particle is not the beam particle

//missing covariance is just sum of covariance matrices of all particles used in the calculation
//because errors are uncorrelated: this doesn't work on kinfit data: just use kinfit matrix from missing particle then
const DParticleComboStep* locParticleComboStep = locParticleCombo->Get_ParticleComboStep(locStepIndex);
auto locReactionStep = locReaction->Get_ReactionStep(locStepIndex);
TMatrixFSym locMissingCovarianceMatrix(3);
locMissingCovarianceMatrix.Zero();

const DKinematicData* locKinematicData = NULL__null;
if(locStepIndex == 0)
{
//initial particle
locKinematicData = locParticleComboStep->Get_InitialParticle_Measured();
TMatrixFSym locParticleCovarianceMatrix = *(locKinematicData->errorMatrix().get());
locParticleCovarianceMatrix.ResizeTo(3, 3);
locMissingCovarianceMatrix += locParticleCovarianceMatrix;
}

auto locParticles = locParticleComboStep->Get_FinalParticles_Measured();
for(size_t loc_j = 0; loc_j < locParticles.size(); ++loc_j)
{
if(int(loc_j) == locReactionStep->Get_MissingParticleIndex())
	continue; //exclude missing particle
if((int(locStepIndex) == locUpToStepIndex) && (locUpThroughIndices.find(loc_j) == locUpThroughIndices.end()))
	continue; //skip it: don't want to include it

int locDecayStepIndex = DAnalysis::Get_DecayStepIndex(locReaction, locStepIndex, loc_j);
if(locDecayStepIndex > 0) //decaying-particle
	locMissingCovarianceMatrix += Calc_MissingP3Covariance(locReaction, locParticleCombo, locDecayStepIndex, locUpToStepIndex, locUpThroughIndices);
else //detected
{
	TMatrixFSym locParticleCovarianceMatrix = *(locParticles[loc_j]->errorMatrix().get());
	locParticleCovarianceMatrix.ResizeTo(3, 3);
	locMissingCovarianceMatrix += locParticleCovarianceMatrix;
}
}

return locMissingCovarianceMatrix;
887}

889DLorentzVector DAnalysisUtilities::Calc_FinalStateP4(const DReaction* locReaction, const DParticleCombo* locParticleCombo, size_t locStepIndex, bool locUseKinFitDataFlag) const
890{
set<pair<const JObject*, unsigned int> > locSourceObjects;
return Calc_FinalStateP4(locReaction, locParticleCombo, locStepIndex, set<size_t>(), locSourceObjects, locUseKinFitDataFlag);
893}

895DLorentzVector DAnalysisUtilities::Calc_FinalStateP4(const DReaction* locReaction, const DParticleCombo* locParticleCombo, size_t locStepIndex, set<pair<const JObject*, unsigned int> >& locSourceObjects, bool locUseKinFitDataFlag) const
896{
return Calc_FinalStateP4(locReaction, locParticleCombo, locStepIndex, set<size_t>(), locSourceObjects, locUseKinFitDataFlag);
898}

900DLorentzVector DAnalysisUtilities::Calc_FinalStateP4(const DReaction* locReaction, const DParticleCombo* locParticleCombo, size_t locStepIndex, set<size_t> locToIncludeIndices, bool locUseKinFitDataFlag) const
901{
set<pair<const JObject*, unsigned int> > locSourceObjects;
return Calc_FinalStateP4(locReaction, locParticleCombo, locStepIndex, locToIncludeIndices, locSourceObjects, locUseKinFitDataFlag);
904}

906DLorentzVector DAnalysisUtilities::Calc_FinalStateP4(const DReaction* locReaction, const DParticleCombo* locParticleCombo, size_t locStepIndex, set<size_t> locToIncludeIndices, set<pair<const JObject*, unsigned int> >& locSourceObjects, bool locUseKinFitDataFlag) const
907{
if(locUseKinFitDataFlag && (locParticleCombo->Get_KinFitResults() == NULL__null))
return Calc_FinalStateP4(locReaction, locParticleCombo, locStepIndex, locToIncludeIndices, locSourceObjects, false); //kinematic fit failed

DLorentzVector locFinalStateP4;
const DParticleComboStep* locParticleComboStep = locParticleCombo->Get_ParticleComboStep(locStepIndex);
if(locParticleComboStep == NULL__null)
return (DLorentzVector());
auto locReactionStep = locReaction->Get_ReactionStep(locStepIndex);

auto locParticles = locUseKinFitDataFlag ? locParticleComboStep->Get_FinalParticles() : locParticleComboStep->Get_FinalParticles_Measured();

//subtract rescattering target if any!!
if(locStepIndex != 0)
{
Particle_t locPID = locReactionStep->Get_TargetPID();
if(locPID != Unknown)
	locFinalStateP4 -= DLorentzVector(DVector3(0.0, 0.0, 0.0), ParticleMass(locPID));
}

bool locDoSubsetFlag = !locToIncludeIndices.empty();
for(size_t loc_i = 0; loc_i < locParticles.size(); ++loc_i)
{
if(locDoSubsetFlag && (locToIncludeIndices.find(loc_i) == locToIncludeIndices.end()))
	continue; //skip it: don't want to include it

if(locReactionStep->Get_MissingParticleIndex() == int(loc_i))
	return (DLorentzVector()); //bad!

int locDecayStepIndex = DAnalysis::Get_DecayStepIndex(locReaction, locStepIndex, loc_i);
if(locDecayStepIndex >= 0) //decaying particle
{
	//measured results, or not constrained by kinfit (either non-fixed mass or excluded from kinfit)
	if((!locUseKinFitDataFlag) || (!IsFixedMass(locReactionStep->Get_FinalPID(loc_i))))
		locFinalStateP4 += Calc_FinalStateP4(locReaction, locParticleCombo, locDecayStepIndex, set<size_t>(), locSourceObjects, locUseKinFitDataFlag);
	else //want kinfit results, and decaying particle p4 is constrained by kinfit
	{
		locFinalStateP4 += locParticles[loc_i]->lorentzMomentum();
		//still need source objects of decay products! dive down anyway, but ignore p4 result
		Calc_FinalStateP4(locReaction, locParticleCombo, locDecayStepIndex, set<size_t>(), locSourceObjects, locUseKinFitDataFlag);
	}
}
else //detected particle
{
	locFinalStateP4 += locParticles[loc_i]->lorentzMomentum();
	locSourceObjects.insert(pair<const JObject*, unsigned int>(locParticleComboStep->Get_FinalParticle_SourceObject(loc_i), abs(PDGtype(locParticles[loc_i]->PID()))));
}
}
return locFinalStateP4;
956}

958int DAnalysisUtilities::Calc_Energy_UnusedShowers(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo, double &locEnergy_UnusedShowers, int &locNumber_UnusedShowers_Quality, double &locEnergy_UnusedShowers_Quality) const
959{
DVector3 locVertex(0.0, 0.0, dTargetZCenter);
const DEventRFBunch* locEventRFBunch = locParticleCombo->Get_EventRFBunch();
double locRFTime = (locEventRFBunch != NULL__null) ? locEventRFBunch->dTime : numeric_limits<double>::quiet_NaN();

vector<const DNeutralShower*> locUnusedNeutralShowers;
Get_UnusedNeutralShowers(locEventLoop, locParticleCombo, locUnusedNeutralShowers);

int locNumber_UnusedShowers = 0;
locEnergy_UnusedShowers = 0.;
locEnergy_UnusedShowers_Quality = 0.; 
locNumber_UnusedShowers_Quality = 0;
for(size_t loc_i = 0; loc_i < locUnusedNeutralShowers.size(); ++loc_i) {
const DNeutralShower* locUnusedNeutralShower = locUnusedNeutralShowers[loc_i];

// requirements on unused showers
double locFlightTime = (locUnusedNeutralShower->dSpacetimeVertex.Vect() - locVertex).Mag()/SPEED_OF_LIGHT29.9792458;
double locDeltaT = locUnusedNeutralShower->dSpacetimeVertex.T() - locFlightTime - locRFTime;
// next line commented out to suppress warnings
// double locDetectorTheta = (locUnusedNeutralShower->dSpacetimeVertex.Vect()-locVertex).Theta()*180./TMath::Pi();
if(fabs(locDeltaT) > 4.)
	continue;		

locEnergy_UnusedShowers += locUnusedNeutralShower->dEnergy;
locNumber_UnusedShowers++;

if(locUnusedNeutralShower->dQuality > 0.5) {
	locEnergy_UnusedShowers_Quality += locUnusedNeutralShower->dEnergy;
               locNumber_UnusedShowers_Quality++;
}
}

return locNumber_UnusedShowers;
992}

994int DAnalysisUtilities::Calc_Momentum_UnusedTracks(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo, double &locSumPMag_UnusedTracks, TVector3 &locSumP3_UnusedTracks) const
995{
vector<const DChargedTrack*> locUnusedChargedTracks;
Get_UnusedChargedTracks(locEventLoop, locParticleCombo, locUnusedChargedTracks);

for(size_t loc_i = 0; loc_i < locUnusedChargedTracks.size(); ++loc_i) {
const DChargedTrack* locUnusedChargedTrack = locUnusedChargedTracks[loc_i];
const DChargedTrackHypothesis *locUnusedChargedTrackHypothesis = locUnusedChargedTrack->Get_BestTrackingFOM();

locSumPMag_UnusedTracks += locUnusedChargedTrackHypothesis->pmag();
locSumP3_UnusedTracks += locUnusedChargedTrackHypothesis->momentum();
}

return (int)locUnusedChargedTracks.size();
1008}

1010double DAnalysisUtilities::Calc_DOCAToVertex(const DKinematicData* locKinematicData, const DVector3& locVertex) const
1011{
DVector3 locPOCA;
return Calc_DOCAToVertex(locKinematicData, locVertex, locPOCA);
1014}

1016double DAnalysisUtilities::Calc_DOCAToVertex(const DKinematicData* locKinematicData, const DVector3& locVertex, DVector3& locPOCA) const
1017{
DVector3 locUnitDir = (1.0/locKinematicData->momentum().Mag())*locKinematicData->momentum();
DVector3 locPosition = locKinematicData->position();
return Calc_DOCAToVertex(locUnitDir, locPosition, locVertex, locPOCA);
1021}

1023double DAnalysisUtilities::Calc_DOCAToVertex(const DKinFitParticle* locKinFitParticle, const DVector3& locVertex) const
1024{
DVector3 locPOCA;
return Calc_DOCAToVertex(locKinFitParticle, locVertex, locPOCA);
1027}

1029double DAnalysisUtilities::Calc_DOCAToVertex(const DKinFitParticle* locKinFitParticle, const DVector3& locVertex, DVector3& locPOCA) const
1030{
DVector3 locUnitDir(locKinFitParticle->Get_Momentum().Unit().X(),locKinFitParticle->Get_Momentum().Unit().Y(),locKinFitParticle->Get_Momentum().Unit().Z());
DVector3 locPosition(locKinFitParticle->Get_Position().X(),locKinFitParticle->Get_Position().Y(),locKinFitParticle->Get_Position().Z());
return Calc_DOCAToVertex(locUnitDir, locPosition, locVertex, locPOCA);
1034}

1036double DAnalysisUtilities::Calc_DOCAToVertex(const DVector3& locUnitDir, const DVector3& locPosition, const DVector3& locVertex) const
1037{
DVector3 locPOCA1, locPOCA2;
return Calc_DOCA(locUnitDir, locUnitDir, locPosition, locVertex, locPOCA1, locPOCA2);
1040}

1042double DAnalysisUtilities::Calc_DOCAToVertex(const DVector3& locUnitDir, const DVector3& locPosition, const DVector3& locVertex, DVector3& locPOCA) const
1043{
DVector3 locPOCA2;
return Calc_DOCA(locUnitDir, locUnitDir, locPosition, locVertex, locPOCA, locPOCA2);
1046}

1048double DAnalysisUtilities::Calc_DOCAVertex(const DKinFitParticle* locKinFitParticle1, const DKinFitParticle* locKinFitParticle2, DVector3& locDOCAVertex) const
1049{
DVector3 locPOCA1, locPOCA2;
double locDOCA;
// Try to use helical approximation if B!=0
if (dIsNoFieldFlag==false && Calc_DOCA(locKinFitParticle1,locKinFitParticle2,
			 locPOCA1,locPOCA2,locDOCA)==NOERROR){
  locDOCAVertex=0.5*(locPOCA1+locPOCA2);
  return locDOCA;
}

// Use straight line approximation
DVector3 locUnitDir1(locKinFitParticle1->Get_Momentum().Unit().X(),locKinFitParticle1->Get_Momentum().Unit().Y(),locKinFitParticle1->Get_Momentum().Unit().Z());
DVector3 locUnitDir2(locKinFitParticle2->Get_Momentum().Unit().X(),locKinFitParticle2->Get_Momentum().Unit().Y(),locKinFitParticle2->Get_Momentum().Unit().Z());
DVector3 locVertex1(locKinFitParticle1->Get_Position().X(),locKinFitParticle1->Get_Position().Y(),locKinFitParticle1->Get_Position().Z());
DVector3 locVertex2(locKinFitParticle2->Get_Position().X(),locKinFitParticle2->Get_Position().Y(),locKinFitParticle2->Get_Position().Z());
return Calc_DOCAVertex(locUnitDir1, locUnitDir2, locVertex1, locVertex2, locDOCAVertex);
1065}

1067double DAnalysisUtilities::Calc_DOCAVertex(const DKinematicData* locKinematicData1, const DKinematicData* locKinematicData2, DVector3& locDOCAVertex) const
1068{
DVector3 locPOCA1, locPOCA2;
double locDOCA;
// Try to use helical approximation if B!=0
if (dIsNoFieldFlag==false && Calc_DOCA(locKinematicData1,locKinematicData2,
			 locPOCA1,locPOCA2,locDOCA)==NOERROR){
  locDOCAVertex=0.5*(locPOCA1+locPOCA2);
  return locDOCA;
}

// Use straight line approximation
DVector3 locUnitDir1 = (1.0/locKinematicData1->momentum().Mag())*locKinematicData1->momentum();
DVector3 locUnitDir2 = (1.0/locKinematicData2->momentum().Mag())*locKinematicData2->momentum();
DVector3 locVertex1 = locKinematicData1->position();
DVector3 locVertex2 = locKinematicData2->position();
return Calc_DOCAVertex(locUnitDir1, locUnitDir2, locVertex1, locVertex2, locDOCAVertex);
1084}

1086double DAnalysisUtilities::Calc_DOCAVertex(const DVector3 &locUnitDir1, const DVector3 &locUnitDir2, const DVector3 &locVertex1, const DVector3 &locVertex2, DVector3& locDOCAVertex) const
1087{
DVector3 locPOCA1, locPOCA2;
double locDOCA = Calc_DOCA(locUnitDir1, locUnitDir2, locVertex1, locVertex2, locPOCA1, locPOCA2);
locDOCAVertex = 0.5*(locPOCA1 + locPOCA2);
return locDOCA;
1092}

1094double DAnalysisUtilities::Calc_DOCA(const DKinFitParticle* locKinFitParticle1, const DKinFitParticle* locKinFitParticle2) const
1095{
DVector3 locPOCA1, locPOCA2;
return Calc_DOCA(locKinFitParticle1, locKinFitParticle2, locPOCA1, locPOCA2);
1098}

1100double DAnalysisUtilities::Calc_DOCA(const DKinematicData* locKinematicData1, const DKinematicData* locKinematicData2) const
1101{
DVector3 locPOCA1, locPOCA2;
return Calc_DOCA(locKinematicData1, locKinematicData2, locPOCA1, locPOCA2);
1104}

1106double DAnalysisUtilities::Calc_DOCA(const DKinFitParticle* locKinFitParticle1, const DKinFitParticle* locKinFitParticle2, DVector3 &locPOCA1, DVector3 &locPOCA2) const
1107{
DVector3 locUnitDir1(locKinFitParticle1->Get_Momentum().Unit().X(),locKinFitParticle1->Get_Momentum().Unit().Y(),locKinFitParticle1->Get_Momentum().Unit().Z());
DVector3 locUnitDir2(locKinFitParticle2->Get_Momentum().Unit().X(),locKinFitParticle2->Get_Momentum().Unit().Y(),locKinFitParticle2->Get_Momentum().Unit().Z());
DVector3 locVertex1(locKinFitParticle1->Get_Position().X(),locKinFitParticle1->Get_Position().Y(),locKinFitParticle1->Get_Position().Z());
DVector3 locVertex2(locKinFitParticle2->Get_Position().X(),locKinFitParticle2->Get_Position().Y(),locKinFitParticle2->Get_Position().Z());
return Calc_DOCA(locUnitDir1, locUnitDir2, locVertex1, locVertex2, locPOCA1, locPOCA2);
1113}

1115double DAnalysisUtilities::Calc_DOCA(const DKinematicData* locKinematicData1, const DKinematicData* locKinematicData2, DVector3 &locPOCA1, DVector3 &locPOCA2) const
1116{
DVector3 locUnitDir1 = (1.0/locKinematicData1->momentum().Mag())*locKinematicData1->momentum();
DVector3 locUnitDir2 = (1.0/locKinematicData2->momentum().Mag())*locKinematicData2->momentum();
DVector3 locVertex1 = locKinematicData1->position();
DVector3 locVertex2 = locKinematicData2->position();
return Calc_DOCA(locUnitDir1, locUnitDir2, locVertex1, locVertex2, locPOCA1, locPOCA2);
1122}

1124double DAnalysisUtilities::Calc_DOCA(const DVector3 &locUnitDir1, const DVector3 &locUnitDir2, const DVector3 &locVertex1, const DVector3 &locVertex2) const
1125{
DVector3 locPOCA1, locPOCA2;
return Calc_DOCA(locUnitDir1, locUnitDir2, locVertex1, locVertex2, locPOCA1, locPOCA2);
1128}

1130double DAnalysisUtilities::Calc_DOCA(const DVector3 &locUnitDir1, const DVector3 &locUnitDir2, const DVector3 &locVertex1, const DVector3 &locVertex2, DVector3 &locPOCA1, DVector3 &locPOCA2) const
1131{
//originated from code by Jorn Langheinrich
//you can use this function to find the DOCA to a fixed point by calling this function with locUnitDir1 and 2 parallel, and the fixed vertex as locVertex2
double locUnitDot = locUnitDir1.Dot(locUnitDir2);
double locDenominator = locUnitDot*locUnitDot - 1.0; // scalar product of directions
double locDistVertToInterDOCA1 = 0.0, locDistVertToInterDOCA2 = 0.0; //distance from vertex to DOCA point

if(fabs(locDenominator) < 1.0e-15) //parallel
{
locDistVertToInterDOCA1 = (locVertex2 - locVertex1).Dot(locUnitDir2)/locUnitDot; //the opposite
locDistVertToInterDOCA2 = (locVertex1 - locVertex2).Dot(locUnitDir1)/locUnitDot;
}
else
{
double locA = (locVertex1 - locVertex2).Dot(locUnitDir1);
double locB = (locVertex1 - locVertex2).Dot(locUnitDir2);
locDistVertToInterDOCA1 = (locA - locUnitDot*locB)/locDenominator;
locDistVertToInterDOCA2 = (locUnitDot*locA - locB)/locDenominator;
}

locPOCA1 = locVertex1 + locDistVertToInterDOCA1*locUnitDir1; //intersection point of DOCA line and track 1
locPOCA2 = locVertex2 + locDistVertToInterDOCA2*locUnitDir2; //intersection point of DOCA line and track 2
return (locPOCA1 - locPOCA2).Mag();
1154}

1156// Routine for steering the code that uses a small arc length approximation to 
1157// a helical trajectory to find the doca between two tracks curving in the 
1158// magnetic field.
1159jerror_t 
1160DAnalysisUtilities::Calc_DOCA(const DKinFitParticle* locKinFitParticle1, 
	      const DKinFitParticle* locKinFitParticle2,
	      DVector3 &pos1_out,DVector3 &pos2_out,
	      double &doca) const{
// Charges for the two input tracks
double q1=locKinFitParticle1->Get_Charge();
double q2=locKinFitParticle2->Get_Charge();

// position info from the two input tracks
DVector3 pos1_in=locKinFitParticle1->Get_Position();
DVector3 pos2_in=locKinFitParticle2->Get_Position();

// momentum info from the two input tracks   
DVector3 mom1_in=locKinFitParticle1->Get_Momentum();
DVector3 mom2_in=locKinFitParticle2->Get_Momentum();

return Calc_DOCA(q1,q2,pos1_in,pos2_in,mom1_in,mom2_in,pos1_out,pos2_out,
   doca);
1178}

1180// Routine for steering the code that uses a small arc length approximation to 
1181// a helical trajectory to find the doca between two tracks curving in the 
1182// magnetic field.
1183jerror_t DAnalysisUtilities::Calc_DOCA(const DKinematicData* kinematicData1,
		       const DKinematicData* kinematicData2,
		       DVector3 &pos1_out,DVector3 &pos2_out,
		       double &doca
		       ) const {
// Charges for the two input tracks
double q1=kinematicData1->charge();
double q2=kinematicData2->charge();

// position info from the two input tracks
DVector3 pos1_in=kinematicData1->position();
DVector3 pos2_in=kinematicData2->position();  

// momentum info from the two input tracks
DVector3 mom1_in=kinematicData1->momentum();
DVector3 mom2_in=kinematicData2->momentum();

return Calc_DOCA(q1,q2,pos1_in,pos2_in,mom1_in,mom2_in,pos1_out,pos2_out,
   doca);
1202}

1204// Use a small arc length approximation to a helical trajectory to find the 
1205// doca between two tracks curving in the magnetic field.
1206jerror_t DAnalysisUtilities::Calc_DOCA(double q1,double q2,
		       const DVector3 &pos1_in,
		       const DVector3 &pos2_in,
		       const DVector3 &mom1_in,
		       const DVector3 &mom2_in,
		       DVector3 &pos1_out,DVector3 &pos2_out,
		       double &doca
		       ) const {
DVector3 avg_pos=0.5*(pos1_in+pos2_in);
DVector3 diff_pos=pos1_in-pos2_in;
double B=fabs(dMagneticFieldMap->GetBz(avg_pos.x(),avg_pos.y(),avg_pos.z()));
double kap1=-0.003*B*q1/(2.*mom1_in.Perp());
double kap2=-0.003*B*q2/(2.*mom2_in.Perp());
double dx=diff_pos.x(),dy=diff_pos.y(),dz=diff_pos.z();
double tan1=tan(M_PI_21.57079632679489661923-mom1_in.Theta());
double tan1sq=tan1*tan1;
double phi1=mom1_in.Phi();
double cos1=cos(phi1),sin1=sin(phi1);
double cos1sq=cos1*cos1,sin1sq=sin1*sin1;
double tan2=tan(M_PI_21.57079632679489661923-mom2_in.Theta());
double tan2sq=tan2*tan2;
double phi2=mom2_in.Phi();
double cos2=cos(phi2),sin2=sin(phi2);
double cos2sq=cos2*cos2,sin2sq=sin2*sin2;
double dx_sq=dx*dx,dy_sq=dy*dy;
double cos2sin1_minus_cos1sin2=cos2*sin1 - cos1*sin2;

double B1=pow(2*dx*kap2*sin1sq*sin2 - 2*dx*kap1*sin1*sin2sq - sin1sq*sin2sq 
+ tan1sq + 2*dx*kap2*sin2*tan1sq - 2*dx*kap1*sin1*tan2sq 
- 2*sin1*sin2*(2*dx_sq*kap1*kap2 + tan1*tan2) 
+ sin1sq*(1 + tan2sq) + cos1sq*(-2*cos2*dy*kap2 
				+ 4*dx*kap2*sin2 + sin2sq 
				+ tan2sq) 
+ 2*cos2*(-2*dy*kap2*sin1sq + dy*kap1*sin1*sin2 
	  - dy*kap2*tan1sq + sin1*(2*dx*dy*kap1*kap2 
				   - dz*kap2*tan1 
				   + dz*kap1*tan2)) 
+ 2*cos1*(cos2sq*dy*kap1 + dy*kap2*sin1*sin2 + dy*kap1*tan2sq 
	  + sin2*(2*dx*dy*kap1*kap2 + dz*kap2*tan1
		  - dz*kap1*tan2) 
	  - cos2*(2*dy_sq*kap1*kap2 + dx*kap2*sin1 
		  + dx*kap1*sin2 + sin1*sin2 + tan1*tan2)),2)
  + 8*(cos2sin1_minus_cos1sin2)*(cos1*kap1*(cos2 - 2*dy*kap2) 
		 + kap2*(-1 + 2*dx*kap1*sin1 - tan1sq)
		 + kap1*(sin1*sin2 + tan1*tan2))*
  (2*cos2*dy_sq*kap2*sin1 + cos2*dx*sin1*sin2 - dz*tan1 
   - 2*dx*dz*kap2*sin2*tan1 + dz*sin1*sin2*tan2 + cos2*dx*tan1*tan2 
   + dy*(-2*dx*kap2*sin1*sin2 + sin1*sin2sq + 2*cos2*dz*kap2*tan1 
  + sin2*tan1*tan2 - sin1*(1 + tan2sq)) 
   + cos1*(cos2*(2*dx*dy*kap2 + dy*sin2 + dz*tan2) 
    - dx*(2*dx*kap2*sin2 + sin2sq + tan2sq)));
if (B1<0.){
  doca=9.9e9;
  return VALUE_OUT_OF_RANGE;
}

double B2=pow(2*dx*kap2*sin1sq*sin2 - 2*dx*kap1*sin1*sin2sq - sin1sq*sin2sq 
+ tan1sq + 2*dx*kap2*sin2*tan1sq - 2*dx*kap1*sin1*tan2sq 
- 2*sin1*sin2*(2*dx_sq*kap1*kap2 + tan1*tan2) 
+ sin1sq*(1 + tan2sq) + cos1sq*(-2*cos2*dy*kap2 
				+ 4*dx*kap2*sin2 + sin2sq 
				+ tan2sq)
+ 2*cos2*(-2*dy*kap2*sin1sq + dy*kap1*sin1*sin2 -dy*kap2*tan1sq
	  + sin1*(2*dx*dy*kap1*kap2 - dz*kap2*tan1 
		  + dz*kap1*tan2)) 
+ 2*cos1*(cos2sq*dy*kap1 + dy*kap2*sin1*sin2 + dy*kap1*tan2sq 
	  + sin2*(2*dx*dy*kap1*kap2 + dz*kap2*tan1 
		  - dz*kap1*tan2) 
	  - cos2*(2*dy_sq*kap1*kap2 + dx*kap2*sin1 
		  + dx*kap1*sin2 + sin1*sin2 + tan1*tan2)),2) 
  + 8*(cos2sin1_minus_cos1sin2)*(cos1*kap1*(cos2 - 2*dy*kap2)
		 + kap2*(-1 + 2*dx*kap1*sin1 - tan1sq)
		 + kap1*(sin1*sin2 + tan1*tan2))
  * (2*cos2*dy_sq*kap2*sin1 + cos2*dx*sin1*sin2 - dz*tan1 
     - 2*dx*dz*kap2*sin2*tan1 + dz*sin1*sin2*tan2 + cos2*dx*tan1*tan2 
     + dy*(-2*dx*kap2*sin1*sin2 + sin1*sin2sq + 2*cos2*dz*kap2*tan1 
    + sin2*tan1*tan2 - sin1*(1 + tan2sq)) 
     + cos1*(cos2*(2*dx*dy*kap2 + dy*sin2 + dz*tan2) 
      - dx*(2*dx*kap2*sin2 + sin2sq + tan2sq)));
if (B2<0){
  doca=9.9e9;
  return VALUE_OUT_OF_RANGE;
}

double A1=2*cos1*cos2sq*dy*kap1 - 2*cos1sq*cos2*dy*kap2 
  - 4*cos1*cos2*dy_sq*kap1*kap2 - 2*cos1*cos2*dx*kap2*sin1 
  + 4*cos2*dx*dy*kap1*kap2*sin1 + cos2sq*sin1sq - 4*cos2*dy*kap2*sin1sq 
  - 2*cos1*cos2*dx*kap1*sin2 + 4*cos1sq*dx*kap2*sin2 
  + 4*cos1*dx*dy*kap1*kap2*sin2 - 2*cos1*cos2*sin1*sin2 
  + 2*cos2*dy*kap1*sin1*sin2 + 2*cos1*dy*kap2*sin1*sin2 
  - 4*dx_sq*kap1*kap2*sin1*sin2 + 2*dx*kap2*sin1sq*sin2 + cos1sq*sin2sq 
  - 2*dx*kap1*sin1*sin2sq - 2*cos2*dz*kap2*sin1*tan1 
  + 2*cos1*dz*kap2*sin2*tan1 + cos2sq*tan1sq - 2*cos2*dy*kap2*tan1sq 
  + 2*dx*kap2*sin2*tan1sq + sin2sq*tan1sq + 2*cos2*dz*kap1*sin1*tan2 
  - 2*cos1*dz*kap1*sin2*tan2 - 2*cos1*cos2*tan1*tan2 - 2*sin1*sin2*tan1*tan2
  + cos1sq*tan2sq + 2*cos1*dy*kap1*tan2sq - 2*dx*kap1*sin1*tan2sq 
  + sin1sq*tan2sq;
double C1=-4.*(cos2sin1_minus_cos1sin2)
  *(cos1*kap1*(cos2 - 2*dy*kap2) + kap2*(-1 + 2*dx*kap1*sin1 - tan1sq)
    + kap1*(sin1*sin2 + tan1*tan2));

double A2=2*cos1*cos2sq*dy*kap1 - 2*cos1sq*cos2*dy*kap2 
  - 4*cos1*cos2*dy_sq*kap1*kap2 - 4*cos2sq*dx*kap1*sin1 
  + 2*cos1*cos2*dx*kap2*sin1 + 4*cos2*dx*dy*kap1*kap2*sin1 + cos2sq*sin1sq 
  + 2*cos1*cos2*dx*kap1*sin2 + 4*cos1*dx*dy*kap1*kap2*sin2 
  - 2*cos1*cos2*sin1*sin2 - 2*cos2*dy*kap1*sin1*sin2 
  - 2*cos1*dy*kap2*sin1*sin2 - 4*dx_sq*kap1*kap2*sin1*sin2 
  + 2*dx*kap2*sin1sq*sin2 + cos1sq*sin2sq + 4*cos1*dy*kap1*sin2sq 
  - 2*dx*kap1*sin1*sin2sq + 2*cos2*dz*kap2*sin1*tan1 
  - 2*cos1*dz*kap2*sin2*tan1 + cos2sq*tan1sq - 2*cos2*dy*kap2*tan1sq 
  + 2*dx*kap2*sin2*tan1sq + sin2sq*tan1sq - 2*cos2*dz*kap1*sin1*tan2 
  + 2*cos1*dz*kap1*sin2*tan2 - 2*cos1*cos2*tan1*tan2 - 2*sin1*sin2*tan1*tan2
  + cos1sq*tan2sq + 2*cos1*dy*kap1*tan2sq - 2*dx*kap1*sin1*tan2sq 
  + sin1sq*tan2sq;  
double C2=4.*(cos2sin1_minus_cos1sin2)
  *(kap2*(cos1*cos2 + sin1*sin2 + tan1*tan2) 
    + kap1*(-1 + 2*cos2*dy*kap2 - 2*dx*kap2*sin2 - tan2sq));

// Arc lengths to doca points
double s1_minus=(A1-sqrt(B1))/C1;
double s2_minus=(A2-sqrt(B2))/C2; 
double s1_plus=(A1+sqrt(B1))/C1;
double s2_plus=(A2+sqrt(B2))/C2;

// Position components
double x1=pos1_in.x();
double y1=pos1_in.y();
double z1=pos1_in.z();
double x2=pos2_in.x();
double y2=pos2_in.y();
double z2=pos2_in.z(); 

// Use solution corresponding to the smaller arc lengths
if (fabs(s1_minus+s2_minus)<fabs(s1_plus+s2_plus)){
  pos1_out.SetXYZ(x1+cos1*s1_minus,y1+sin1*s1_minus,z1+tan1*s1_minus);
  pos2_out.SetXYZ(x2+cos2*s2_minus,y2+sin2*s2_minus,z2+tan2*s2_minus);
}
else{
  pos1_out.SetXYZ(x1+cos1*s1_plus,y1+sin1*s1_plus,z1+tan1*s1_plus);
  pos2_out.SetXYZ(x2+cos2*s2_plus,y2+sin2*s2_plus,z2+tan2*s2_plus);
}

doca=(pos1_out-pos2_out).Mag();

if (DEBUG_LEVEL>0){
  DVector3 myvertex=0.5*(pos1_out+pos2_out);
  printf("Vertex position (x,y,z)=(%3.2f, %3.2f, %3.2f)\n",myvertex.x(),
  myvertex.y(),myvertex.z());
}

return NOERROR;
1357}



1361double DAnalysisUtilities::Calc_CrudeTime(const vector<const DKinematicData*>& locParticles, const DVector3& locCommonVertex) const
1362{
//crudely propagate the track times to the common vertex and return the average track time
DVector3 locPOCA;
DVector3 locDeltaVertex;
DVector3 locMomentum;
double locAverageTime = 0.0;
for(size_t loc_i = 0; loc_i < locParticles.size(); ++loc_i)
{
Calc_DOCAToVertex(locParticles[loc_i], locCommonVertex, locPOCA);
locDeltaVertex = locPOCA - locParticles[loc_i]->position();
locMomentum = locParticles[loc_i]->momentum();
double locTime = locParticles[loc_i]->time() + locDeltaVertex.Dot(locMomentum)*locParticles[loc_i]->energy()/(29.9792458*locMomentum.Mag2());
locAverageTime += locTime;
}
return locAverageTime/(double(locParticles.size()));
1377}

1379double DAnalysisUtilities::Calc_CrudeTime(const vector<DKinFitParticle*>& locParticles, const DVector3& locCommonVertex) const
1380{
//crudely propagate the track times to the common vertex and return the average track time
DVector3 locPOCA;
DVector3 locDeltaVertex;
double locAverageTime = 0.0;
for(size_t loc_i = 0; loc_i < locParticles.size(); ++loc_i)
{
double locTime = 0.0;
double locE = locParticles[loc_i]->Get_ShowerEnergy();
if((locParticles[loc_i]->Get_Charge() == 0) && (locE > 0.0))
{
	double locMass = locParticles[loc_i]->Get_Mass();
	double locPMag = sqrt(locE*locE - locMass*locMass);
	DVector3 locPosition = locParticles[loc_i]->Get_Position();
	DVector3 locDPosition(locPosition.X(), locPosition.Y(), locPosition.Z());
	DVector3 locDeltaVertex = locDPosition - locCommonVertex;
	locTime = locParticles[loc_i]->Get_Time() + locDeltaVertex.Mag()*locE/(29.9792458*locPMag);
}
else
{
	Calc_DOCAToVertex(locParticles[loc_i], locCommonVertex, locPOCA);
	locDeltaVertex = locPOCA - DVector3(locParticles[loc_i]->Get_Position().X(),locParticles[loc_i]->Get_Position().Y(),locParticles[loc_i]->Get_Position().Z());
	DVector3 locMomentum(locParticles[loc_i]->Get_Momentum().X(),locParticles[loc_i]->Get_Momentum().Y(),locParticles[loc_i]->Get_Momentum().Z());
	locTime = locParticles[loc_i]->Get_Time() + locDeltaVertex.Dot(locMomentum)*locParticles[loc_i]->Get_Energy()/(29.9792458*locMomentum.Mag2());
}
locAverageTime += locTime;
}
return locAverageTime/(double(locParticles.size()));
1408}

1410DVector3 DAnalysisUtilities::Calc_CrudeVertex(const vector<const DTrackTimeBased*>& locParticles) const
1411{
//assumes tracks are straight lines
//uses the midpoint of the smallest DOCA line
DVector3 locVertex(0.0, 0.0, dTargetZCenter);

if(locParticles.size() == 0)
return locVertex;
if(locParticles.size() == 1)
return locParticles[0]->position();

double locDOCA, locSmallestDOCA;
DVector3 locTempVertex;

locSmallestDOCA = 9.9E9;
for(int loc_j = 0; loc_j < (int(locParticles.size()) - 1); ++loc_j)
{
for(size_t loc_k = loc_j + 1; loc_k < locParticles.size(); ++loc_k)
{
    locDOCA = Calc_DOCAVertex(locParticles[loc_j], locParticles[loc_k], locTempVertex);

	if(locDOCA < locSmallestDOCA)
	{
		locSmallestDOCA = locDOCA;
		locVertex = locTempVertex;
	}
}
}
return locVertex;
1439}

1441DVector3 DAnalysisUtilities::Calc_CrudeVertex(const vector<const DChargedTrackHypothesis*>& locParticles) const
1442{
//assumes tracks are straight lines
//uses the midpoint of the smallest DOCA line
DVector3 locVertex(0.0, 0.0, dTargetZCenter);

if(locParticles.size() == 0)
return locVertex;
if(locParticles.size() == 1)
return locParticles[0]->position();

double locDOCA, locSmallestDOCA;
DVector3 locTempVertex;

locSmallestDOCA = 9.9E9;
for(int loc_j = 0; loc_j < (int(locParticles.size()) - 1); ++loc_j)
{
for(size_t loc_k = loc_j + 1; loc_k < locParticles.size(); ++loc_k)
{
	locDOCA = Calc_DOCAVertex(locParticles[loc_j], locParticles[loc_k], locTempVertex);
	if(locDOCA < locSmallestDOCA)
	{
		locSmallestDOCA = locDOCA;
		locVertex = locTempVertex;
	}
}
}
return locVertex;
1469}

1471DVector3 DAnalysisUtilities::Calc_CrudeVertex(const vector<const DKinematicData*>& locParticles) const
1472{
//assumes tracks are straight lines
//uses the midpoint of the smallest DOCA line
DVector3 locVertex(0.0, 0.0, dTargetZCenter);

if(locParticles.size() == 0)
return locVertex;
if(locParticles.size() == 1)
return locParticles[0]->position();

double locDOCA, locSmallestDOCA;
DVector3 locTempVertex;

locSmallestDOCA = 9.9E9;
for(int loc_j = 0; loc_j < (int(locParticles.size()) - 1); ++loc_j)
{
for(size_t loc_k = loc_j + 1; loc_k < locParticles.size(); ++loc_k)
{
	locDOCA = Calc_DOCAVertex(locParticles[loc_j], locParticles[loc_k], locTempVertex);
	if(locDOCA < locSmallestDOCA)
	{
		locSmallestDOCA = locDOCA;
		locVertex = locTempVertex;
	}
}
}
return locVertex;
1499}

1501DVector3 DAnalysisUtilities::Calc_CrudeVertex(const vector<shared_ptr<DKinFitParticle>>& locParticles) const
1502{
//assumes tracks are straight lines
//uses the midpoint of the smallest DOCA line
DVector3 locVertex(0.0, 0.0, dTargetZCenter);

if(locParticles.size() == 0)
return locVertex;

if(locParticles.size() == 1)
 return DVector3(locParticles[0]->Get_Position().X(), locParticles[0]->Get_Position().Y(), locParticles[0]->Get_Position().Z());

double locDOCA, locSmallestDOCA;
DVector3 locTempVertex;

locSmallestDOCA = 9.9E9;
for(int loc_j = 0; loc_j < (int(locParticles.size()) - 1); ++loc_j)
{
for(size_t loc_k = loc_j + 1; loc_k < locParticles.size(); ++loc_k)
{
	locDOCA = Calc_DOCAVertex(locParticles[loc_j].get(), locParticles[loc_k].get(), locTempVertex);
	if(locDOCA < locSmallestDOCA)
	{
		locSmallestDOCA = locDOCA;
		locVertex = locTempVertex;
	}
}
}
return locVertex;
1530}

1532set<set<size_t> > DAnalysisUtilities::Build_IndexCombos(const DReactionStep* locReactionStep, deque<Particle_t> locToIncludePIDs) const
1533{
//if locToIncludePIDs is empty, will return one set with all (except missing)
set<set<size_t> > locCombos;

auto locReactionStepPIDs = locReactionStep->Get_FinalPIDs();
int locMissingParticleIndex = locReactionStep->Get_MissingParticleIndex();

if(locToIncludePIDs.empty())
{
set<size_t> locCombo;
for(size_t loc_i = 0; loc_i < locReactionStepPIDs.size(); ++loc_i)
{
	if(int(loc_i) == locMissingParticleIndex)
		continue;
	locCombo.insert(loc_i);
}
locCombos.insert(locCombo);
return locCombos;
}

//deque indices corresponds to locToIncludePIDs, and each set is what could pass for it
deque<deque<size_t> > locPossibilities(locToIncludePIDs.size(), deque<size_t>());
deque<int> locResumeAtIndices(locToIncludePIDs.size(), 0);

//build possibilities: loop over reaction PIDs
for(size_t loc_i = 0; loc_i < locReactionStepPIDs.size(); ++loc_i)
{
if(int(loc_i) == locMissingParticleIndex)
	continue;
Particle_t locPID = locReactionStepPIDs[loc_i];

//register where this is a valid option: loop over to-include PIDs
for(size_t loc_j = 0; loc_j < locToIncludePIDs.size(); ++loc_j)
{
	if(locToIncludePIDs[loc_j] == locPID)
		locPossibilities[loc_j].push_back(loc_i);
}
}

//build combos
int locParticleIndex = 0;
deque<size_t> locComboDeque;
while(true)
{
if(locParticleIndex == int(locPossibilities.size())) //end of combo: save it
{
	set<size_t> locComboSet; //convert set to deque
	for(size_t loc_i = 0; loc_i < locComboDeque.size(); ++loc_i)
		locComboSet.insert(locComboDeque[loc_i]);
	locCombos.insert(locComboSet); //saved

	if(!Handle_Decursion(locParticleIndex, locComboDeque, locResumeAtIndices, locPossibilities))
		break;
	continue;
}

int& locResumeAtIndex = locResumeAtIndices[locParticleIndex];

//if two identical pids: locResumeAtIndex must always be >= the previous locResumeAtIndex (prevents duplicates) e.g. g, p -> p, pi0, pi0
//search for same pid previously in this step
Particle_t locToIncludePID = locToIncludePIDs[locParticleIndex];
for(int loc_i = locParticleIndex - 1; loc_i >= 0; --loc_i)
{
	if(locToIncludePIDs[loc_i] == locToIncludePID)
	{
		if(locResumeAtIndex < locResumeAtIndices[loc_i])
			locResumeAtIndex = locResumeAtIndices[loc_i];
		break; //dupe type in step: resume-at advanced to next
	}
}

if(locResumeAtIndex >= int(locPossibilities[locParticleIndex].size()))
{
	if(!Handle_Decursion(locParticleIndex, locComboDeque, locResumeAtIndices, locPossibilities))
		break;
	continue;
}

// pid found
locComboDeque.push_back(locPossibilities[locParticleIndex][locResumeAtIndex]);
++locResumeAtIndex;
++locParticleIndex;
}

return locCombos;
1618}

1620bool DAnalysisUtilities::Handle_Decursion(int& locParticleIndex, deque<size_t>& locComboDeque, deque<int>& locResumeAtIndices, deque<deque<size_t> >& locPossibilities) const
1621{
do
{
if(locParticleIndex < int(locResumeAtIndices.size())) //else just saved a combo
	locResumeAtIndices[locParticleIndex] = 0; //finding this particle failed: reset

--locParticleIndex; //go to previous particle
if(locParticleIndex < 0)
	return false; //end of particles: end of finding all combos

locComboDeque.pop_back(); //reset this index
}
while(locResumeAtIndices[locParticleIndex] == int(locPossibilities[locParticleIndex].size()));

return true;
1636}

1638//The POCA cannot technically be solved analytically, but we can approximate it pretty accurately
//First, propagate the track to somewhere very close to the true POCA
//Then, the equation can be solved by substituting cos(x) = 1, sin(x) = x
1641double DAnalysisUtilities::Propagate_Track(int locCharge, const DVector3& locPropagateToPoint, DLorentzVector& locMeasuredX4, DLorentzVector& locMeasuredP4, TMatrixFSym* locCovarianceMatrix) const
1642{
//ASSUMES THAT THE B-FIELD IS IN THE +Z DIRECTION!!!!!!!!!!!!!!!
double locDistance = (locMeasuredX4.Vect() - locPropagateToPoint).Mag();
if(!Get_IsBFieldNearBeamline() || (locCharge == 0) || (locDistance < dMinDistanceForStraightTrack))
{
//use simpler methods
DVector3 locPOCA;
auto locP3 = locMeasuredP4.Vect();
Calc_DOCAToVertex(locP3.Unit(), locMeasuredX4.Vect(), locPropagateToPoint, locPOCA);
locMeasuredX4.SetVect(locPOCA);
auto locDistanceVector = locPOCA - locMeasuredX4.Vect();
//negative: if you had to propagate the track forwards, that means the path length is LESS than what you thought it was
auto locDeltaPathLength = (locP3.Dot(locDistanceVector) > 0.0) ? -1.0*locDistanceVector.Mag() : locDistanceVector.Mag();
locMeasuredX4.SetT(locMeasuredX4.T() - locDeltaPathLength/(locMeasuredP4.Beta()*SPEED_OF_LIGHT29.9792458)); //v = s/t, t = s/v = s/(beta*c) = s*E/(p*c) =
return locDeltaPathLength;
}

//propagate the track to the same z as the vertex (if pz != 0)
double locTotalDeltaPathLength = 0.0;
DLorentzVector locTempPropagatedPosition = locMeasuredX4;
DLorentzVector locTempPropagatedMomentum = locMeasuredP4;
if(fabs(locMeasuredP4.Pz()) > 0.0)
{
locTotalDeltaPathLength += (locPropagateToPoint.Z() - locMeasuredX4.Z())*locMeasuredP4.P()/locMeasuredP4.Pz(); //s = (z - z0)*p/p0z
Propagate_Track(locTotalDeltaPathLength, locCharge, locTempPropagatedPosition, locTempPropagatedMomentum, NULL__null);
}

//now, step along the track until we get close to the POCA
locTotalDeltaPathLength += Calc_PathLength_Step(locCharge, locPropagateToPoint, locTempPropagatedPosition, locTempPropagatedMomentum);

//now, the path length is very close to accurate. get the rest of the way there
locTotalDeltaPathLength += Calc_PathLength_FineGrained(locCharge, locPropagateToPoint, locTempPropagatedPosition.Vect(), locTempPropagatedMomentum.Vect());

//Finally, propagate the track the given distance and return
1676//cout << "FINAL: xyz, distance = " << locMeasuredX4.X() << ", " << locMeasuredX4.Y() << ", " << locMeasuredX4.Z() << ", " << locTotalDeltaPathLength << endl;
Propagate_Track(locTotalDeltaPathLength, locCharge, locMeasuredX4, locMeasuredP4, locCovarianceMatrix);

//negative: if you had to propagate the track forwards, that means the path length is LESS than what you thought it was
return -1.0*locTotalDeltaPathLength;
1681}

1683double DAnalysisUtilities::Calc_PathLength_Step(int locCharge, const DVector3& locPropagateToPoint, DLorentzVector& locMeasuredX4, DLorentzVector& locMeasuredP4) const
1684{
//now, step slowly along the track (in path length), trying to get closer to the POCA
//for the final calculation to work, rho*s must be small: cos(rho*s) -> 1
DVector3 locBField = Get_BField(locPropagateToPoint);
double locA = -0.00299792458*(double(locCharge))*locBField.Mag();

double locPMag = locMeasuredP4.P();
double locPathLengthOneRotation = 2.0*TMath::Pi()*locPMag/locA;
double locStepDistance = locPathLengthOneRotation/12.0; //30 degrees //e.g. for 90 degree tracks

//guard against long steps in z (e.g. for 2 degree tracks)
double locStepDeltaZ = locMeasuredP4.Pz()*locStepDistance/locPMag;
double locTargetDeltaZ = 1.0;
if(fabs(locStepDeltaZ) > 5.0)
locStepDistance = locTargetDeltaZ*locPMag/locMeasuredP4.Pz();

//First, try stepping in +pvec direction
double locStepDirection = 1.0;
double locTotalDeltaPathLength = 0.0;
double locDeltaX3 = (locMeasuredX4.Vect() - locPropagateToPoint).Mag();
double locRhoS = locStepDistance*locA/locPMag;
double locRhoSCutOff = 0.1;
while(locRhoS > locRhoSCutOff) //at this point, cos(rho*s) = 0.995: close enough
{
double locDeltaPathLength = locStepDirection*locStepDistance;
locTotalDeltaPathLength += locDeltaPathLength;
1710//cout << "STEP LOOP: xyz, distance = " << locMeasuredX4.X() << ", " << locMeasuredX4.Y() << ", " << locMeasuredX4.Z() << ", " << locDeltaPathLength << endl;
Propagate_Track(locDeltaPathLength, locCharge, locMeasuredX4, locMeasuredP4, NULL__null);

double locNewDeltaX3 = (locMeasuredX4.Vect() - locPropagateToPoint).Mag();
bool locGettingCloserFlag = (locNewDeltaX3 < locDeltaX3);
locDeltaX3 = locNewDeltaX3;

if(locGettingCloserFlag)
	continue; //stepping in correct direction, can still try to get closer

//are now farther away than before: reverse direction, decrease step size
locStepDirection *= -1.0;
locStepDistance /= 2.0;
locRhoS = locStepDistance*locA/locPMag;

//if about to break, revert to older, better position
if(locRhoS < locRhoSCutOff)
{
	locTotalDeltaPathLength -= locDeltaPathLength;
	Propagate_Track(-1.0*locDeltaPathLength, locCharge, locMeasuredX4, locMeasuredP4, NULL__null);
}
}

return locTotalDeltaPathLength;
1734}

1736double DAnalysisUtilities::Calc_PathLength_FineGrained(int locCharge, const DVector3& locPropagateToPoint, DVector3 locMeasuredPosition, DVector3 locMeasuredMomentum) const
1737{
//ASSUMES B-FIELD IS IN +Z DIRECTION!!

//distance = sqrt((delta-x)^2 + (delta-y)^2 + (delta-z)^2)
//find the delta-path-length s that minimizes the distance equation
//take derivative of the distance equation, set = 0

//Mathematica yields: F*(B*s + C) + D*cos(rho*s) + E*sin(rho*s) = 0
//(where BCDEF are various, non-s-dependent terms)
//This is unsolvable analytically. However, we know that s is small, since we're almost there

//So, substitute cos(x) = 1, sin(x) = x:
//F*(B*s + C) + D + E*rho*s = 0
//Solving gives: s = - (F*C + D) / (F*B + E*rho)
//Note that B = F

DVector3 locBField = Get_BField(locMeasuredPosition);
double locA = -0.00299792458*(double(locCharge))*locBField.Mag();
double locPMag = locMeasuredMomentum.Mag();
double locRho = locA/locPMag;

DVector3 locDeltaX3 = locMeasuredPosition - locPropagateToPoint;
double locF = locMeasuredMomentum.Pz()/locPMag;
double locC = locDeltaX3.Z();
double locD = locRho*(locMeasuredMomentum.Px()*locDeltaX3.X() + locMeasuredMomentum.Py()*locDeltaX3.Y())/locA;
double locPerpPSq = locMeasuredMomentum.Px()*locMeasuredMomentum.Px() + locMeasuredMomentum.Py()*locMeasuredMomentum.Py();
double locE = locRho*(locPerpPSq/locA - locMeasuredMomentum.Py()*locDeltaX3.X() + locMeasuredMomentum.Px()*locDeltaX3.Y())/locA;
return -1.0*(locF*locC + locD)/(locF*locF + locE*locRho);
1765}

1767void DAnalysisUtilities::Propagate_Track(double locDeltaPathLength, int locCharge, DLorentzVector& locX4, DLorentzVector& locP4, TMatrixFSym* locCovarianceMatrix) const
1768{
//ASSUMES THAT THE B-FIELD IS IN THE +Z DIRECTION!!!!!!!!!!!!!!!

DVector3 locBField = Get_BField(locX4.Vect());
if(!(locBField.Mag() > 0.0))
return;
DVector3 locH = locBField.Unit();
double locA = -0.00299792458*(double(locCharge))*locBField.Mag();

double locPMag = locP4.P();
double locRhoS = locDeltaPathLength*locA/locPMag;

DLorentzVector locDeltaX4; //x - x0
locDeltaX4.SetX(locP4.Px()*sin(locRhoS)/locA - locP4.Py()*(1.0 - cos(locRhoS))/locA);
locDeltaX4.SetY(locP4.Py()*sin(locRhoS)/locA + locP4.Px()*(1.0 - cos(locRhoS))/locA);
locDeltaX4.SetZ(locP4.Pz()*locDeltaPathLength/locPMag);
locDeltaX4.SetT(locDeltaPathLength/(locP4.Beta()*SPEED_OF_LIGHT29.9792458)); //v = s/t, t = s/v = s/(beta*c) = s*E/(p*c) =
locX4 += locDeltaX4;

if(locCovarianceMatrix == NULL__null)
{
locP4.SetVect(locP4.Vect() - locDeltaX4.Vect().Cross(locA*locH));
return; //don't update error matrix
}

//transform(i, j) = di/dj, i = new, j = old //pxyz, xyz, t
TMatrixF locTransformMatrix(7, 7);
locTransformMatrix.Zero();

double locPCubed = locPMag*locPMag*locPMag;
double locSOverP3 = locDeltaPathLength/locPCubed;

double locSOverP3_PxPx = locSOverP3*locP4.Px()*locP4.Px();
double locSOverP3_PxPy = locSOverP3*locP4.Px()*locP4.Py();
double locSOverP3_PxPz = locSOverP3*locP4.Px()*locP4.Pz();

double locSOverP3_PyPy = locSOverP3*locP4.Py()*locP4.Py();
double locSOverP3_PyPz = locSOverP3*locP4.Py()*locP4.Pz();
double locSOverP3_PzPz = locSOverP3*locP4.Pz()*locP4.Pz();

//dpx
locTransformMatrix(0, 0) = (1.0 + locA*locSOverP3_PxPy)*cos(locRhoS) + locA*locSOverP3_PxPx*sin(locRhoS);
locTransformMatrix(0, 1) = (-1.0 + locA*locSOverP3_PxPy)*sin(locRhoS) + locA*locSOverP3_PyPy*cos(locRhoS);
locTransformMatrix(0, 2) = locA*locSOverP3_PyPz*cos(locRhoS) + locA*locSOverP3_PxPz*sin(locRhoS);

//dpy
locTransformMatrix(1, 0) = -1.0*locA*locSOverP3_PxPx*cos(locRhoS) + (1.0 + locA*locSOverP3_PxPy)*sin(locRhoS);
locTransformMatrix(1, 1) = (1.0 - locA*locSOverP3_PxPy)*cos(locRhoS) + locA*locSOverP3_PyPy*sin(locRhoS);
locTransformMatrix(1, 2) = locA*locSOverP3_PyPz*sin(locRhoS) - locA*locSOverP3_PxPz*cos(locRhoS);

//dpz
locTransformMatrix(2, 2) = 1.0;

//dx
locTransformMatrix(3, 0) = (1.0/locA + locSOverP3_PxPy)*sin(locRhoS) - locSOverP3_PxPx*cos(locRhoS);
locTransformMatrix(3, 1) = -1.0/locA + (1.0/locA - locSOverP3_PxPy)*cos(locRhoS) + locSOverP3_PyPy*sin(locRhoS);
locTransformMatrix(3, 2) = locSOverP3_PyPz*sin(locRhoS) - locSOverP3_PxPz*cos(locRhoS);
locTransformMatrix(3, 3) = 1.0;

//dy
locTransformMatrix(4, 0) = 1.0/locA - (1.0/locA + locSOverP3_PxPy)*cos(locRhoS) + locSOverP3_PxPx*sin(locRhoS);
locTransformMatrix(4, 1) = (1.0/locA - locSOverP3_PxPy)*sin(locRhoS) - locSOverP3_PyPy*cos(locRhoS);
locTransformMatrix(4, 2) = -1.0*locSOverP3_PyPz*cos(locRhoS) - locSOverP3_PxPz*sin(locRhoS);
locTransformMatrix(4, 4) = 1.0;

//dz
locTransformMatrix(5, 0) = -1.0*locSOverP3_PxPz;
locTransformMatrix(5, 1) = -1.0*locSOverP3_PyPz;
locTransformMatrix(5, 2) = locDeltaPathLength/locPMag - locSOverP3_PzPz;
locTransformMatrix(5, 5) = 1.0;

//dt
locTransformMatrix(6, 0) = -1.0*locP4.M2()*locSOverP3*locP4.Px()/(SPEED_OF_LIGHT29.9792458*locP4.E());
locTransformMatrix(6, 1) = -1.0*locP4.M2()*locSOverP3*locP4.Py()/(SPEED_OF_LIGHT29.9792458*locP4.E());
locTransformMatrix(6, 2) = -1.0*locP4.M2()*locSOverP3*locP4.Pz()/(SPEED_OF_LIGHT29.9792458*locP4.E());
locTransformMatrix(6, 6) = 1.0;

//transform!!
locCovarianceMatrix->Similarity(locTransformMatrix);

//update p3 //must do at the end!!
locP4.SetVect(locP4.Vect() - locDeltaX4.Vect().Cross(locA*locH));
1850}

←

/w/halld-scifs17exp/halld2/home/sdobbs/Software/jana/jana_0.8.2/Linux_CentOS7.7-x86_64-gcc4.8.5/include/JANA/JEventLoop.h

→

1// $Id: JEventLoop.h 1763 2006-05-10 14:29:25Z davidl $
2//
3//    File: JEventLoop.h
4// Created: Wed Jun  8 12:30:51 EDT 2005
5// Creator: davidl (on Darwin wire129.jlab.org 7.8.0 powerpc)
6//

8#ifndef _JEventLoop_
9#define _JEventLoop_

11#include <sys/time.h>

13#include <vector>
14#include <list>
15#include <string>
16#include <utility>
17#include <typeinfo>
18#include <string.h>
19#include <map>
20#include <utility>
21using std::vector;
22using std::list;
23using std::string;
24using std::type_info;

26#include <JANA/jerror.h>
27#include <JANA/JObject.h>
28#include <JANA/JException.h>
29#include <JANA/JEvent.h>
30#include <JANA/JThread.h>
31#include <JANA/JFactory_base.h>
32#include <JANA/JCalibration.h>
33#include <JANA/JGeometry.h>
34#include <JANA/JResourceManager.h>
35#include <JANA/JStreamLog.h>

37// The following is here just so we can use ROOT's THtml class to generate documentation.
38#include "cint.h"


41// Place everything in JANA namespace
42namespace jana{


45template<class T> class JFactory;
46class JApplication;
47class JEventProcessor;


50class JEventLoop{
public:

friend class JApplication;

enum data_source_t{
	DATA_NOT_AVAILABLE = 1,
	DATA_FROM_CACHE,
	DATA_FROM_SOURCE,
	DATA_FROM_FACTORY
};

typedef struct{
	string caller_name;
	string caller_tag;
	string callee_name;
	string callee_tag;
	double start_time;
	double end_time;
	data_source_t data_source;
}call_stack_t;

typedef struct{
	const char* factory_name;
	string tag;
	const char* filename;
	int line;
}error_call_stack_t;

                                  JEventLoop(JApplication *app); ///< Constructor
					   virtual ~JEventLoop(); ////< Destructor
		   virtual const char* className(void){return static_className();}
			static const char* static_className(void){return "JEventLoop";}

				 JApplication* GetJApplication(void) const {return app;} ///< Get pointer to the JApplication object
                                void RefreshProcessorListFromJApplication(void); ///< Re-copy the list of JEventProcessors from JApplication
                    virtual jerror_t AddFactory(JFactory_base* factory); ///< Add a factory
					  jerror_t RemoveFactory(JFactory_base* factory); ///< Remove a factory
                      JFactory_base* GetFactory(const string data_name, const char *tag="", bool allow_deftag=true); ///< Get a specific factory pointer
              vector<JFactory_base*> GetFactories(void) const {return factories;} ///< Get all factory pointers
                                void GetFactoryNames(vector<string> &factorynames); ///< Get names of all factories in name:tag format
                                void GetFactoryNames(map<string,string> &factorynames); ///< Get names of all factories in map with key=name, value=tag
                  map<string,string> GetDefaultTags(void) const {return default_tags;}
                            jerror_t ClearFactories(void); ///< Reset all factories in preparation for next event.
					  jerror_t PrintFactories(int sparsify=0); ///< Print a list of all factories.
                            jerror_t Print(const string data_name, const char *tag=""); ///< Print the data of the given type

				 JCalibration* GetJCalibration();
              template<class T> bool GetCalib(string namepath, map<string,T> &vals);
              template<class T> bool GetCalib(string namepath, vector<T> &vals);
              template<class T> bool GetCalib(string namepath, T &val);

                          JGeometry* GetJGeometry();
              template<class T> bool GetGeom(string namepath, map<string,T> &vals);
	    template<class T> bool GetGeom(string namepath, T &val);

                   JResourceManager* GetJResourceManager(void);
                              string GetResource(string namepath);
              template<class T> bool GetResource(string namepath, T vals, int event_number=0);

                                void Initialize(void); ///< Do initializations just before event processing starts
                            jerror_t Loop(void); ///< Loop over events
                            jerror_t OneEvent(uint64_t event_number); ///< Process a specific single event (if source supports it)
                            jerror_t OneEvent(void); ///< Process a single event
                         inline void Pause(void){pause = 1;} ///< Pause event processing
                         inline void Resume(void){pause = 0;} ///< Resume event processing
                         inline void Quit(void){quit = 1;} ///< Clean up and exit the event loop
                         inline bool GetQuit(void) const {return quit;}
                                void QuitProgram(void);

                               // Support for random access of events
                          bool HasRandomAccess(void);
                          void AddToEventQueue(uint64_t event_number){ next_events_to_process.push_back(event_number); }
                          void AddToEventQueue(list<uint64_t> &event_numbers) { next_events_to_process.insert(next_events_to_process.end(), event_numbers.begin(), event_numbers.end()); }
                list<uint64_t> GetEventQueue(void){ return next_events_to_process; }
                          void ClearEventQueue(void){ next_events_to_process.clear(); }

      template<class T> JFactory<T>* GetSingle(const T* &t, const char *tag="", bool exception_if_not_one=true); ///< Get pointer to first data object from (source or factory).
      template<class T> JFactory<T>* Get(vector<const T*> &t, const char *tag="", bool allow_deftag=true); ///< Get data object pointers from (source or factory)
      template<class T> JFactory<T>* GetFromFactory(vector<const T*> &t, const char *tag="", data_source_t &data_source=null_data_source, bool allow_deftag=true); ///< Get data object pointers from factory
          template<class T> jerror_t GetFromSource(vector<const T*> &t, JFactory_base *factory=NULL__null); ///< Get data object pointers from source.
                      inline JEvent& GetJEvent(void){return event;} ///< Get pointer to the current JEvent object.
                         inline void SetJEvent(JEvent *event){this->event = *event;} ///< Set the JEvent pointer.
                         inline void SetAutoFree(int auto_free){this->auto_free = auto_free;} ///< Set the Auto-Free flag on/off
                    inline pthread_t GetPThreadID(void) const {return pthread_id;} ///< Get the pthread of the thread to which this JEventLoop belongs
                              double GetInstantaneousRate(void) const {return rate_instantaneous;} ///< Get the current event processing rate
                              double GetIntegratedRate(void) const {return rate_integrated;} ///< Get the current event processing rate
                              double GetLastEventProcessingTime(void) const {return delta_time_single;}
                        unsigned int GetNevents(void) const {return Nevents;}

                         inline bool CheckEventBoundary(uint64_t event_numberA, uint64_t event_numberB);

                         inline bool GetCallStackRecordingStatus(void){ return record_call_stack; }
                         inline void DisableCallStackRecording(void){ record_call_stack = false; }
                         inline void EnableCallStackRecording(void){ record_call_stack = true; }
                     inline void CallStackStart(JEventLoop::call_stack_t &cs, const string &caller_name, const string &caller_tag, const string callee_name, const string callee_tag);
                     inline void CallStackEnd(JEventLoop::call_stack_t &cs);
         inline vector<call_stack_t> GetCallStack(void){return call_stack;} ///< Get the current factory call stack
                         inline void AddToCallStack(call_stack_t &cs){if(record_call_stack) call_stack.push_back(cs);} ///< Add specified item to call stack record but only if record_call_stack is true
                         inline void AddToErrorCallStack(error_call_stack_t &cs){error_call_stack.push_back(cs);} ///< Add layer to the factory call stack
   inline vector<error_call_stack_t> GetErrorCallStack(void){return error_call_stack;} ///< Get the current factory error call stack
                                void PrintErrorCallStack(void); ///< Print the current factory call stack

                      const JObject* FindByID(JObject::oid_t id); ///< Find a data object by its identifier.
          template<class T> const T* FindByID(JObject::oid_t id); ///< Find a data object by its type and identifier
                      JFactory_base* FindOwner(const JObject *t); ///< Find the factory that owns a data object by pointer
                      JFactory_base* FindOwner(JObject::oid_t id); ///< Find a factory that owns a data object by identifier

                                     // User defined references
              template<class T> void SetRef(T *t);    ///< Add a user reference to this JEventLoop (must be a pointer)
                template<class T> T* GetRef(void);   ///< Get a user-defined reference of a specific type
        template<class T> vector<T*> GetRefsT(void); ///< Get all user-defined refrences of a specicif type
   vector<pair<const char*, void*> > GetRefs(void){ return user_refs; }  ///< Get copy of full list of user-defined references
              template<class T> void RemoveRef(T *t); ///< Remove user reference from list

							   // Convenience methods wrapping JEvent methods of same name
                      uint64_t GetStatus(void){return event.GetStatus();}
                          bool GetStatusBit(uint32_t bit){return event.GetStatusBit(bit);}
                          bool SetStatusBit(uint32_t bit, bool val=true){return event.SetStatusBit(bit, val);}
                          bool ClearStatusBit(uint32_t bit){return event.ClearStatusBit(bit);}
                          void ClearStatus(void){event.ClearStatus();}
                          void SetStatusBitDescription(uint32_t bit, string description){event.SetStatusBitDescription(bit, description);}
                        string GetStatusBitDescription(uint32_t bit){return event.GetStatusBitDescription(bit);}
                          void GetStatusBitDescriptions(map<uint32_t, string> &status_bit_descriptions){return event.GetStatusBitDescriptions(status_bit_descriptions);}
                              string StatusWordToString(void);

private:
JEvent event;
vector<JFactory_base*> factories;
vector<JEventProcessor*> processors;
vector<error_call_stack_t> error_call_stack;
vector<call_stack_t> call_stack;
JApplication *app;
JThread *jthread;
bool initialized;
bool print_parameters_called;
int pause;
int quit;
int auto_free;
pthread_t pthread_id;
map<string, string> default_tags;
vector<pair<string,string> > auto_activated_factories;
bool record_call_stack;
string caller_name;
string caller_tag;
vector<uint64_t> event_boundaries;
int32_t event_boundaries_run; ///< Run number boundaries were retrieved from (possbily 0)
list<uint64_t> next_events_to_process;

uint64_t Nevents;			      ///< Total events processed (this thread)
uint64_t Nevents_rate;		   ///< Num. events accumulated for "instantaneous" rate
double delta_time_single;		///< Time spent processing last event
double delta_time_rate;			///< Integrated time accumulated "instantaneous" rate (partial number of events)
double delta_time;				///< Total time spent processing events (this thread)
double rate_instantaneous;		///< Latest instantaneous rate
double rate_integrated;			///< Rate integrated over all events
 
static data_source_t null_data_source;

vector<pair<const char*, void*> > user_refs;
210};


213// The following is here just so we can use ROOT's THtml class to generate documentation.
214#ifdef G__DICTIONARY
215typedef JEventLoop::call_stack_t call_stack_t;
216typedef JEventLoop::error_call_stack_t error_call_stack_t;
217#endif

219#if !defined(__CINT__) && !defined(__CLING__)

221//-------------
222// GetSingle
223//-------------
224template<class T>
225JFactory<T>* JEventLoop::GetSingle(const T* &t, const char *tag, bool exception_if_not_one)
226{
/// This is a convenience method that can be used to get a pointer to the single
/// object of type T from the specified factory. It simply calls the Get(vector<...>) method
/// and copies the first pointer into "t" (or NULL if something other than 1 object is returned).
/// 
/// This is intended to address the common situation in which there is an interest
/// in the event if and only if there is exactly 1 object of type T. If the event
/// has no objects of that type or more than 1 object of that type (for the specified
/// factory) then an exception of type "unsigned long" is thrown with the value
/// being the number of objects of type T. You can supress the exception by setting
/// exception_if_not_one to false. In that case, you will have to check if t==NULL to
/// know if the call succeeded.
vector<const T*> v;
JFactory<T> *fac = Get(v, tag);

if(v.size()!=1){
t = NULL__null;
if(exception_if_not_one) throw v.size();
}

t = v[0];

return fac;
249}

251//-------------
252// Get
253//-------------
254template<class T> 
255JFactory<T>* JEventLoop::Get(vector<const T*> &t, const char *tag, bool allow_deftag)
256{
/// Retrieve or generate the array of objects of
/// type T for the curent event being processed
/// by this thread.
///
/// By default, preference is given to reading the
/// objects from the data source(e.g. file) before generating
/// them in the factory. A flag exists in the factory
/// however to change this so that the factory is
/// given preference.
///
/// Note that regardless of the setting of this flag,
/// the data are only either read in or generated once.
/// Ownership of the objects will always be with the
/// factory so subsequent calls will always return pointers to
/// the same data.
///
/// If the factory is called on to generate the data,
/// it is done by calling the factory's Get() method
/// which, in turn, calls the evnt() method. 
/// 
/// First, we just call the GetFromFactory() method.
/// It will make the initial decision as to whether
/// it should look in the source first or not. If
/// it returns NULL, then the factory couldn't be
/// found so we automatically try the file.
///
/// Note that if no factory exists to hold the objects
/// from the file, one can be created automatically
/// providing the <i>JANA:AUTOFACTORYCREATE</i>
/// configuration parameter is set.

// Check if a tag was specified for this data type to use for the
// default.
const char *mytag = tag1.1
'tag' is not equal to NULL
1
'tag' is not equal to NULL
1
'tag' is not equal to NULL
==NULL__null ? "":tag; // protection against NULL tags
2
←
'?' condition is false→
if(strlen(mytag)==0 && allow_deftag2.1
'allow_deftag' is true
1
'allow_deftag' is true
1
'allow_deftag' is true
){
3
←
Taking true branch→
map<string, string>::const_iterator iter = default_tags.find(T::static_className());
if(iter!=default_tags.end())tag = iter->second.c_str();
4
←
Assuming the condition is true→
5
←
Taking true branch→
6
←
Value assigned to 'tag'→
}


// If we are trying to keep track of the call stack then we
// need to add a new call_stack_t object to the the list
// and initialize it with the start time and caller/callee
// info.
call_stack_t cs;

// Optionally record starting info of call stack entry
if(record_call_stack) CallStackStart(cs, caller_name, caller_tag, T::static_className(), tag);
7
←
Assuming field 'record_call_stack' is false→
8
←
Taking false branch→

// Get the data (or at least try to)
JFactory<T>* factory=NULL__null;
try{
factory = GetFromFactory(t, tag, cs.data_source, allow_deftag);
9
←
Passing value via 2nd parameter 'tag'→
10
←
Calling 'JEventLoop::GetFromFactory'→
if(!factory){
	// No factory exists for this type and tag. It's possible
	// that the source may be able to provide the objects
	// but it will need a place to put them. We can create a
	// dumb JFactory just to hold the data in case the source
	// can provide the objects. Before we do though, make sure
	// the user condones this via the presence of the
	// "JANA:AUTOFACTORYCREATE" config parameter.
	string p;
	try{
		gPARMS->GetParameter("JANA:AUTOFACTORYCREATE", p);
	}catch(...){}
	if(p.size()==0){
		jout<<std::endl;
		_DBG__std::cerr<<"/w/halld-scifs17exp/halld2/home/sdobbs/Software/jana/jana_0.8.2/Linux_CentOS7.7-x86_64-gcc4.8.5/include/JANA/JEventLoop.h"
<<":"<<324<<std::endl;
		jout<<"No factory of type \""<<T::static_className()<<"\" with tag \""<<tag<<"\" exists."<<std::endl;
		jout<<"If you are reading objects from a file, I can auto-create a factory"<<std::endl;
		jout<<"of the appropriate type to hold the objects, but this feature is turned"<<std::endl;
		jout<<"off by default. To turn it on, set the \"JANA:AUTOFACTORYCREATE\""<<std::endl;
		jout<<"configuration parameter. This can usually be done by passing the"<<std::endl;
		jout<<"following argument to the program from the command line:"<<std::endl;
		jout<<std::endl;
		jout<<"   -PJANA:AUTOFACTORYCREATE=1"<<std::endl;
		jout<<std::endl;
		jout<<"Note that since the most commonly expected occurance of this situation."<<std::endl;
		jout<<"is an error, the program will now throw an exception so that the factory."<<std::endl;
		jout<<"call stack can be printed."<<std::endl;
		jout<<std::endl;
		throw exception();
	}else{
		AddFactory(new JFactory<T>(tag));
		jout<<__FILE__"/w/halld-scifs17exp/halld2/home/sdobbs/Software/jana/jana_0.8.2/Linux_CentOS7.7-x86_64-gcc4.8.5/include/JANA/JEventLoop.h"<<":"<<__LINE__341<<" Auto-created "<<T::static_className()<<":"<<tag<<" factory"<<std::endl;
	
		// Now try once more. The GetFromFactory method will call
		// GetFromSource since it's empty.
		factory = GetFromFactory(t, tag, cs.data_source, allow_deftag);
	}
}
}catch(exception &e){
// Uh-oh, an exception was thrown. Add us to the call stack
// and re-throw the exception
error_call_stack_t ecs;
ecs.factory_name = T::static_className();
ecs.tag = tag;
ecs.filename = NULL__null;
error_call_stack.push_back(ecs);
throw e;
}

// If recording the call stack, update the end_time field and add to stack
if(record_call_stack) CallStackEnd(cs);

return factory;
363}

365//-------------
366// GetFromFactory
367//-------------
368template<class T> 
369JFactory<T>* JEventLoop::GetFromFactory(vector<const T*> &t, const char *tag, data_source_t &data_source, bool allow_deftag)
370{
// We need to find the factory providing data type T with
// tag given by "tag". 
vector<JFactory_base*>::iterator iter=factories.begin();
JFactory<T> *factory = NULL__null;
string className(T::static_className());

// Check if a tag was specified for this data type to use for the
// default.
const char *mytag = tag==NULL__null ? "":tag; // protection against NULL tags
11
←
Assuming 'tag' is equal to NULL→
12
←
Assuming pointer value is null→
13
←
'?' condition is true→
if(strlen(mytag)==0 && allow_deftag13.1
'allow_deftag' is true
1
'allow_deftag' is true
1
'allow_deftag' is true
){
14
←
Taking true branch→
map<string, string>::const_iterator iter = default_tags.find(className);
if(iter!=default_tags.end())tag = iter->second.c_str();
15
←
Assuming the condition is false→
16
←
Taking false branch→
}

for(; iter!=factories.end(); iter++){
17
←
Calling 'operator!=<jana::JFactory_base **, std::vector<jana::JFactory_base *>>'→
20
←
Returning from 'operator!=<jana::JFactory_base **, std::vector<jana::JFactory_base *>>'→
21
←
Loop condition is true.  Entering loop body→
// It turns out a long standing bug in g++ makes dynamic_cast return
// zero improperly when used on objects created on one side of
// a dynamically shared object (DSO) and the cast occurs on the 
// other side. I saw bug reports ranging from 2001 to 2004. I saw
// saw it first-hand on LinuxEL4 using g++ 3.4.5. This is too bad
// since it is much more elegant (and safe) to use dynamic_cast.
// To avoid this problem which can occur with plugins, we check
// the name of the data classes are the same. (sigh)
//factory = dynamic_cast<JFactory<T> *>(*iter);
if(className == (*iter)->GetDataClassName())factory = (JFactory<T>*)(*iter);
22
←
Taking true branch→
if(factory == NULL__null)continue;
23
←
Assuming 'factory' is not equal to NULL→
24
←
Taking false branch→
const char *factag = factory->Tag()==NULL__null ? "":factory->Tag();
25
←
Assuming the condition is true→
26
←
'?' condition is true→
if(!strcmp(factag, tag)){
27
←
Null pointer passed to 2nd parameter expecting 'nonnull'
	break;
}else{
	factory=NULL__null;
}
}

// If factory not found, just return now
if(!factory){
data_source = DATA_NOT_AVAILABLE;
return NULL__null;
}

// OK, we found the factory. If the evnt() routine has already
// been called, then just call the factory's Get() routine
// to return a copy of the existing data
if(factory->evnt_was_called()){
factory->CopyFrom(t);
data_source = DATA_FROM_CACHE;
return factory;
}

// Next option is to get the objects from the data source
if(factory->GetCheckSourceFirst()){
// If the object type/tag is found in the source, it
// will return NOERROR, even if there are zero instances
// of it. If it is not available in the source then it
// will return OBJECT_NOT_AVAILABLE.

jerror_t err = GetFromSource(t, factory);
if(err == NOERROR){
	// A return value of NOERROR means the source had the objects
	// even if there were zero of them.(If the source had no
	// information about the objects OBJECT_NOT_AVAILABLE would 
	// have been returned.)
	// The GetFromSource() call will eventually lead to a call to
	// the GetObjects() method of the concrete class derived
	// from JEventSource. That routine should copy the object
	// pointers into the factory using the factory's CopyTo()
	// method which also sets the evnt_called flag for the factory.
	// Note also that the "t" vector is then filled with a call
	// to the factory's CopyFrom() method in JEvent::GetObjects().
	// All we need to do now is just set the factory pointers in
	// the newly generated JObjects and return the factory pointer.

	factory->SetFactoryPointers();
	data_source = DATA_FROM_SOURCE;

	return factory;
}
}

// OK. It looks like we have to have the factory make this.
// Get pointers to data from the factory.
factory->Get(t);
factory->SetFactoryPointers();
data_source = DATA_FROM_FACTORY;

return factory;
457}

459//-------------
460// GetFromSource
461//-------------
462template<class T> 
463jerror_t JEventLoop::GetFromSource(vector<const T*> &t, JFactory_base *factory)
464{
/// This tries to get objects from the event source.
/// "factory" must be a valid pointer to a JFactory
/// object since that will take ownership of the objects
/// created by the source.
/// This should usually be called from JEventLoop::GetFromFactory
/// which is called from JEventLoop::Get. The latter will
/// create a dummy JFactory of the proper flavor and tag if
/// one does not already exist so if objects exist in the
/// file without a corresponding factory to create them, they
/// can still be used.
if(!factory)throw OBJECT_NOT_AVAILABLE;

return event.GetObjects(t, factory);
478}

480//-------------
481// CallStackStart
482//-------------
483inline void JEventLoop::CallStackStart(JEventLoop::call_stack_t &cs, const string &caller_name, const string &caller_tag, const string callee_name, const string callee_tag)
484{
/// This is used to fill initial info into a call_stack_t stucture
/// for recording the call stack. It should be matched with a call
/// to CallStackEnd. It is normally called from the Get() method
/// above, but may also be used by external actors to manipulate
/// the call stack (presumably for good and not evil).

struct itimerval tmr;
getitimer(ITIMER_PROFITIMER_PROF, &tmr);

cs.caller_name    = this->caller_name;
cs.caller_tag     = this->caller_tag;
this->caller_name = cs.callee_name = callee_name;
this->caller_tag  = cs.callee_tag  = callee_tag;
cs.start_time = tmr.it_value.tv_sec + tmr.it_value.tv_usec/1.0E6;
499}

501//-------------
502// CallStackEnd
503//-------------
504inline void JEventLoop::CallStackEnd(JEventLoop::call_stack_t &cs)
505{
/// Complete a call stack entry. This should be matched
/// with a previous call to CallStackStart which was
/// used to fill the cs structure.

struct itimerval tmr;
getitimer(ITIMER_PROFITIMER_PROF, &tmr);
cs.end_time = tmr.it_value.tv_sec + tmr.it_value.tv_usec/1.0E6;
caller_name = cs.caller_name;
caller_tag  = cs.caller_tag;
call_stack.push_back(cs);
516}

518//-------------
519// CheckEventBoundary
520//-------------
521inline bool JEventLoop::CheckEventBoundary(uint64_t event_numberA, uint64_t event_numberB)
522{
/// Check whether the two event numbers span one or more boundaries
/// in the calibration/conditions database for the current run number.
/// Return true if they do and false if they don't. The first parameter
/// "event_numberA" is also checked if it lands on a boundary in which
/// case true is also returned. If event_numberB lands on a boundary,
/// but event_numberA does not, then false is returned.
///
/// This method is not expected to be called by a user. It is, however called,
/// everytime a JFactory's Get() method is called.

// Make sure our copy of the boundaries is up to date
if(event.GetRunNumber()!=event_boundaries_run){
event_boundaries.clear(); // in case we can't get the JCalibration pointer
JCalibration *jcalib = GetJCalibration();
if(jcalib)jcalib->GetEventBoundaries(event_boundaries);
event_boundaries_run = event.GetRunNumber();
}

// Loop over boundaries
for(unsigned int i=0; i<event_boundaries.size(); i++){
uint64_t eb = event_boundaries[i];
if((eb - event_numberA)*(eb - event_numberB) < 0.0 || eb==event_numberA){ // think about it ....
	// events span a boundary or is on a boundary. Return true
	return true;
}
}

return false;
551}

553//-------------
554// FindByID
555//-------------
556template<class T> 
557const T* JEventLoop::FindByID(JObject::oid_t id)
558{
/// This is a templated method that can be used in place
/// of the non-templated FindByID(oid_t) method if one knows
/// the class of the object with the specified id.
/// This method is faster than calling the non-templated
/// FindByID and dynamic_cast-ing the JObject since
/// this will only search the objects of factories that
/// produce the desired data type.
/// This method will cast the JObject pointer to one
/// of the specified type. To use this method,
/// a type is specified in the call as follows:
///
/// const DMyType *t = loop->FindByID<DMyType>(id);

// Loop over factories looking for ones that provide
// specified data type.
for(unsigned int i=0; i<factories.size(); i++){
if(factories[i]->GetDataClassName() != T::static_className())continue;

// This factory provides data of type T. Search it for
// the object with the specified id.
const JObject *my_obj = factories[i]->GetByID(id);
if(my_obj)return dynamic_cast<const T*>(my_obj);
}

return NULL__null;
584}

586//-------------
587// GetCalib (map)
588//-------------
589template<class T>
590bool JEventLoop::GetCalib(string namepath, map<string,T> &vals)
591{
/// Get the JCalibration object from JApplication for the run number of
/// the current event and call its Get() method to get the constants.

// Note that we could do this by making "vals" a generic type T thus, combining
// this with the vector version below. However, doing this explicitly will make
// it easier for the user to understand how to call us.

vals.clear();

JCalibration *calib = GetJCalibration();
if(!calib){
_DBG_std::cerr<<"/w/halld-scifs17exp/halld2/home/sdobbs/Software/jana/jana_0.8.2/Linux_CentOS7.7-x86_64-gcc4.8.5/include/JANA/JEventLoop.h"
<<":"<<603<<" "<<"Unable to get JCalibration object for run "<<event.GetRunNumber()<<std::endl;
return true;
}

return calib->Get(namepath, vals, event.GetEventNumber());
608}

610//-------------
611// GetCalib (vector)
612//-------------
613template<class T> bool JEventLoop::GetCalib(string namepath, vector<T> &vals)
614{
/// Get the JCalibration object from JApplication for the run number of
/// the current event and call its Get() method to get the constants.

vals.clear();

JCalibration *calib = GetJCalibration();
if(!calib){
_DBG_std::cerr<<"/w/halld-scifs17exp/halld2/home/sdobbs/Software/jana/jana_0.8.2/Linux_CentOS7.7-x86_64-gcc4.8.5/include/JANA/JEventLoop.h"
<<":"<<622<<" "<<"Unable to get JCalibration object for run "<<event.GetRunNumber()<<std::endl;
return true;
}

return calib->Get(namepath, vals, event.GetEventNumber());
627}

629//-------------
630// GetCalib (single)
631//-------------
632template<class T> bool JEventLoop::GetCalib(string namepath, T &val)
633{
/// This is a convenience method for getting a single entry. It
/// simply calls the vector version and returns the first entry.
/// It returns true if the vector version returns true AND there
/// is at least one entry in the vector. No check is made for there
/// there being more than one entry in the vector.

vector<T> vals;
bool ret = GetCalib(namepath, vals);
if(vals.empty()) return true;
val = vals[0];

return ret;
646}

648//-------------
649// GetGeom (map)
650//-------------
651template<class T>
652bool JEventLoop::GetGeom(string namepath, map<string,T> &vals)
653{
/// Get the JGeometry object from JApplication for the run number of
/// the current event and call its Get() method to get the constants.

// Note that we could do this by making "vals" a generic type T thus, combining
// this with the vector version below. However, doing this explicitly will make
// it easier for the user to understand how to call us.

vals.clear();

JGeometry *geom = GetJGeometry();
if(!geom){
_DBG_std::cerr<<"/w/halld-scifs17exp/halld2/home/sdobbs/Software/jana/jana_0.8.2/Linux_CentOS7.7-x86_64-gcc4.8.5/include/JANA/JEventLoop.h"
<<":"<<665<<" "<<"Unable to get JGeometry object for run "<<event.GetRunNumber()<<std::endl;
return true;
}

return geom->Get(namepath, vals);
670}

672//-------------
673// GetGeom (atomic)
674//-------------
675template<class T> bool JEventLoop::GetGeom(string namepath, T &val)
676{
/// Get the JCalibration object from JApplication for the run number of
/// the current event and call its Get() method to get the constants.

JGeometry *geom = GetJGeometry();
if(!geom){
_DBG_std::cerr<<"/w/halld-scifs17exp/halld2/home/sdobbs/Software/jana/jana_0.8.2/Linux_CentOS7.7-x86_64-gcc4.8.5/include/JANA/JEventLoop.h"
<<":"<<682<<" "<<"Unable to get JGeometry object for run "<<event.GetRunNumber()<<std::endl;
return true;
}

return geom->Get(namepath, val);
687}

689//-------------
690// SetRef
691//-------------
692template<class T>
693void JEventLoop::SetRef(T *t)
694{
pair<const char*, void*> p(typeid(T).name(), (void*)t);
user_refs.push_back(p);
697}

699//-------------
700// GetResource
701//-------------
702template<class T> bool JEventLoop::GetResource(string namepath, T vals, int event_number)
703{
JResourceManager *resource_manager = GetJResourceManager();
if(!resource_manager){
string mess = string("Unable to get the JResourceManager object (namepath=\"")+namepath+"\")";
throw JException(mess);
}

return resource_manager->Get(namepath, vals, event_number);
711}

713//-------------
714// GetRef
715//-------------
716template<class T>
717T* JEventLoop::GetRef(void)
718{
/// Get a user-defined reference (a pointer)
for(unsigned int i=0; i<user_refs.size(); i++){
if(user_refs[i].first == typeid(T).name()) return (T*)user_refs[i].second;
}

return NULL__null;
725}

727//-------------
728// GetRefsT
729//-------------
730template<class T>
731vector<T*> JEventLoop::GetRefsT(void)
732{
vector<T*> refs;
for(unsigned int i=0; i<user_refs.size(); i++){
if(user_refs[i].first == typeid(T).name()){
	refs.push_back((T*)user_refs[i].second);
}
}

return refs;
741}

743//-------------
744// RemoveRef
745//-------------
746template<class T>
747void JEventLoop::RemoveRef(T *t)
748{
vector<pair<const char*, void*> >::iterator iter;
for(iter=user_refs.begin(); iter!= user_refs.end(); iter++){
if(iter->second == (void*)t){
	user_refs.erase(iter);
	return;
}
}
_DBG_std::cerr<<"/w/halld-scifs17exp/halld2/home/sdobbs/Software/jana/jana_0.8.2/Linux_CentOS7.7-x86_64-gcc4.8.5/include/JANA/JEventLoop.h"
<<":"<<756<<" "<<" Attempt to remove user reference not in event loop!" << std::endl;
757}


760#endif //__CINT__  __CLING__

762} // Close JANA namespace



766#endif // _JEventLoop_


←

/usr/lib/gcc/x86_64-redhat-linux/4.8.5/../../../../include/c++/4.8.5/bits/stl_iterator.h

1// Iterators -*- C++ -*-

3// Copyright (C) 2001-2013 Free Software Foundation, Inc.
4//
5// This file is part of the GNU ISO C++ Library.  This library is free
6// software; you can redistribute it and/or modify it under the
7// terms of the GNU General Public License as published by the
8// Free Software Foundation; either version 3, or (at your option)
9// any later version.

11// This library is distributed in the hope that it will be useful,
12// but WITHOUT ANY WARRANTY; without even the implied warranty of
13// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
14// GNU General Public License for more details.

16// Under Section 7 of GPL version 3, you are granted additional
17// permissions described in the GCC Runtime Library Exception, version
18// 3.1, as published by the Free Software Foundation.

20// You should have received a copy of the GNU General Public License and
21// a copy of the GCC Runtime Library Exception along with this program;
22// see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
23// <http://www.gnu.org/licenses/>.

25/*
*
* Copyright (c) 1994
* Hewlett-Packard Company
*
* Permission to use, copy, modify, distribute and sell this software
* and its documentation for any purpose is hereby granted without fee,
* provided that the above copyright notice appear in all copies and
* that both that copyright notice and this permission notice appear
* in supporting documentation.  Hewlett-Packard Company makes no
* representations about the suitability of this software for any
* purpose.  It is provided "as is" without express or implied warranty.
*
*
* Copyright (c) 1996-1998
* Silicon Graphics Computer Systems, Inc.
*
* Permission to use, copy, modify, distribute and sell this software
* and its documentation for any purpose is hereby granted without fee,
* provided that the above copyright notice appear in all copies and
* that both that copyright notice and this permission notice appear
* in supporting documentation.  Silicon Graphics makes no
* representations about the suitability of this software for any
* purpose.  It is provided "as is" without express or implied warranty.
*/

51/** @file bits/stl_iterator.h
*  This is an internal header file, included by other library headers.
*  Do not attempt to use it directly. @headername{iterator}
*
*  This file implements reverse_iterator, back_insert_iterator,
*  front_insert_iterator, insert_iterator, __normal_iterator, and their
*  supporting functions and overloaded operators.
*/

60#ifndef _STL_ITERATOR_H1
61#define _STL_ITERATOR_H1 1

63#include <bits/cpp_type_traits.h>
64#include <ext/type_traits.h>
65#include <bits/move.h>

67namespace std _GLIBCXX_VISIBILITY(default)__attribute__ ((__visibility__ ("default")))
68{
69_GLIBCXX_BEGIN_NAMESPACE_VERSION

/**
 * @addtogroup iterators
 * @{
 */

// 24.4.1 Reverse iterators
/**
 *  Bidirectional and random access iterators have corresponding reverse
 *  %iterator adaptors that iterate through the data structure in the
 *  opposite direction.  They have the same signatures as the corresponding
 *  iterators.  The fundamental relation between a reverse %iterator and its
 *  corresponding %iterator @c i is established by the identity:
 *  @code
 *      &*(reverse_iterator(i)) == &*(i - 1)
 *  @endcode
 *
 *  <em>This mapping is dictated by the fact that while there is always a
 *  pointer past the end of an array, there might not be a valid pointer
 *  before the beginning of an array.</em> [24.4.1]/1,2
 *
 *  Reverse iterators can be tricky and surprising at first.  Their
 *  semantics make sense, however, and the trickiness is a side effect of
 *  the requirement that the iterators must be safe.
*/
template<typename _Iterator>
  class reverse_iterator
  : public iterator<typename iterator_traits<_Iterator>::iterator_category,
      typename iterator_traits<_Iterator>::value_type,
      typename iterator_traits<_Iterator>::difference_type,
      typename iterator_traits<_Iterator>::pointer,
                    typename iterator_traits<_Iterator>::reference>
  {
  protected:
    _Iterator current;

    typedef iterator_traits<_Iterator>		__traits_type;

  public:
    typedef _Iterator					iterator_type;
    typedef typename __traits_type::difference_type	difference_type;
    typedef typename __traits_type::pointer		pointer;
    typedef typename __traits_type::reference		reference;

    /**
     *  The default constructor value-initializes member @p current.
     *  If it is a pointer, that means it is zero-initialized.
    */
    // _GLIBCXX_RESOLVE_LIB_DEFECTS
    // 235 No specification of default ctor for reverse_iterator
    reverse_iterator() : current() { }

    /**
     *  This %iterator will move in the opposite direction that @p x does.
    */
    explicit
    reverse_iterator(iterator_type __x) : current(__x) { }

    /**
     *  The copy constructor is normal.
    */
    reverse_iterator(const reverse_iterator& __x)
    : current(__x.current) { }

    /**
     *  A %reverse_iterator across other types can be copied if the
     *  underlying %iterator can be converted to the type of @c current.
    */
    template<typename _Iter>
      reverse_iterator(const reverse_iterator<_Iter>& __x)
: current(__x.base()) { }

    /**
     *  @return  @c current, the %iterator used for underlying work.
    */
    iterator_type
    base() const
    { return current; }

    /**
     *  @return  A reference to the value at @c --current
     *
     *  This requires that @c --current is dereferenceable.
     *
     *  @warning This implementation requires that for an iterator of the
     *           underlying iterator type, @c x, a reference obtained by
     *           @c *x remains valid after @c x has been modified or
     *           destroyed. This is a bug: http://gcc.gnu.org/PR51823
    */
    reference
    operator*() const
    {
_Iterator __tmp = current;
return *--__tmp;
    }

    /**
     *  @return  A pointer to the value at @c --current
     *
     *  This requires that @c --current is dereferenceable.
    */
    pointer
    operator->() const
    { return &(operator*()); }

    /**
     *  @return  @c *this
     *
     *  Decrements the underlying iterator.
    */
    reverse_iterator&
    operator++()
    {
--current;
return *this;
    }

    /**
     *  @return  The original value of @c *this
     *
     *  Decrements the underlying iterator.
    */
    reverse_iterator
    operator++(int)
    {
reverse_iterator __tmp = *this;
--current;
return __tmp;
    }

    /**
     *  @return  @c *this
     *
     *  Increments the underlying iterator.
    */
    reverse_iterator&
    operator--()
    {
++current;
return *this;
    }

    /**
     *  @return  A reverse_iterator with the previous value of @c *this
     *
     *  Increments the underlying iterator.
    */
    reverse_iterator
    operator--(int)
    {
reverse_iterator __tmp = *this;
++current;
return __tmp;
    }

    /**
     *  @return  A reverse_iterator that refers to @c current - @a __n
     *
     *  The underlying iterator must be a Random Access Iterator.
    */
    reverse_iterator
    operator+(difference_type __n) const
    { return reverse_iterator(current - __n); }

    /**
     *  @return  *this
     *
     *  Moves the underlying iterator backwards @a __n steps.
     *  The underlying iterator must be a Random Access Iterator.
    */
    reverse_iterator&
    operator+=(difference_type __n)
    {
current -= __n;
return *this;
    }

    /**
     *  @return  A reverse_iterator that refers to @c current - @a __n
     *
     *  The underlying iterator must be a Random Access Iterator.
    */
    reverse_iterator
    operator-(difference_type __n) const
    { return reverse_iterator(current + __n); }

    /**
     *  @return  *this
     *
     *  Moves the underlying iterator forwards @a __n steps.
     *  The underlying iterator must be a Random Access Iterator.
    */
    reverse_iterator&
    operator-=(difference_type __n)
    {
current += __n;
return *this;
    }

    /**
     *  @return  The value at @c current - @a __n - 1
     *
     *  The underlying iterator must be a Random Access Iterator.
    */
    reference
    operator[](difference_type __n) const
    { return *(*this + __n); }
  };

//@{
/**
 *  @param  __x  A %reverse_iterator.
 *  @param  __y  A %reverse_iterator.
 *  @return  A simple bool.
 *
 *  Reverse iterators forward many operations to their underlying base()
 *  iterators.  Others are implemented in terms of one another.
 *
*/
template<typename _Iterator>
  inline bool
  operator==(const reverse_iterator<_Iterator>& __x,
      const reverse_iterator<_Iterator>& __y)
  { return __x.base() == __y.base(); }

template<typename _Iterator>
  inline bool
  operator<(const reverse_iterator<_Iterator>& __x,
     const reverse_iterator<_Iterator>& __y)
  { return __y.base() < __x.base(); }

template<typename _Iterator>
  inline bool
  operator!=(const reverse_iterator<_Iterator>& __x,
      const reverse_iterator<_Iterator>& __y)
  { return !(__x == __y); }

template<typename _Iterator>
  inline bool
  operator>(const reverse_iterator<_Iterator>& __x,
     const reverse_iterator<_Iterator>& __y)
  { return __y < __x; }

template<typename _Iterator>
  inline bool
  operator<=(const reverse_iterator<_Iterator>& __x,
      const reverse_iterator<_Iterator>& __y)
  { return !(__y < __x); }

template<typename _Iterator>
  inline bool
  operator>=(const reverse_iterator<_Iterator>& __x,
      const reverse_iterator<_Iterator>& __y)
  { return !(__x < __y); }

template<typename _Iterator>
  inline typename reverse_iterator<_Iterator>::difference_type
  operator-(const reverse_iterator<_Iterator>& __x,
     const reverse_iterator<_Iterator>& __y)
  { return __y.base() - __x.base(); }

template<typename _Iterator>
  inline reverse_iterator<_Iterator>
  operator+(typename reverse_iterator<_Iterator>::difference_type __n,
     const reverse_iterator<_Iterator>& __x)
  { return reverse_iterator<_Iterator>(__x.base() - __n); }

// _GLIBCXX_RESOLVE_LIB_DEFECTS
// DR 280. Comparison of reverse_iterator to const reverse_iterator.
template<typename _IteratorL, typename _IteratorR>
  inline bool
  operator==(const reverse_iterator<_IteratorL>& __x,
      const reverse_iterator<_IteratorR>& __y)
  { return __x.base() == __y.base(); }

template<typename _IteratorL, typename _IteratorR>
  inline bool
  operator<(const reverse_iterator<_IteratorL>& __x,
     const reverse_iterator<_IteratorR>& __y)
  { return __y.base() < __x.base(); }

template<typename _IteratorL, typename _IteratorR>
  inline bool
  operator!=(const reverse_iterator<_IteratorL>& __x,
      const reverse_iterator<_IteratorR>& __y)
  { return !(__x == __y); }

template<typename _IteratorL, typename _IteratorR>
  inline bool
  operator>(const reverse_iterator<_IteratorL>& __x,
     const reverse_iterator<_IteratorR>& __y)
  { return __y < __x; }

template<typename _IteratorL, typename _IteratorR>
  inline bool
  operator<=(const reverse_iterator<_IteratorL>& __x,
      const reverse_iterator<_IteratorR>& __y)
  { return !(__y < __x); }

template<typename _IteratorL, typename _IteratorR>
  inline bool
  operator>=(const reverse_iterator<_IteratorL>& __x,
      const reverse_iterator<_IteratorR>& __y)
  { return !(__x < __y); }

template<typename _IteratorL, typename _IteratorR>
376#if __cplusplus201103L >= 201103L
  // DR 685.
  inline auto
  operator-(const reverse_iterator<_IteratorL>& __x,
     const reverse_iterator<_IteratorR>& __y)
  -> decltype(__y.base() - __x.base())
382#else
  inline typename reverse_iterator<_IteratorL>::difference_type
  operator-(const reverse_iterator<_IteratorL>& __x,
     const reverse_iterator<_IteratorR>& __y)
386#endif
  { return __y.base() - __x.base(); }
//@}

// 24.4.2.2.1 back_insert_iterator
/**
 *  @brief  Turns assignment into insertion.
 *
 *  These are output iterators, constructed from a container-of-T.
 *  Assigning a T to the iterator appends it to the container using
 *  push_back.
 *
 *  Tip:  Using the back_inserter function to create these iterators can
 *  save typing.
*/
template<typename _Container>
  class back_insert_iterator
  : public iterator<output_iterator_tag, void, void, void, void>
  {
  protected:
    _Container* container;

  public:
    /// A nested typedef for the type of whatever container you used.
    typedef _Container          container_type;

    /// The only way to create this %iterator is with a container.
    explicit
    back_insert_iterator(_Container& __x) : container(&__x) { }

    /**
     *  @param  __value  An instance of whatever type
     *                 container_type::const_reference is; presumably a
     *                 reference-to-const T for container<T>.
     *  @return  This %iterator, for chained operations.
     *
     *  This kind of %iterator doesn't really have a @a position in the
     *  container (you can think of the position as being permanently at
     *  the end, if you like).  Assigning a value to the %iterator will
     *  always append the value to the end of the container.
    */
427#if __cplusplus201103L < 201103L
    back_insert_iterator&
    operator=(typename _Container::const_reference __value)
    {
container->push_back(__value);
return *this;
    }
434#else
    back_insert_iterator&
    operator=(const typename _Container::value_type& __value)
    {
container->push_back(__value);
return *this;
    }

    back_insert_iterator&
    operator=(typename _Container::value_type&& __value)
    {
container->push_back(std::move(__value));
return *this;
    }
448#endif

    /// Simply returns *this.
    back_insert_iterator&
    operator*()
    { return *this; }

    /// Simply returns *this.  (This %iterator does not @a move.)
    back_insert_iterator&
    operator++()
    { return *this; }

    /// Simply returns *this.  (This %iterator does not @a move.)
    back_insert_iterator
    operator++(int)
    { return *this; }
  };

/**
 *  @param  __x  A container of arbitrary type.
 *  @return  An instance of back_insert_iterator working on @p __x.
 *
 *  This wrapper function helps in creating back_insert_iterator instances.
 *  Typing the name of the %iterator requires knowing the precise full
 *  type of the container, which can be tedious and impedes generic
 *  programming.  Using this function lets you take advantage of automatic
 *  template parameter deduction, making the compiler match the correct
 *  types for you.
*/
template<typename _Container>
  inline back_insert_iterator<_Container>
  back_inserter(_Container& __x)
  { return back_insert_iterator<_Container>(__x); }

/**
 *  @brief  Turns assignment into insertion.
 *
 *  These are output iterators, constructed from a container-of-T.
 *  Assigning a T to the iterator prepends it to the container using
 *  push_front.
 *
 *  Tip:  Using the front_inserter function to create these iterators can
 *  save typing.
*/
template<typename _Container>
  class front_insert_iterator
  : public iterator<output_iterator_tag, void, void, void, void>
  {
  protected:
    _Container* container;

  public:
    /// A nested typedef for the type of whatever container you used.
    typedef _Container          container_type;

    /// The only way to create this %iterator is with a container.
    explicit front_insert_iterator(_Container& __x) : container(&__x) { }

    /**
     *  @param  __value  An instance of whatever type
     *                 container_type::const_reference is; presumably a
     *                 reference-to-const T for container<T>.
     *  @return  This %iterator, for chained operations.
     *
     *  This kind of %iterator doesn't really have a @a position in the
     *  container (you can think of the position as being permanently at
     *  the front, if you like).  Assigning a value to the %iterator will
     *  always prepend the value to the front of the container.
    */
517#if __cplusplus201103L < 201103L
    front_insert_iterator&
    operator=(typename _Container::const_reference __value)
    {
container->push_front(__value);
return *this;
    }
524#else
    front_insert_iterator&
    operator=(const typename _Container::value_type& __value)
    {
container->push_front(__value);
return *this;
    }

    front_insert_iterator&
    operator=(typename _Container::value_type&& __value)
    {
container->push_front(std::move(__value));
return *this;
    }
538#endif

    /// Simply returns *this.
    front_insert_iterator&
    operator*()
    { return *this; }

    /// Simply returns *this.  (This %iterator does not @a move.)
    front_insert_iterator&
    operator++()
    { return *this; }

    /// Simply returns *this.  (This %iterator does not @a move.)
    front_insert_iterator
    operator++(int)
    { return *this; }
  };

/**
 *  @param  __x  A container of arbitrary type.
 *  @return  An instance of front_insert_iterator working on @p x.
 *
 *  This wrapper function helps in creating front_insert_iterator instances.
 *  Typing the name of the %iterator requires knowing the precise full
 *  type of the container, which can be tedious and impedes generic
 *  programming.  Using this function lets you take advantage of automatic
 *  template parameter deduction, making the compiler match the correct
 *  types for you.
*/
template<typename _Container>
  inline front_insert_iterator<_Container>
  front_inserter(_Container& __x)
  { return front_insert_iterator<_Container>(__x); }

/**
 *  @brief  Turns assignment into insertion.
 *
 *  These are output iterators, constructed from a container-of-T.
 *  Assigning a T to the iterator inserts it in the container at the
 *  %iterator's position, rather than overwriting the value at that
 *  position.
 *
 *  (Sequences will actually insert a @e copy of the value before the
 *  %iterator's position.)
 *
 *  Tip:  Using the inserter function to create these iterators can
 *  save typing.
*/
template<typename _Container>
  class insert_iterator
  : public iterator<output_iterator_tag, void, void, void, void>
  {
  protected:
    _Container* container;
    typename _Container::iterator iter;

  public:
    /// A nested typedef for the type of whatever container you used.
    typedef _Container          container_type;

    /**
     *  The only way to create this %iterator is with a container and an
     *  initial position (a normal %iterator into the container).
    */
    insert_iterator(_Container& __x, typename _Container::iterator __i)
    : container(&__x), iter(__i) {}

    /**
     *  @param  __value  An instance of whatever type
     *                 container_type::const_reference is; presumably a
     *                 reference-to-const T for container<T>.
     *  @return  This %iterator, for chained operations.
     *
     *  This kind of %iterator maintains its own position in the
     *  container.  Assigning a value to the %iterator will insert the
     *  value into the container at the place before the %iterator.
     *
     *  The position is maintained such that subsequent assignments will
     *  insert values immediately after one another.  For example,
     *  @code
     *     // vector v contains A and Z
     *
     *     insert_iterator i (v, ++v.begin());
     *     i = 1;
     *     i = 2;
     *     i = 3;
     *
     *     // vector v contains A, 1, 2, 3, and Z
     *  @endcode
    */
628#if __cplusplus201103L < 201103L
    insert_iterator&
    operator=(typename _Container::const_reference __value)
    {
iter = container->insert(iter, __value);
++iter;
return *this;
    }
636#else
    insert_iterator&
    operator=(const typename _Container::value_type& __value)
    {
iter = container->insert(iter, __value);
++iter;
return *this;
    }

    insert_iterator&
    operator=(typename _Container::value_type&& __value)
    {
iter = container->insert(iter, std::move(__value));
++iter;
return *this;
    }
652#endif

    /// Simply returns *this.
    insert_iterator&
    operator*()
    { return *this; }

    /// Simply returns *this.  (This %iterator does not @a move.)
    insert_iterator&
    operator++()
    { return *this; }

    /// Simply returns *this.  (This %iterator does not @a move.)
    insert_iterator&
    operator++(int)
    { return *this; }
  };

/**
 *  @param __x  A container of arbitrary type.
 *  @return  An instance of insert_iterator working on @p __x.
 *
 *  This wrapper function helps in creating insert_iterator instances.
 *  Typing the name of the %iterator requires knowing the precise full
 *  type of the container, which can be tedious and impedes generic
 *  programming.  Using this function lets you take advantage of automatic
 *  template parameter deduction, making the compiler match the correct
 *  types for you.
*/
template<typename _Container, typename _Iterator>
  inline insert_iterator<_Container>
  inserter(_Container& __x, _Iterator __i)
  {
    return insert_iterator<_Container>(__x,
			 typename _Container::iterator(__i));
  }

// @} group iterators

691_GLIBCXX_END_NAMESPACE_VERSION
692} // namespace

694namespace __gnu_cxx _GLIBCXX_VISIBILITY(default)__attribute__ ((__visibility__ ("default")))
695{
696_GLIBCXX_BEGIN_NAMESPACE_VERSION

// This iterator adapter is @a normal in the sense that it does not
// change the semantics of any of the operators of its iterator
// parameter.  Its primary purpose is to convert an iterator that is
// not a class, e.g. a pointer, into an iterator that is a class.
// The _Container parameter exists solely so that different containers
// using this template can instantiate different types, even if the
// _Iterator parameter is the same.
using std::iterator_traits;
using std::iterator;
template<typename _Iterator, typename _Container>
  class __normal_iterator
  {
  protected:
    _Iterator _M_current;

    typedef iterator_traits<_Iterator>		__traits_type;

  public:
    typedef _Iterator					iterator_type;
    typedef typename __traits_type::iterator_category iterator_category;
    typedef typename __traits_type::value_type  	value_type;
    typedef typename __traits_type::difference_type 	difference_type;
    typedef typename __traits_type::reference 	reference;
    typedef typename __traits_type::pointer   	pointer;

    _GLIBCXX_CONSTEXPRconstexpr __normal_iterator() : _M_current(_Iterator()) { }

    explicit
    __normal_iterator(const _Iterator& __i) : _M_current(__i) { }

    // Allow iterator to const_iterator conversion
    template<typename _Iter>
      __normal_iterator(const __normal_iterator<_Iter,
	  typename __enable_if<
    	       (std::__are_same<_Iter, typename _Container::pointer>::__value),
      _Container>::__type>& __i)
      : _M_current(__i.base()) { }

    // Forward iterator requirements
    reference
    operator*() const
    { return *_M_current; }

    pointer
    operator->() const
    { return _M_current; }

    __normal_iterator&
    operator++()
    {
++_M_current;
return *this;
    }

    __normal_iterator
    operator++(int)
    { return __normal_iterator(_M_current++); }

    // Bidirectional iterator requirements
    __normal_iterator&
    operator--()
    {
--_M_current;
return *this;
    }

    __normal_iterator
    operator--(int)
    { return __normal_iterator(_M_current--); }

    // Random access iterator requirements
    reference
    operator[](const difference_type& __n) const
    { return _M_current[__n]; }

    __normal_iterator&
    operator+=(const difference_type& __n)
    { _M_current += __n; return *this; }

    __normal_iterator
    operator+(const difference_type& __n) const
    { return __normal_iterator(_M_current + __n); }

    __normal_iterator&
    operator-=(const difference_type& __n)
    { _M_current -= __n; return *this; }

    __normal_iterator
    operator-(const difference_type& __n) const
    { return __normal_iterator(_M_current - __n); }

    const _Iterator&
    base() const
    { return _M_current; }
  };

// Note: In what follows, the left- and right-hand-side iterators are
// allowed to vary in types (conceptually in cv-qualification) so that
// comparison between cv-qualified and non-cv-qualified iterators be
// valid.  However, the greedy and unfriendly operators in std::rel_ops
// will make overload resolution ambiguous (when in scope) if we don't
// provide overloads whose operands are of the same type.  Can someone
// remind me what generic programming is about? -- Gaby

// Forward iterator requirements
template<typename _IteratorL, typename _IteratorR, typename _Container>
  inline bool
  operator==(const __normal_iterator<_IteratorL, _Container>& __lhs,
      const __normal_iterator<_IteratorR, _Container>& __rhs)
  { return __lhs.base() == __rhs.base(); }

template<typename _Iterator, typename _Container>
  inline bool
  operator==(const __normal_iterator<_Iterator, _Container>& __lhs,
      const __normal_iterator<_Iterator, _Container>& __rhs)
  { return __lhs.base() == __rhs.base(); }

template<typename _IteratorL, typename _IteratorR, typename _Container>
  inline bool
  operator!=(const __normal_iterator<_IteratorL, _Container>& __lhs,
      const __normal_iterator<_IteratorR, _Container>& __rhs)
  { return __lhs.base() != __rhs.base(); }

template<typename _Iterator, typename _Container>
  inline bool
  operator!=(const __normal_iterator<_Iterator, _Container>& __lhs,
      const __normal_iterator<_Iterator, _Container>& __rhs)
  { return __lhs.base() != __rhs.base(); }
18
←
Assuming the condition is true→
19
←
Returning the value 1, which participates in a condition later→

// Random access iterator requirements
template<typename _IteratorL, typename _IteratorR, typename _Container>
  inline bool
  operator<(const __normal_iterator<_IteratorL, _Container>& __lhs,
     const __normal_iterator<_IteratorR, _Container>& __rhs)
  { return __lhs.base() < __rhs.base(); }

template<typename _Iterator, typename _Container>
  inline bool
  operator<(const __normal_iterator<_Iterator, _Container>& __lhs,
     const __normal_iterator<_Iterator, _Container>& __rhs)
  { return __lhs.base() < __rhs.base(); }

template<typename _IteratorL, typename _IteratorR, typename _Container>
  inline bool
  operator>(const __normal_iterator<_IteratorL, _Container>& __lhs,
     const __normal_iterator<_IteratorR, _Container>& __rhs)
  { return __lhs.base() > __rhs.base(); }

template<typename _Iterator, typename _Container>
  inline bool
  operator>(const __normal_iterator<_Iterator, _Container>& __lhs,
     const __normal_iterator<_Iterator, _Container>& __rhs)
  { return __lhs.base() > __rhs.base(); }

template<typename _IteratorL, typename _IteratorR, typename _Container>
  inline bool
  operator<=(const __normal_iterator<_IteratorL, _Container>& __lhs,
      const __normal_iterator<_IteratorR, _Container>& __rhs)
  { return __lhs.base() <= __rhs.base(); }

template<typename _Iterator, typename _Container>
  inline bool
  operator<=(const __normal_iterator<_Iterator, _Container>& __lhs,
      const __normal_iterator<_Iterator, _Container>& __rhs)
  { return __lhs.base() <= __rhs.base(); }

template<typename _IteratorL, typename _IteratorR, typename _Container>
  inline bool
  operator>=(const __normal_iterator<_IteratorL, _Container>& __lhs,
      const __normal_iterator<_IteratorR, _Container>& __rhs)
  { return __lhs.base() >= __rhs.base(); }

template<typename _Iterator, typename _Container>
  inline bool
  operator>=(const __normal_iterator<_Iterator, _Container>& __lhs,
      const __normal_iterator<_Iterator, _Container>& __rhs)
  { return __lhs.base() >= __rhs.base(); }

// _GLIBCXX_RESOLVE_LIB_DEFECTS
// According to the resolution of DR179 not only the various comparison
// operators but also operator- must accept mixed iterator/const_iterator
// parameters.
template<typename _IteratorL, typename _IteratorR, typename _Container>
881#if __cplusplus201103L >= 201103L
  // DR 685.
  inline auto
  operator-(const __normal_iterator<_IteratorL, _Container>& __lhs,
     const __normal_iterator<_IteratorR, _Container>& __rhs)
  -> decltype(__lhs.base() - __rhs.base())
887#else
  inline typename __normal_iterator<_IteratorL, _Container>::difference_type
  operator-(const __normal_iterator<_IteratorL, _Container>& __lhs,
     const __normal_iterator<_IteratorR, _Container>& __rhs)
891#endif
  { return __lhs.base() - __rhs.base(); }

template<typename _Iterator, typename _Container>
  inline typename __normal_iterator<_Iterator, _Container>::difference_type
  operator-(const __normal_iterator<_Iterator, _Container>& __lhs,
     const __normal_iterator<_Iterator, _Container>& __rhs)
  { return __lhs.base() - __rhs.base(); }

template<typename _Iterator, typename _Container>
  inline __normal_iterator<_Iterator, _Container>
  operator+(typename __normal_iterator<_Iterator, _Container>::difference_type
     __n, const __normal_iterator<_Iterator, _Container>& __i)
  { return __normal_iterator<_Iterator, _Container>(__i.base() + __n); }

906_GLIBCXX_END_NAMESPACE_VERSION
907} // namespace

909#if __cplusplus201103L >= 201103L

911namespace std _GLIBCXX_VISIBILITY(default)__attribute__ ((__visibility__ ("default")))
912{
913_GLIBCXX_BEGIN_NAMESPACE_VERSION

/**
 * @addtogroup iterators
 * @{
 */

// 24.4.3  Move iterators
/**
 *  Class template move_iterator is an iterator adapter with the same
 *  behavior as the underlying iterator except that its dereference
 *  operator implicitly converts the value returned by the underlying
 *  iterator's dereference operator to an rvalue reference.  Some
 *  generic algorithms can be called with move iterators to replace
 *  copying with moving.
 */
template<typename _Iterator>
  class move_iterator
  {
  protected:
    _Iterator _M_current;

    typedef iterator_traits<_Iterator>		__traits_type;

  public:
    typedef _Iterator					iterator_type;
    typedef typename __traits_type::iterator_category iterator_category;
    typedef typename __traits_type::value_type  	value_type;
    typedef typename __traits_type::difference_type	difference_type;
    // NB: DR 680.
    typedef _Iterator					pointer;
    typedef value_type&&				reference;

    move_iterator()
    : _M_current() { }

    explicit
    move_iterator(iterator_type __i)
    : _M_current(__i) { }

    template<typename _Iter>
move_iterator(const move_iterator<_Iter>& __i)
: _M_current(__i.base()) { }

    iterator_type
    base() const
    { return _M_current; }

    reference
    operator*() const
    { return std::move(*_M_current); }

    pointer
    operator->() const
    { return _M_current; }

    move_iterator&
    operator++()
    {
++_M_current;
return *this;
    }

    move_iterator
    operator++(int)
    {
move_iterator __tmp = *this;
++_M_current;
return __tmp;
    }

    move_iterator&
    operator--()
    {
--_M_current;
return *this;
    }

    move_iterator
    operator--(int)
    {
move_iterator __tmp = *this;
--_M_current;
return __tmp;
    }

    move_iterator
    operator+(difference_type __n) const
    { return move_iterator(_M_current + __n); }

    move_iterator&
    operator+=(difference_type __n)
    {
_M_current += __n;
return *this;
    }

    move_iterator
    operator-(difference_type __n) const
    { return move_iterator(_M_current - __n); }
  
    move_iterator&
    operator-=(difference_type __n)
    { 
_M_current -= __n;
return *this;
    }

    reference
    operator[](difference_type __n) const
    { return std::move(_M_current[__n]); }
  };

// Note: See __normal_iterator operators note from Gaby to understand
// why there are always 2 versions for most of the move_iterator
// operators.
template<typename _IteratorL, typename _IteratorR>
  inline bool
  operator==(const move_iterator<_IteratorL>& __x,
      const move_iterator<_IteratorR>& __y)
  { return __x.base() == __y.base(); }

template<typename _Iterator>
  inline bool
  operator==(const move_iterator<_Iterator>& __x,
      const move_iterator<_Iterator>& __y)
  { return __x.base() == __y.base(); }

template<typename _IteratorL, typename _IteratorR>
  inline bool
  operator!=(const move_iterator<_IteratorL>& __x,
      const move_iterator<_IteratorR>& __y)
  { return !(__x == __y); }

template<typename _Iterator>
  inline bool
  operator!=(const move_iterator<_Iterator>& __x,
      const move_iterator<_Iterator>& __y)
  { return !(__x == __y); }

template<typename _IteratorL, typename _IteratorR>
  inline bool
  operator<(const move_iterator<_IteratorL>& __x,
     const move_iterator<_IteratorR>& __y)
  { return __x.base() < __y.base(); }

template<typename _Iterator>
  inline bool
  operator<(const move_iterator<_Iterator>& __x,
     const move_iterator<_Iterator>& __y)
  { return __x.base() < __y.base(); }

template<typename _IteratorL, typename _IteratorR>
  inline bool
  operator<=(const move_iterator<_IteratorL>& __x,
      const move_iterator<_IteratorR>& __y)
  { return !(__y < __x); }

template<typename _Iterator>
  inline bool
  operator<=(const move_iterator<_Iterator>& __x,
      const move_iterator<_Iterator>& __y)
  { return !(__y < __x); }

template<typename _IteratorL, typename _IteratorR>
  inline bool
  operator>(const move_iterator<_IteratorL>& __x,
     const move_iterator<_IteratorR>& __y)
  { return __y < __x; }

template<typename _Iterator>
  inline bool
  operator>(const move_iterator<_Iterator>& __x,
     const move_iterator<_Iterator>& __y)
  { return __y < __x; }

template<typename _IteratorL, typename _IteratorR>
  inline bool
  operator>=(const move_iterator<_IteratorL>& __x,
      const move_iterator<_IteratorR>& __y)
  { return !(__x < __y); }

template<typename _Iterator>
  inline bool
  operator>=(const move_iterator<_Iterator>& __x,
      const move_iterator<_Iterator>& __y)
  { return !(__x < __y); }

// DR 685.
template<typename _IteratorL, typename _IteratorR>
  inline auto
  operator-(const move_iterator<_IteratorL>& __x,
     const move_iterator<_IteratorR>& __y)
  -> decltype(__x.base() - __y.base())
  { return __x.base() - __y.base(); }

template<typename _Iterator>
  inline auto
  operator-(const move_iterator<_Iterator>& __x,
     const move_iterator<_Iterator>& __y)
  -> decltype(__x.base() - __y.base())
  { return __x.base() - __y.base(); }

template<typename _Iterator>
  inline move_iterator<_Iterator>
  operator+(typename move_iterator<_Iterator>::difference_type __n,
     const move_iterator<_Iterator>& __x)
  { return __x + __n; }

template<typename _Iterator>
  inline move_iterator<_Iterator>
  make_move_iterator(_Iterator __i)
  { return move_iterator<_Iterator>(__i); }

template<typename _Iterator, typename _ReturnType
  = typename conditional<__move_if_noexcept_cond
    <typename iterator_traits<_Iterator>::value_type>::value,
              _Iterator, move_iterator<_Iterator>>::type>
  inline _ReturnType
  __make_move_if_noexcept_iterator(_Iterator __i)
  { return _ReturnType(__i); }

// @} group iterators

1137_GLIBCXX_END_NAMESPACE_VERSION
1138} // namespace

1140#define _GLIBCXX_MAKE_MOVE_ITERATOR(_Iter)std::make_move_iterator(_Iter) std::make_move_iterator(_Iter)
1141#define _GLIBCXX_MAKE_MOVE_IF_NOEXCEPT_ITERATOR(_Iter)std::__make_move_if_noexcept_iterator(_Iter) \
std::__make_move_if_noexcept_iterator(_Iter)
1143#else
1144#define _GLIBCXX_MAKE_MOVE_ITERATOR(_Iter)std::make_move_iterator(_Iter) (_Iter)
1145#define _GLIBCXX_MAKE_MOVE_IF_NOEXCEPT_ITERATOR(_Iter)std::__make_move_if_noexcept_iterator(_Iter) (_Iter)
1146#endif // C++11

1148#endif