/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 DEventWriterROOT.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/DEventWriterROOT.cc

libraries/ANALYSIS/DEventWriterROOT.cc

→

1#include "DEventWriterROOT.h"

3static bool BCAL_VERBOSE_OUTPUT = false;
4static bool FCAL_VERBOSE_OUTPUT = false;
5static bool CCAL_VERBOSE_OUTPUT = false;
6static bool DIRC_OUTPUT = true;

8static bool STORE_PULL_INFO = false;
9static bool STORE_ERROR_MATRIX_INFO = false;


12void DEventWriterROOT::Initialize(JEventLoop* locEventLoop)
13{
dInitNumThrownArraySize = 20;
dInitNumBeamArraySize = 20;
dInitNumTrackArraySize = 50;
dInitNumNeutralArraySize = 15;
dInitNumComboArraySize = 100;
dThrownTreeInterface = NULL__null;

locEventLoop->GetSingle(dAnalysisUtilities);

auto locReactions = DAnalysis::Get_Reactions(locEventLoop);

//CREATE & INITIALIZE ANALYSIS ACTIONS
for(size_t loc_i = 0; loc_i < locReactions.size(); ++loc_i)
{
if(!locReactions[loc_i]->Get_EnableTTreeOutputFlag())
	continue;

dCutActionMap_ThrownTopology[locReactions[loc_i]] = new DCutAction_ThrownTopology(locReactions[loc_i], true);
dCutActionMap_ThrownTopology[locReactions[loc_i]]->Initialize(locEventLoop);

dCutActionMap_TrueCombo[locReactions[loc_i]] = new DCutAction_TrueCombo(locReactions[loc_i], -1.0, true);
dCutActionMap_TrueCombo[locReactions[loc_i]]->Initialize(locEventLoop);

dCutActionMap_BDTSignalCombo[locReactions[loc_i]] = new DCutAction_BDTSignalCombo(locReactions[loc_i], 5.73303E-7, true, true); //+/- 5sigma
dCutActionMap_BDTSignalCombo[locReactions[loc_i]]->Initialize(locEventLoop);
}

//CREATE TREES
vector<const DMCThrown*> locMCThrowns;
locEventLoop->Get(locMCThrowns);

vector<const DReactionVertexInfo*> locVertexInfos;
locEventLoop->Get(locVertexInfos);

//Get Target Center Z
DApplication* locApplication = dynamic_cast<DApplication*>(locEventLoop->GetJApplication());
DGeometry* locGeometry = locApplication->GetDGeometry(locEventLoop->GetJEvent().GetRunNumber());
dTargetCenterZ = 65.0;
locGeometry->GetTargetZ(dTargetCenterZ);

//CREATE TTREES
for(auto& locVertexInfo : locVertexInfos)
{
auto locVertexReactions = locVertexInfo->Get_Reactions();
for(auto& locReaction : locVertexReactions)
{
	dVertexInfoMap.emplace(locReaction, locVertexInfo);
	if(locReaction->Get_EnableTTreeOutputFlag())
		Create_DataTree(locReaction, locEventLoop, !locMCThrowns.empty());
}
}
65}

67void DEventWriterROOT::Run_Update(JEventLoop* locEventLoop)
68{
locEventLoop->GetSingle(dAnalysisUtilities);
1
Calling 'JEventLoop::GetSingle'→

//Get Target Center Z
DApplication* locApplication = dynamic_cast<DApplication*>(locEventLoop->GetJApplication());
DGeometry* locGeometry = locApplication->GetDGeometry(locEventLoop->GetJEvent().GetRunNumber());
dTargetCenterZ = 65.0;
locGeometry->GetTargetZ(dTargetCenterZ);

// update run-dependent info/objects
for(auto& locMapPair : dCutActionMap_ThrownTopology)
locMapPair.second->Run_Update(locEventLoop);
for(auto& locMapPair : dCutActionMap_TrueCombo)
locMapPair.second->Run_Update(locEventLoop);
for(auto& locMapPair : dCutActionMap_BDTSignalCombo)
locMapPair.second->Run_Update(locEventLoop);
84}

86DEventWriterROOT::~DEventWriterROOT(void)
87{
//Delete tree interface objects
for(auto& locMapPair : dTreeInterfaceMap)
delete locMapPair.second;
if(dThrownTreeInterface != NULL__null)
delete dThrownTreeInterface;

for(auto& locMapPair : dTreeFillDataMap)
delete locMapPair.second;

//Delete actions
for(auto& locMapPair : dCutActionMap_ThrownTopology)
delete locMapPair.second;
for(auto& locMapPair : dCutActionMap_TrueCombo)
delete locMapPair.second;
for(auto& locMapPair : dCutActionMap_BDTSignalCombo)
delete locMapPair.second;
104}

106void DEventWriterROOT::Create_ThrownTree(JEventLoop* locEventLoop, string locOutputFileName) const
107{
if(dThrownTreeInterface != nullptr)
return; //Already setup for this thread!
dThrownTreeInterface = DTreeInterface::Create_DTreeInterface("Thrown_Tree", locOutputFileName); //set up this thread
if(dThrownTreeInterface->Get_BranchesCreatedFlag())
return; //branches already created: return

//TTREE BRANCHES
DTreeBranchRegister locBranchRegister;

//Get target PID
const DMCReaction* locMCReaction = NULL__null;
locEventLoop->GetSingle(locMCReaction);
Particle_t locTargetPID = locMCReaction->target.PID();

//setup target info
Create_UserTargetInfo(locBranchRegister, locTargetPID);

//create basic/misc. tree branches (run#, event#, etc.)
locBranchRegister.Register_Single<UInt_t>("RunNumber");
locBranchRegister.Register_Single<ULong64_t>("EventNumber");

//Thrown Data
Create_Branches_Thrown(locBranchRegister, true);

//CUSTOM
Create_CustomBranches_ThrownTree(locBranchRegister, locEventLoop);

//CREATE BRANCHES
dThrownTreeInterface->Create_Branches(locBranchRegister);
dThrownTreeInterface->Set_TreeIndexBranchNames("RunNumber", "EventNumber");
138}

140void DEventWriterROOT::Create_DataTree(const DReaction* locReaction, JEventLoop* locEventLoop, bool locIsMCDataFlag)
141{
string locReactionName = locReaction->Get_ReactionName();
string locOutputFileName = locReaction->Get_TTreeOutputFileName();
string locTreeName = locReactionName + string("_Tree");

//create fill object
dTreeFillDataMap[locReaction] = new DTreeFillData();

//create tree interface
DTreeInterface* locTreeInterface = DTreeInterface::Create_DTreeInterface(locTreeName, locOutputFileName);
dTreeInterfaceMap[locReaction] = locTreeInterface;
if(locTreeInterface->Get_BranchesCreatedFlag())
return; //branches already created, then return

//Branch register
DTreeBranchRegister locBranchRegister;

//fill maps
TMap* locPositionToNameMap = Create_UserInfoMaps(locBranchRegister, locEventLoop, locReaction);

161/******************************************************************** Create Branches ********************************************************************/

//create basic/misc. tree branches (run#, event#, etc.)
locBranchRegister.Register_Single<UInt_t>("RunNumber");
locBranchRegister.Register_Single<ULong64_t>("EventNumber");
locBranchRegister.Register_Single<UInt_t>("L1TriggerBits");

//create X4_Production
locBranchRegister.Register_Single<TLorentzVector>("X4_Production");

//create thrown branches
if(locIsMCDataFlag)
{
Create_Branches_Thrown(locBranchRegister, false);
locBranchRegister.Register_Single<Bool_t>("IsThrownTopology");
}

bool locBeamUsedFlag = DAnalysis::Get_IsFirstStepBeam(locReaction);

//create branches for final-state particle hypotheses
if(locBeamUsedFlag)
Create_Branches_Beam(locBranchRegister, locIsMCDataFlag);
Create_Branches_NeutralHypotheses(locBranchRegister, locIsMCDataFlag);
Create_Branches_ChargedHypotheses(locBranchRegister, locIsMCDataFlag);

//create branches for combos
locBranchRegister.Register_Single<UChar_t>("NumUnusedTracks");
Create_Branches_Combo(locBranchRegister, locReaction, locIsMCDataFlag, locPositionToNameMap);

//Kinematic fit data (pulls and covariance matrices)
Create_Branches_KinFitData(locBranchRegister, locEventLoop, locReaction, locIsMCDataFlag);

//Custom branches
Create_CustomBranches_DataTree(locBranchRegister, locEventLoop, locReaction, locIsMCDataFlag);

//Create branches
locTreeInterface->Create_Branches(locBranchRegister);
locTreeInterface->Set_TreeIndexBranchNames("RunNumber", "EventNumber");
199}

201TMap* DEventWriterROOT::Create_UserInfoMaps(DTreeBranchRegister& locBranchRegister, JEventLoop* locEventLoop, const DReaction* locReaction) const
202{
auto locReactionVertexInfo = dVertexInfoMap.find(locReaction)->second;

//kinfit type
DKinFitType locKinFitType = locReaction->Get_KinFitType();

//create & add reaction identification maps
TList* locUserInfo = locBranchRegister.Get_UserInfo();
TMap* locNameToPIDMap = new TMap();
locNameToPIDMap->SetName("NameToPIDMap");
locUserInfo->Add(locNameToPIDMap);

TMap* locNameToPositionMap = new TMap(); //not filled for target or initial particles
locNameToPositionMap->SetName("NameToPositionMap");
locUserInfo->Add(locNameToPositionMap);

TMap* locPositionToNameMap = new TMap(); //not filled for target or initial particles
locPositionToNameMap->SetName("PositionToNameMap");
locUserInfo->Add(locPositionToNameMap);

TMap* locPositionToPIDMap = new TMap();
locPositionToPIDMap->SetName("PositionToPIDMap");
locUserInfo->Add(locPositionToPIDMap);

TMap* locDecayProductMap = new TMap(); //excludes resonances!!! //excludes intermediate decays (e.g. if xi- -> pi-, lambda -> pi-, pi-, p: will be xi- -> pi-, pi-, p and no lambda decay present)
locDecayProductMap->SetName("DecayProductMap"); //parent name string -> tlist of decay product name strings
locUserInfo->Add(locDecayProductMap);

TMap* locMiscInfoMap = new TMap();
locMiscInfoMap->SetName("MiscInfoMap");
locUserInfo->Add(locMiscInfoMap);

TList* locParticleNameList = new TList();
locParticleNameList->SetName("ParticleNameList");
locUserInfo->Add(locParticleNameList);

//Set some misc info
ostringstream locKinFitTypeStream;
locKinFitTypeStream << locKinFitType;
locMiscInfoMap->Add(new TObjString("KinFitType"), new TObjString(locKinFitTypeStream.str().c_str()));

string ANALYSIS_VERSION_STRING = "";
if(gPARMS->Exists("ANALYSIS:DATAVERSIONSTRING"))
gPARMS->GetParameter("ANALYSIS:DATAVERSIONSTRING", ANALYSIS_VERSION_STRING);
if(ANALYSIS_VERSION_STRING != "")
locMiscInfoMap->Add(new TObjString("ANALYSIS:DATAVERSIONSTRING"), new TObjString(ANALYSIS_VERSION_STRING.c_str()));

string HDDM_DATA_VERSION_STRING = "";
if(gPARMS->Exists("REST:DATAVERSIONSTRING"))
gPARMS->GetParameter("REST:DATAVERSIONSTRING", HDDM_DATA_VERSION_STRING);
if(HDDM_DATA_VERSION_STRING != "")
locMiscInfoMap->Add(new TObjString("REST:DATAVERSIONSTRING"), new TObjString(HDDM_DATA_VERSION_STRING.c_str()));

string REST_JANA_CALIB_CONTEXT = "";
if(gPARMS->Exists("REST:JANACALIBCONTEXT"))
gPARMS->GetParameter("REST:JANACALIBCONTEXT", REST_JANA_CALIB_CONTEXT);
if(REST_JANA_CALIB_CONTEXT == "")
gPARMS->GetParameter("JANA_CALIB_CONTEXT", REST_JANA_CALIB_CONTEXT);
if(REST_JANA_CALIB_CONTEXT != "")
locMiscInfoMap->Add(new TObjString("REST:JANACALIBCONTEXT"), new TObjString(REST_JANA_CALIB_CONTEXT.c_str()));


if(gPARMS->Exists("ANALYSIS:BCAL_VERBOSE_ROOT_OUTPUT"))
gPARMS->GetParameter("ANALYSIS:BCAL_VERBOSE_ROOT_OUTPUT", BCAL_VERBOSE_OUTPUT);
if(gPARMS->Exists("ANALYSIS:FCAL_VERBOSE_ROOT_OUTPUT"))
gPARMS->GetParameter("ANALYSIS:FCAL_VERBOSE_ROOT_OUTPUT", FCAL_VERBOSE_OUTPUT);
if(gPARMS->Exists("ANALYSIS:CCAL_VERBOSE_ROOT_OUTPUT"))
gPARMS->GetParameter("ANALYSIS:CCAL_VERBOSE_ROOT_OUTPUT", CCAL_VERBOSE_OUTPUT);
if(gPARMS->Exists("ANALYSIS:DIRC_ROOT_OUTPUT"))
gPARMS->GetParameter("ANALYSIS:DIRC_ROOT_OUTPUT", DIRC_OUTPUT);
if(gPARMS->Exists("ANALYSIS:STORE_PULL_INFO"))
   gPARMS->GetParameter("ANALYSIS:STORE_PULL_INFO",STORE_PULL_INFO);
if(gPARMS->Exists("ANALYSIS:STORE_ERROR_MATRIX_INFO"))
  	gPARMS->GetParameter("ANALYSIS:STORE_ERROR_MATRIX_INFO",STORE_ERROR_MATRIX_INFO);

if(locKinFitType != d_NoFit)
{
DKinFitUtils_GlueX locKinFitUtils(locEventLoop);
size_t locNumConstraints = 0, locNumUnknowns = 0;
string locConstraintString = locKinFitUtils.Get_ConstraintInfo(locReactionVertexInfo, locReaction, locNumConstraints, locNumUnknowns);
locMiscInfoMap->Add(new TObjString("KinFitConstraints"), new TObjString(locConstraintString.c_str()));

ostringstream locKinFitInfoStream;
locKinFitInfoStream << locNumConstraints;
locMiscInfoMap->Add(new TObjString("NumKinFitConstraints"), new TObjString(locKinFitInfoStream.str().c_str()));

locKinFitInfoStream.str("");
locKinFitInfoStream << locNumUnknowns;
locMiscInfoMap->Add(new TObjString("NumKinFitUnknowns"), new TObjString(locKinFitInfoStream.str().c_str()));
}

//find the # particles of each pid
map<Particle_t, unsigned int> locParticleNumberMap;
map<Particle_t, unsigned int> locDecayingParticleNumberMap;
map<Particle_t, unsigned int> locTargetParticleNumberMap;
for(size_t loc_i = 0; loc_i < locReaction->Get_NumReactionSteps(); ++loc_i)
{
const DReactionStep* locReactionStep = locReaction->Get_ReactionStep(loc_i);
auto locFinalParticleIDs = locReactionStep->Get_FinalPIDs();

auto locTargetPID = locReactionStep->Get_TargetPID();
if(locTargetPID != Unknown)
{
	if(locTargetParticleNumberMap.find(locTargetPID) == locTargetParticleNumberMap.end())
		locTargetParticleNumberMap[locTargetPID] = 1;
	else
		++locTargetParticleNumberMap[locTargetPID];
}

for(size_t loc_j = 0; loc_j < locFinalParticleIDs.size(); ++loc_j)
{
	if(locReactionStep->Get_MissingParticleIndex() == int(loc_j)) //missing particle
		continue;
	Particle_t locPID = locFinalParticleIDs[loc_j];

	//see if decays in a future step
	int locDecayStepIndex = DAnalysis::Get_DecayStepIndex(locReaction, loc_i, loc_j);
	if(locDecayStepIndex >= 0) //decaying
	{
		if(locDecayingParticleNumberMap.find(locPID) == locDecayingParticleNumberMap.end())
			locDecayingParticleNumberMap[locPID] = 1;
		else
			++locDecayingParticleNumberMap[locPID];
	}
	else //detected, not decaying
	{
		if(locParticleNumberMap.find(locPID) == locParticleNumberMap.end())
			locParticleNumberMap[locPID] = 1;
		else
			++locParticleNumberMap[locPID];
	}
}
}

//Create map objects
map<Particle_t, unsigned int> locParticleNumberMap_Current, locDecayingParticleNumberMap_Current, locTargetParticleNumberMap_Current;
Particle_t locTargetPID = Unknown;
TObjString *locObjString_PID, *locObjString_Position, *locObjString_ParticleName;
map<int, string> locDecayingParticleNames; //key is step index where they decay
for(size_t loc_i = 0; loc_i < locReaction->Get_NumReactionSteps(); ++loc_i)
{
const DReactionStep* locReactionStep = locReaction->Get_ReactionStep(loc_i);

//initial particle
{
	ostringstream locPIDStream, locPositionStream, locParticleNameStream;
	Particle_t locPID = locReactionStep->Get_InitialPID();
	locPIDStream << PDGtype(locPID);
	locObjString_PID = new TObjString(locPIDStream.str().c_str());

	locPositionStream << loc_i << "_" << -1;
	locObjString_Position = new TObjString(locPositionStream.str().c_str());

	locPositionToPIDMap->Add(locObjString_Position, locObjString_PID);
	if((loc_i == 0) && ((locPID == Gamma) || (locPID == Electron) || (locPID == Positron)))
	{
		locParticleNameStream << "ComboBeam";
		locObjString_ParticleName = new TObjString(locParticleNameStream.str().c_str());
		locNameToPIDMap->Add(locObjString_ParticleName, locObjString_PID);
		locParticleNameList->AddLast(locObjString_ParticleName);
	}
	else //decaying particle
	{
		if(loc_i == 0)
			locParticleNameStream << "Decaying" << Convert_ToBranchName(ParticleType(locPID));
		else //name already created for final particle: use it
			locParticleNameStream << locDecayingParticleNames[loc_i];
		locObjString_ParticleName = new TObjString(locParticleNameStream.str().c_str());
		if(loc_i == 0) //in first step
		{
			locNameToPIDMap->Add(locObjString_ParticleName, locObjString_PID);
			locParticleNameList->AddLast(locObjString_ParticleName);
		}
	}
	locPositionToNameMap->Add(locObjString_Position, locObjString_ParticleName);
	locNameToPositionMap->Add(locObjString_ParticleName, locObjString_Position);
}

//target particle
Particle_t locTempTargetPID = locReactionStep->Get_TargetPID();
if(locTempTargetPID != Unknown)
{
	locTargetPID = locTempTargetPID;

	if(locTargetParticleNumberMap_Current.find(locTargetPID) == locTargetParticleNumberMap_Current.end())
		locTargetParticleNumberMap_Current[locTargetPID] = 1;
	else
		++locTargetParticleNumberMap_Current[locTargetPID];

	ostringstream locPIDStream, locPositionStream, locParticleNameStream;
	locPIDStream << PDGtype(locTargetPID);
	locObjString_PID = new TObjString(locPIDStream.str().c_str());

	locPositionStream << loc_i << "_" << -2;
	locObjString_Position = new TObjString(locPositionStream.str().c_str());

	locPositionToPIDMap->Add(locObjString_Position, locObjString_PID);

	locParticleNameStream << "Target";
	if(locTargetParticleNumberMap[locTargetPID] > 1)
		locParticleNameStream << locTargetParticleNumberMap_Current[locTargetPID];
	locObjString_ParticleName = new TObjString(locParticleNameStream.str().c_str());

	locNameToPositionMap->Add(locObjString_ParticleName, locObjString_Position);
	locNameToPIDMap->Add(locObjString_ParticleName, locObjString_PID);
	locPositionToNameMap->Add(locObjString_Position, locObjString_ParticleName);

	locParticleNameList->AddLast(locObjString_ParticleName);
}

//final particles
auto locFinalParticleIDs = locReactionStep->Get_FinalPIDs();
for(size_t loc_j = 0; loc_j < locFinalParticleIDs.size(); ++loc_j)
{
	ostringstream locPIDStream, locPositionStream;
	Particle_t locPID = locFinalParticleIDs[loc_j];
	locPIDStream << PDGtype(locPID);
	locObjString_PID = new TObjString(locPIDStream.str().c_str());

	locPositionStream << loc_i << "_" << loc_j;
	locObjString_Position = new TObjString(locPositionStream.str().c_str());

	if(locReactionStep->Get_MissingParticleIndex() == int(loc_j)) //missing particle
	{
		ostringstream locParticleNameStream;
		locParticleNameStream << "Missing" << Convert_ToBranchName(ParticleType(locPID));
		locObjString_ParticleName = new TObjString(locParticleNameStream.str().c_str());
		locNameToPositionMap->Add(locObjString_ParticleName, locObjString_Position);
		locPositionToNameMap->Add(locObjString_Position, locObjString_ParticleName);
		locNameToPIDMap->Add(locObjString_ParticleName, locObjString_PID);
		string locPIDName = locParticleNameStream.str() + string("__PID");
		locMiscInfoMap->Add(new TObjString(locPIDName.c_str()), locObjString_PID);
		ostringstream locMassStream;
		locMassStream << ParticleMass(locPID);
		string locMassName = locParticleNameStream.str() + string("__Mass");
		locMiscInfoMap->Add(new TObjString(locMassName.c_str()), new TObjString(locMassStream.str().c_str()));
		locParticleNameList->AddLast(locObjString_ParticleName);
		continue;
	}

	//build name
	ostringstream locParticleNameStream;
	//see if decays in a future step
	int locDecayStepIndex = DAnalysis::Get_DecayStepIndex(locReaction, loc_i, loc_j);
	if(locDecayStepIndex >= 0) //decays
	{
		if(locDecayingParticleNumberMap_Current.find(locPID) == locDecayingParticleNumberMap_Current.end())
			locDecayingParticleNumberMap_Current[locPID] = 1;
		else
			++locDecayingParticleNumberMap_Current[locPID];

		locParticleNameStream << "Decaying" << Convert_ToBranchName(ParticleType(locPID));
		if(locDecayingParticleNumberMap[locPID] > 1)
			locParticleNameStream << locDecayingParticleNumberMap_Current[locPID];
		locDecayingParticleNames[locDecayStepIndex] = locParticleNameStream.str();
	}
	else
	{
		if(locParticleNumberMap_Current.find(locPID) == locParticleNumberMap_Current.end())
			locParticleNumberMap_Current[locPID] = 1;
		else
			++locParticleNumberMap_Current[locPID];

		locParticleNameStream << Convert_ToBranchName(ParticleType(locPID));
		if(locParticleNumberMap[locPID] > 1)
			locParticleNameStream << locParticleNumberMap_Current[locPID];
	}

	locObjString_ParticleName = new TObjString(locParticleNameStream.str().c_str());
	locParticleNameList->AddLast(locObjString_ParticleName);

	locPositionToPIDMap->Add(locObjString_Position, locObjString_PID);
	locNameToPositionMap->Add(locObjString_ParticleName, locObjString_Position);
	locPositionToNameMap->Add(locObjString_Position, locObjString_ParticleName);
	locNameToPIDMap->Add(locObjString_ParticleName, locObjString_PID);
	if(locDecayStepIndex >= 0)
	{
		ostringstream locMassStream;
		locMassStream << ParticleMass(locPID);
		string locMassName = locParticleNameStream.str() + string("__Mass");
		locMiscInfoMap->Add(new TObjString(locMassName.c_str()), new TObjString(locMassStream.str().c_str()));
	}
}
}

//setup target info
Create_UserTargetInfo(locBranchRegister, locTargetPID);

//fill decay product map
deque<size_t> locSavedSteps;
for(size_t loc_i = 0; loc_i < locReaction->Get_NumReactionSteps(); ++loc_i)
{
const DReactionStep* locReactionStep = locReaction->Get_ReactionStep(loc_i);

//initial particle
Particle_t locPID = locReactionStep->Get_InitialPID();
if(loc_i == 0)
	continue;
if(!IsFixedMass(locPID))
	continue; //don't save resonance decays to the map

//check to see if this decay has already been saved
bool locStepAlreadySavedFlag = false;
for(size_t loc_j = 0; loc_j < locSavedSteps.size(); ++loc_j)
{
	if(locSavedSteps[loc_j] != loc_i)
		continue;
	locStepAlreadySavedFlag = true;
	break;
}
if(locStepAlreadySavedFlag)
	continue;

//construct the name 
ostringstream locParticleNameStream;
locParticleNameStream << locDecayingParticleNames[loc_i];
locObjString_ParticleName = new TObjString(locParticleNameStream.str().c_str());

TList* locDecayProductNames = NULL__null;
Get_DecayProductNames(locReaction, loc_i, locPositionToNameMap, locDecayProductNames, locSavedSteps);
locDecayProductMap->Add(locObjString_ParticleName, locDecayProductNames); //parent name string -> tobjarray of decay product name strings		
}

return locPositionToNameMap;
526}

528void DEventWriterROOT::Get_DecayProductNames(const DReaction* locReaction, size_t locReactionStepIndex, TMap* locPositionToNameMap, TList*& locDecayProductNames, deque<size_t>& locSavedSteps) const
529{
const DReactionStep* locReactionStep = locReaction->Get_ReactionStep(locReactionStepIndex);

if(locDecayProductNames == NULL__null)
locDecayProductNames = new TList();

auto locFinalParticleIDs = locReactionStep->Get_FinalPIDs();
for(size_t loc_j = 0; loc_j < locFinalParticleIDs.size(); ++loc_j)
{
//see if decays in a future step //save and continue if doesn't decay
int locDecayStepIndex = DAnalysis::Get_DecayStepIndex(locReaction, locReactionStepIndex, loc_j);
if(locDecayStepIndex < 0)
{
	ostringstream locPositionStream;
	locPositionStream << locReactionStepIndex << "_" << loc_j;
	locDecayProductNames->AddLast(locPositionToNameMap->GetValue(locPositionStream.str().c_str()));
	continue;
}

//add decay products
Get_DecayProductNames(locReaction, locDecayStepIndex, locPositionToNameMap, locDecayProductNames, locSavedSteps);
}

locSavedSteps.push_back(locReactionStepIndex);
553}

555void DEventWriterROOT::Create_UserTargetInfo(DTreeBranchRegister& locBranchRegister, Particle_t locTargetPID) const
556{
TList* locUserInfo = locBranchRegister.Get_UserInfo();
TMap* locMiscInfoMap = (TMap*)locUserInfo->FindObject("MiscInfoMap");
if(locMiscInfoMap == NULL__null)
{
locMiscInfoMap = new TMap();
locMiscInfoMap->SetName("MiscInfoMap");
locUserInfo->Add(locMiscInfoMap);
}

//PID
ostringstream locPIDStream;
locPIDStream << PDGtype(locTargetPID);
locMiscInfoMap->Add(new TObjString("Target__PID"), new TObjString(locPIDStream.str().c_str()));

//Mass
ostringstream locMassStream;
locMassStream << ParticleMass(locTargetPID);
locMiscInfoMap->Add(new TObjString("Target__Mass"), new TObjString(locMassStream.str().c_str()));

//X, Y
ostringstream locZeroStream;
locZeroStream << 0.0;
TObjString* locObjString_Zero = new TObjString(locZeroStream.str().c_str());
locMiscInfoMap->Add(new TObjString("Target__CenterX"), locObjString_Zero);
locMiscInfoMap->Add(new TObjString("Target__CenterY"), locObjString_Zero);

//Z
ostringstream locPositionStream;
locPositionStream << dTargetCenterZ;
TObjString* locObjString_Position = new TObjString(locPositionStream.str().c_str());
locMiscInfoMap->Add(new TObjString("Target__CenterZ"), locObjString_Position);
588}

590void DEventWriterROOT::Create_Branches_Thrown(DTreeBranchRegister& locBranchRegister, bool locIsOnlyThrownFlag) const
591{
//BEAM
locBranchRegister.Register_Single<Int_t>(Build_BranchName("ThrownBeam", "PID"));
locBranchRegister.Register_Single<TLorentzVector>(Build_BranchName("ThrownBeam", "X4")); //reported at target center
locBranchRegister.Register_Single<TLorentzVector>(Build_BranchName("ThrownBeam", "P4"));
locBranchRegister.Register_Single<Float_t>(Build_BranchName("ThrownBeam", "GeneratedEnergy"));

//EVENT-WIDE INFO
locBranchRegister.Register_Single<ULong64_t>("NumPIDThrown_FinalState"); //19 digits
locBranchRegister.Register_Single<ULong64_t>("PIDThrown_Decaying"); //19 digits
locBranchRegister.Register_Single<Float_t>("MCWeight");

//PRODUCTS
Create_Branches_ThrownParticles(locBranchRegister, locIsOnlyThrownFlag);
605}

607void DEventWriterROOT::Create_Branches_ThrownParticles(DTreeBranchRegister& locBranchRegister, bool locIsOnlyThrownFlag) const
608{
string locParticleBranchName = "Thrown";

string locArraySizeString = "NumThrown";
locBranchRegister.Register_Single<UInt_t>(locArraySizeString);

//IDENTIFIERS
locBranchRegister.Register_FundamentalArray<Int_t>(Build_BranchName(locParticleBranchName, "ParentIndex"), locArraySizeString, dInitNumThrownArraySize);
locBranchRegister.Register_FundamentalArray<Int_t>(Build_BranchName(locParticleBranchName, "PID"), locArraySizeString, dInitNumThrownArraySize);
if(!locIsOnlyThrownFlag)
{
locBranchRegister.Register_FundamentalArray<Int_t>(Build_BranchName(locParticleBranchName, "MatchID"), locArraySizeString, dInitNumThrownArraySize);
locBranchRegister.Register_FundamentalArray<Float_t>(Build_BranchName(locParticleBranchName, "MatchFOM"), locArraySizeString, dInitNumThrownArraySize);
}

//KINEMATICS: THROWN //at the production vertex
locBranchRegister.Register_ClonesArray<TLorentzVector>(Build_BranchName(locParticleBranchName, "X4"), dInitNumThrownArraySize);
locBranchRegister.Register_ClonesArray<TLorentzVector>(Build_BranchName(locParticleBranchName, "P4"), dInitNumThrownArraySize);
626}

628void DEventWriterROOT::Create_Branches_Beam(DTreeBranchRegister& locBranchRegister, bool locIsMCDataFlag) const
629{
string locArraySizeString = "NumBeam";
locBranchRegister.Register_Single<UInt_t>(locArraySizeString);

string locParticleBranchName = "Beam";

//IDENTIFIER
locBranchRegister.Register_FundamentalArray<Int_t>(Build_BranchName(locParticleBranchName, "PID"), locArraySizeString, dInitNumBeamArraySize);
if(locIsMCDataFlag)
locBranchRegister.Register_FundamentalArray<Bool_t>(Build_BranchName(locParticleBranchName, "IsGenerator"), locArraySizeString, dInitNumBeamArraySize);

//KINEMATICS: MEASURED //at the production vertex
locBranchRegister.Register_ClonesArray<TLorentzVector>(Build_BranchName(locParticleBranchName, "X4_Measured"), dInitNumThrownArraySize);
locBranchRegister.Register_ClonesArray<TLorentzVector>(Build_BranchName(locParticleBranchName, "P4_Measured"), dInitNumThrownArraySize);
643}

645void DEventWriterROOT::Create_Branches_ChargedHypotheses(DTreeBranchRegister& locBranchRegister, bool locIsMCDataFlag) const
646{
string locArraySizeString = "NumChargedHypos";
locBranchRegister.Register_Single<UInt_t>(locArraySizeString);

string locParticleBranchName = "ChargedHypo";

//IDENTIFIERS / MATCHING
locBranchRegister.Register_FundamentalArray<Int_t>(Build_BranchName(locParticleBranchName, "TrackID"), locArraySizeString, dInitNumTrackArraySize);
locBranchRegister.Register_FundamentalArray<Int_t>(Build_BranchName(locParticleBranchName, "PID"), locArraySizeString, dInitNumTrackArraySize);
if(locIsMCDataFlag)
locBranchRegister.Register_FundamentalArray<Int_t>(Build_BranchName(locParticleBranchName, "ThrownIndex"), locArraySizeString, dInitNumTrackArraySize);

//KINEMATICS //at the production vertex
locBranchRegister.Register_ClonesArray<TLorentzVector>(Build_BranchName(locParticleBranchName, "X4_Measured"), dInitNumTrackArraySize);
locBranchRegister.Register_ClonesArray<TLorentzVector>(Build_BranchName(locParticleBranchName, "P4_Measured"), dInitNumTrackArraySize);

// Global PID
locBranchRegister.Register_FundamentalArray<Float_t>(Build_BranchName(locParticleBranchName, "PIDFOM"), locArraySizeString, dInitNumTrackArraySize);

//TRACKING INFO
locBranchRegister.Register_FundamentalArray<UInt_t>(Build_BranchName(locParticleBranchName, "NDF_Tracking"), locArraySizeString, dInitNumTrackArraySize);
locBranchRegister.Register_FundamentalArray<Float_t>(Build_BranchName(locParticleBranchName, "ChiSq_Tracking"), locArraySizeString, dInitNumTrackArraySize);
locBranchRegister.Register_FundamentalArray<UInt_t>(Build_BranchName(locParticleBranchName, "NDF_DCdEdx"), locArraySizeString, dInitNumTrackArraySize);
locBranchRegister.Register_FundamentalArray<Float_t>(Build_BranchName(locParticleBranchName, "ChiSq_DCdEdx"), locArraySizeString, dInitNumTrackArraySize);
locBranchRegister.Register_FundamentalArray<Float_t>(Build_BranchName(locParticleBranchName, "dEdx_CDC"), locArraySizeString, dInitNumTrackArraySize);
locBranchRegister.Register_FundamentalArray<Float_t>(Build_BranchName(locParticleBranchName, "dEdx_CDC_integral"), locArraySizeString, dInitNumTrackArraySize);
locBranchRegister.Register_FundamentalArray<Float_t>(Build_BranchName(locParticleBranchName, "dEdx_FDC"), locArraySizeString, dInitNumTrackArraySize);

//TIMING INFO
locBranchRegister.Register_FundamentalArray<Float_t>(Build_BranchName(locParticleBranchName, "HitTime"), locArraySizeString, dInitNumTrackArraySize);
locBranchRegister.Register_FundamentalArray<Float_t>(Build_BranchName(locParticleBranchName, "RFDeltaTVar"), locArraySizeString, dInitNumTrackArraySize);

//PID QUALITY
locBranchRegister.Register_FundamentalArray<Float_t>(Build_BranchName(locParticleBranchName, "Beta_Timing"), locArraySizeString, dInitNumTrackArraySize);
locBranchRegister.Register_FundamentalArray<Float_t>(Build_BranchName(locParticleBranchName, "ChiSq_Timing"), locArraySizeString, dInitNumTrackArraySize);
locBranchRegister.Register_FundamentalArray<UInt_t>(Build_BranchName(locParticleBranchName, "NDF_Timing"), locArraySizeString, dInitNumTrackArraySize);

//HIT ENERGY
locBranchRegister.Register_FundamentalArray<Float_t>(Build_BranchName(locParticleBranchName, "dEdx_TOF"), locArraySizeString, dInitNumTrackArraySize);
locBranchRegister.Register_FundamentalArray<Float_t>(Build_BranchName(locParticleBranchName, "dEdx_ST"), locArraySizeString, dInitNumTrackArraySize);
locBranchRegister.Register_FundamentalArray<Float_t>(Build_BranchName(locParticleBranchName, "Energy_BCAL"), locArraySizeString, dInitNumNeutralArraySize);
locBranchRegister.Register_FundamentalArray<Float_t>(Build_BranchName(locParticleBranchName, "Energy_BCALPreshower"), locArraySizeString, dInitNumNeutralArraySize);
if(BCAL_VERBOSE_OUTPUT) {
locBranchRegister.Register_FundamentalArray<Float_t>(Build_BranchName(locParticleBranchName, "Energy_BCALLayer2"), locArraySizeString, dInitNumNeutralArraySize);
locBranchRegister.Register_FundamentalArray<Float_t>(Build_BranchName(locParticleBranchName, "Energy_BCALLayer3"), locArraySizeString, dInitNumNeutralArraySize);
locBranchRegister.Register_FundamentalArray<Float_t>(Build_BranchName(locParticleBranchName, "Energy_BCALLayer4"), locArraySizeString, dInitNumNeutralArraySize);
}
locBranchRegister.Register_FundamentalArray<Float_t>(Build_BranchName(locParticleBranchName, "SigLong_BCAL"), locArraySizeString, dInitNumNeutralArraySize);
locBranchRegister.Register_FundamentalArray<Float_t>(Build_BranchName(locParticleBranchName, "SigTheta_BCAL"), locArraySizeString, dInitNumNeutralArraySize);
locBranchRegister.Register_FundamentalArray<Float_t>(Build_BranchName(locParticleBranchName, "SigTrans_BCAL"), locArraySizeString, dInitNumNeutralArraySize);
locBranchRegister.Register_FundamentalArray<Float_t>(Build_BranchName(locParticleBranchName, "RMSTime_BCAL"), locArraySizeString, dInitNumNeutralArraySize);
locBranchRegister.Register_FundamentalArray<Float_t>(Build_BranchName(locParticleBranchName, "Energy_FCAL"), locArraySizeString, dInitNumNeutralArraySize);
locBranchRegister.Register_FundamentalArray<Float_t>(Build_BranchName(locParticleBranchName, "E1E9_FCAL"), locArraySizeString, dInitNumNeutralArraySize);
locBranchRegister.Register_FundamentalArray<Float_t>(Build_BranchName(locParticleBranchName, "E9E25_FCAL"), locArraySizeString, dInitNumNeutralArraySize);
locBranchRegister.Register_FundamentalArray<Float_t>(Build_BranchName(locParticleBranchName, "SumU_FCAL"), locArraySizeString, dInitNumNeutralArraySize);
locBranchRegister.Register_FundamentalArray<Float_t>(Build_BranchName(locParticleBranchName, "SumV_FCAL"), locArraySizeString, dInitNumNeutralArraySize);

locBranchRegister.Register_FundamentalArray<Float_t>(Build_BranchName(locParticleBranchName, "Energy_CCAL"), locArraySizeString, dInitNumNeutralArraySize);


//SHOWER MATCHING:
locBranchRegister.Register_FundamentalArray<Float_t>(Build_BranchName(locParticleBranchName, "TrackBCAL_DeltaPhi"), locArraySizeString, dInitNumTrackArraySize);
locBranchRegister.Register_FundamentalArray<Float_t>(Build_BranchName(locParticleBranchName, "TrackBCAL_DeltaZ"), locArraySizeString, dInitNumTrackArraySize);
locBranchRegister.Register_FundamentalArray<Float_t>(Build_BranchName(locParticleBranchName, "TrackFCAL_DOCA"), locArraySizeString, dInitNumTrackArraySize);

//DIRC:
if(DIRC_OUTPUT) {
locBranchRegister.Register_FundamentalArray<Int_t>(Build_BranchName(locParticleBranchName, "NumPhotons_DIRC"), locArraySizeString, dInitNumTrackArraySize);
locBranchRegister.Register_FundamentalArray<Float_t>(Build_BranchName(locParticleBranchName, "ExtrapolatedX_DIRC"), locArraySizeString, dInitNumTrackArraySize);
locBranchRegister.Register_FundamentalArray<Float_t>(Build_BranchName(locParticleBranchName, "ExtrapolatedY_DIRC"), locArraySizeString, dInitNumTrackArraySize);
locBranchRegister.Register_FundamentalArray<Float_t>(Build_BranchName(locParticleBranchName, "ThetaC_DIRC"), locArraySizeString, dInitNumTrackArraySize);
locBranchRegister.Register_FundamentalArray<Float_t>(Build_BranchName(locParticleBranchName, "Lele_DIRC"), locArraySizeString, dInitNumTrackArraySize);
locBranchRegister.Register_FundamentalArray<Float_t>(Build_BranchName(locParticleBranchName, "Lpi_DIRC"), locArraySizeString, dInitNumTrackArraySize);
locBranchRegister.Register_FundamentalArray<Float_t>(Build_BranchName(locParticleBranchName, "Lk_DIRC"), locArraySizeString, dInitNumTrackArraySize);
locBranchRegister.Register_FundamentalArray<Float_t>(Build_BranchName(locParticleBranchName, "Lp_DIRC"), locArraySizeString, dInitNumTrackArraySize);
}
722}

724void DEventWriterROOT::Create_Branches_NeutralHypotheses(DTreeBranchRegister& locBranchRegister, bool locIsMCDataFlag) const
725{
string locArraySizeString = "NumNeutralHypos";
string locParticleBranchName = "NeutralHypo";
locBranchRegister.Register_Single<UInt_t>(locArraySizeString);

//IDENTIFIERS / MATCHING
locBranchRegister.Register_FundamentalArray<Int_t>(Build_BranchName(locParticleBranchName, "NeutralID"), locArraySizeString, dInitNumNeutralArraySize);
locBranchRegister.Register_FundamentalArray<Int_t>(Build_BranchName(locParticleBranchName, "PID"), locArraySizeString, dInitNumNeutralArraySize);
if(locIsMCDataFlag)
locBranchRegister.Register_FundamentalArray<Int_t>(Build_BranchName(locParticleBranchName, "ThrownIndex"), locArraySizeString, dInitNumNeutralArraySize);

//KINEMATICS //is combo-dependent: P4 defined by X4, X4 defined by other tracks
locBranchRegister.Register_ClonesArray<TLorentzVector>(Build_BranchName(locParticleBranchName, "X4_Measured"), dInitNumNeutralArraySize);
locBranchRegister.Register_ClonesArray<TLorentzVector>(Build_BranchName(locParticleBranchName, "P4_Measured"), dInitNumNeutralArraySize);

//PID QUALITY
locBranchRegister.Register_FundamentalArray<Float_t>(Build_BranchName(locParticleBranchName, "Beta_Timing"), locArraySizeString, dInitNumNeutralArraySize);
locBranchRegister.Register_FundamentalArray<Float_t>(Build_BranchName(locParticleBranchName, "ChiSq_Timing"), locArraySizeString, dInitNumNeutralArraySize);
locBranchRegister.Register_FundamentalArray<UInt_t>(Build_BranchName(locParticleBranchName, "NDF_Timing"), locArraySizeString, dInitNumNeutralArraySize);

//SHOWER INFO
locBranchRegister.Register_ClonesArray<TLorentzVector>(Build_BranchName(locParticleBranchName, "X4_Shower"), dInitNumNeutralArraySize);
locBranchRegister.Register_FundamentalArray<Float_t>(Build_BranchName(locParticleBranchName, "ShowerQuality"), locArraySizeString, dInitNumNeutralArraySize);
locBranchRegister.Register_FundamentalArray<Float_t>(Build_BranchName(locParticleBranchName, "Energy_BCAL"), locArraySizeString, dInitNumNeutralArraySize);
locBranchRegister.Register_FundamentalArray<Float_t>(Build_BranchName(locParticleBranchName, "Energy_BCALPreshower"), locArraySizeString, dInitNumNeutralArraySize);
if(BCAL_VERBOSE_OUTPUT) {
locBranchRegister.Register_FundamentalArray<Float_t>(Build_BranchName(locParticleBranchName, "Energy_BCALLayer2"), locArraySizeString, dInitNumNeutralArraySize);
locBranchRegister.Register_FundamentalArray<Float_t>(Build_BranchName(locParticleBranchName, "Energy_BCALLayer3"), locArraySizeString, dInitNumNeutralArraySize);
locBranchRegister.Register_FundamentalArray<Float_t>(Build_BranchName(locParticleBranchName, "Energy_BCALLayer4"), locArraySizeString, dInitNumNeutralArraySize);
}
locBranchRegister.Register_FundamentalArray<Float_t>(Build_BranchName(locParticleBranchName, "SigLong_BCAL"), locArraySizeString, dInitNumNeutralArraySize);
locBranchRegister.Register_FundamentalArray<Float_t>(Build_BranchName(locParticleBranchName, "SigTheta_BCAL"), locArraySizeString, dInitNumNeutralArraySize);
locBranchRegister.Register_FundamentalArray<Float_t>(Build_BranchName(locParticleBranchName, "SigTrans_BCAL"), locArraySizeString, dInitNumNeutralArraySize);
locBranchRegister.Register_FundamentalArray<Float_t>(Build_BranchName(locParticleBranchName, "RMSTime_BCAL"), locArraySizeString, dInitNumNeutralArraySize);
locBranchRegister.Register_FundamentalArray<Float_t>(Build_BranchName(locParticleBranchName, "Energy_FCAL"), locArraySizeString, dInitNumNeutralArraySize);
if(FCAL_VERBOSE_OUTPUT) {
locBranchRegister.Register_FundamentalArray<Float_t>(Build_BranchName(locParticleBranchName, "E1E9_FCAL"), locArraySizeString, dInitNumNeutralArraySize);
locBranchRegister.Register_FundamentalArray<Float_t>(Build_BranchName(locParticleBranchName, "E9E25_FCAL"), locArraySizeString, dInitNumNeutralArraySize);
locBranchRegister.Register_FundamentalArray<Float_t>(Build_BranchName(locParticleBranchName, "SumU_FCAL"), locArraySizeString, dInitNumNeutralArraySize);
locBranchRegister.Register_FundamentalArray<Float_t>(Build_BranchName(locParticleBranchName, "SumV_FCAL"), locArraySizeString, dInitNumNeutralArraySize);
}
locBranchRegister.Register_FundamentalArray<Float_t>(Build_BranchName(locParticleBranchName, "Energy_CCAL"), locArraySizeString, dInitNumNeutralArraySize);
//NEARBY TRACKS
locBranchRegister.Register_FundamentalArray<Float_t>(Build_BranchName(locParticleBranchName, "TrackBCAL_DeltaPhi"), locArraySizeString, dInitNumNeutralArraySize);
locBranchRegister.Register_FundamentalArray<Float_t>(Build_BranchName(locParticleBranchName, "TrackBCAL_DeltaZ"), locArraySizeString, dInitNumNeutralArraySize);
locBranchRegister.Register_FundamentalArray<Float_t>(Build_BranchName(locParticleBranchName, "TrackFCAL_DOCA"), locArraySizeString, dInitNumNeutralArraySize);

//PHOTON PID INFO
//Computed using DVertex (best estimate of reaction vertex using all "good" tracks)
//Can be used to compute timing chisq //is invalid for non-photons, so computed assuming photon PID
//Variance of X4_Measured.T() - RFTime, regardless of which RF bunch is chosen. //RF bunch is combo-depende
locBranchRegister.Register_FundamentalArray<Float_t>(Build_BranchName(locParticleBranchName, "PhotonRFDeltaTVar"), locArraySizeString, dInitNumNeutralArraySize);
777}

779void DEventWriterROOT::Create_Branches_Combo(DTreeBranchRegister& locBranchRegister, const DReaction* locReaction, bool locIsMCDataFlag, TMap* locPositionToNameMap) const
780{
auto locReactionVertexInfo = dVertexInfoMap.find(locReaction)->second;
string locNumComboString = "NumCombos";
locBranchRegister.Register_Single<UInt_t>(locNumComboString);

//kinfit type
DKinFitType locKinFitType = locReaction->Get_KinFitType();
bool locKinFitFlag = (locKinFitType != d_NoFit);
bool locVertexKinFitFlag = locKinFitFlag && (locKinFitType != d_P4Fit);

//Is-cut
locBranchRegister.Register_FundamentalArray<Bool_t>("IsComboCut", locNumComboString, dInitNumComboArraySize);

//create combo-dependent, particle-independent branches
if(locIsMCDataFlag)
{
locBranchRegister.Register_FundamentalArray<Bool_t>("IsTrueCombo", locNumComboString, dInitNumComboArraySize);
locBranchRegister.Register_FundamentalArray<Bool_t>("IsBDTSignalCombo", locNumComboString, dInitNumComboArraySize);
}

locBranchRegister.Register_FundamentalArray<Float_t>("RFTime_Measured", locNumComboString, dInitNumComboArraySize);
if(locKinFitFlag)
{
locBranchRegister.Register_FundamentalArray<Float_t>("ChiSq_KinFit", locNumComboString, dInitNumComboArraySize);
locBranchRegister.Register_FundamentalArray<UInt_t>("NDF_KinFit", locNumComboString, dInitNumComboArraySize);
if((locKinFitType == d_SpacetimeFit) || (locKinFitType == d_P4AndSpacetimeFit))
	locBranchRegister.Register_FundamentalArray<Float_t>("RFTime_KinFit", locNumComboString, dInitNumComboArraySize);
}
locBranchRegister.Register_FundamentalArray<UChar_t>("NumUnusedShowers", locNumComboString, dInitNumComboArraySize);
locBranchRegister.Register_FundamentalArray<Float_t>("Energy_UnusedShowers", locNumComboString, dInitNumComboArraySize);
locBranchRegister.Register_FundamentalArray<UChar_t>("NumUnusedShowers_Quality", locNumComboString, dInitNumComboArraySize);
locBranchRegister.Register_FundamentalArray<Float_t>("Energy_UnusedShowers_Quality", locNumComboString, dInitNumComboArraySize);
locBranchRegister.Register_FundamentalArray<Float_t>("SumPMag_UnusedTracks", locNumComboString, dInitNumComboArraySize);
locBranchRegister.Register_ClonesArray<TVector3>("SumP3_UnusedTracks", dInitNumComboArraySize);

map<Particle_t, unsigned int> locParticleNumberMap_Current;
for(size_t loc_i = 0; loc_i < locReaction->Get_NumReactionSteps(); ++loc_i)
{
const DReactionStep* locReactionStep = locReaction->Get_ReactionStep(loc_i);

//initial particle
Particle_t locInitialPID = locReactionStep->Get_InitialPID();
//should check to make sure the beam particle isn't missing...
if((loc_i == 0) && (locReactionStep->Get_InitialPID() != Unknown))
	Create_Branches_BeamComboParticle(locBranchRegister, locInitialPID, locKinFitType);
else //decaying
{
	//get the branch name
	ostringstream locPositionStream;
	locPositionStream << loc_i << "_-1";
	TObjString* locObjString = (TObjString*)locPositionToNameMap->GetValue(locPositionStream.str().c_str());
	string locParticleBranchName = (const char*)(locObjString->GetString());

	if(IsFixedMass(locInitialPID) && locReactionStep->Get_KinFitConstrainInitMassFlag() && ((locKinFitType == d_P4Fit) || (locKinFitType == d_P4AndVertexFit) || (locKinFitType == d_P4AndSpacetimeFit)))
		locBranchRegister.Register_ClonesArray<TLorentzVector>(Build_BranchName(locParticleBranchName, "P4_KinFit"), dInitNumComboArraySize);
	if((loc_i == 0) || IsDetachedVertex(locInitialPID))
		locBranchRegister.Register_ClonesArray<TLorentzVector>(Build_BranchName(locParticleBranchName, "X4"), dInitNumComboArraySize);

	auto locStepVertexInfo = locReactionVertexInfo->Get_StepVertexInfo(loc_i);
	auto locParentVertexInfo = locStepVertexInfo->Get_ParentVertexInfo();
	if(IsDetachedVertex(locInitialPID) && locVertexKinFitFlag && (locParentVertexInfo != nullptr) && locStepVertexInfo->Get_FittableVertexFlag() && locParentVertexInfo->Get_FittableVertexFlag())
		locBranchRegister.Register_FundamentalArray<Float_t>(Build_BranchName(locParticleBranchName, "PathLengthSigma"), locNumComboString, dInitNumComboArraySize);
}

//final particles
auto locFinalParticleIDs = locReactionStep->Get_FinalPIDs();
for(size_t loc_j = 0; loc_j < locFinalParticleIDs.size(); ++loc_j)
{
	int locDecayStepIndex = DAnalysis::Get_DecayStepIndex(locReaction, loc_i, loc_j);
	if(locDecayStepIndex >= 0)
		continue; //decaying particle

	//get the branch name
	ostringstream locPositionStream;
	locPositionStream << loc_i << "_" << loc_j;
	TObjString* locObjString = (TObjString*)locPositionToNameMap->GetValue(locPositionStream.str().c_str());
	string locParticleBranchName = (const char*)(locObjString->GetString());

	//missing particle
	if(locReactionStep->Get_MissingParticleIndex() == int(loc_j))
	{
		// missing particle
		if((locKinFitType == d_P4Fit) || (locKinFitType == d_P4AndVertexFit) || (locKinFitType == d_P4AndSpacetimeFit))
			locBranchRegister.Register_ClonesArray<TLorentzVector>(Build_BranchName(locParticleBranchName, "P4_KinFit"), dInitNumComboArraySize);
		continue;
	}

	Particle_t locPID = locFinalParticleIDs[loc_j];
	if(ParticleCharge(locPID) == 0)
		Create_Branches_ComboNeutral(locBranchRegister, locParticleBranchName, locKinFitType);
	else
		Create_Branches_ComboTrack(locBranchRegister, locParticleBranchName, locKinFitType);
}
}
874}

876void DEventWriterROOT::Create_Branches_BeamComboParticle(DTreeBranchRegister& locBranchRegister, Particle_t locBeamPID, DKinFitType locKinFitType) const
877{
string locParticleBranchName = "ComboBeam";
string locArraySizeString = "NumCombos";

//IDENTIFIER
locBranchRegister.Register_FundamentalArray<Int_t>(Build_BranchName(locParticleBranchName, "BeamIndex"), locArraySizeString, dInitNumComboArraySize);

//KINEMATICS: KINFIT //at the interaction vertex
if(locKinFitType != d_NoFit)
{
if(((locKinFitType != d_VertexFit) && (locKinFitType != d_SpacetimeFit)) || (ParticleCharge(locBeamPID) != 0))
	locBranchRegister.Register_ClonesArray<TLorentzVector>(Build_BranchName(locParticleBranchName, "P4_KinFit"), dInitNumComboArraySize);
if(locKinFitType != d_P4Fit)
	locBranchRegister.Register_ClonesArray<TLorentzVector>(Build_BranchName(locParticleBranchName, "X4_KinFit"), dInitNumComboArraySize);
}
892}

894void DEventWriterROOT::Create_Branches_ComboTrack(DTreeBranchRegister& locBranchRegister, string locParticleBranchName, DKinFitType locKinFitType) const
895{
string locArraySizeString = "NumCombos";

//IDENTIFIER
locBranchRegister.Register_FundamentalArray<Int_t>(Build_BranchName(locParticleBranchName, "ChargedIndex"), locArraySizeString, dInitNumComboArraySize);

//MEASURED PID
locBranchRegister.Register_FundamentalArray<Float_t>(Build_BranchName(locParticleBranchName, "Beta_Timing_Measured"), locArraySizeString, dInitNumComboArraySize);
locBranchRegister.Register_FundamentalArray<Float_t>(Build_BranchName(locParticleBranchName, "ChiSq_Timing_Measured"), locArraySizeString, dInitNumComboArraySize);

//KINFIT PID
if((locKinFitType != d_NoFit) && (locKinFitType != d_SpacetimeFit) && (locKinFitType != d_P4AndSpacetimeFit))
{
locBranchRegister.Register_FundamentalArray<Float_t>(Build_BranchName(locParticleBranchName, "Beta_Timing_KinFit"), locArraySizeString, dInitNumComboArraySize);
locBranchRegister.Register_FundamentalArray<Float_t>(Build_BranchName(locParticleBranchName, "ChiSq_Timing_KinFit"), locArraySizeString, dInitNumComboArraySize);
}

//KINFIT INFO //at the production vertex
if(locKinFitType != d_NoFit)
{
//update p4 even if vertex-only fit, because charged momentum propagated through b-field
locBranchRegister.Register_ClonesArray<TLorentzVector>(Build_BranchName(locParticleBranchName, "P4_KinFit"), dInitNumComboArraySize);
if(locKinFitType != d_P4Fit)
	locBranchRegister.Register_ClonesArray<TLorentzVector>(Build_BranchName(locParticleBranchName, "X4_KinFit"), dInitNumComboArraySize);
}
920}

922void DEventWriterROOT::Create_Branches_ComboNeutral(DTreeBranchRegister& locBranchRegister, string locParticleBranchName, DKinFitType locKinFitType) const
923{
string locArraySizeString = "NumCombos";

//IDENTIFIER
locBranchRegister.Register_FundamentalArray<Int_t>(Build_BranchName(locParticleBranchName, "NeutralIndex"), locArraySizeString, dInitNumComboArraySize);

//KINEMATICS //is combo-dependent: P4 defined by X4, X4 defined by other tracks
locBranchRegister.Register_ClonesArray<TLorentzVector>(Build_BranchName(locParticleBranchName, "X4_Measured"), dInitNumComboArraySize);
locBranchRegister.Register_ClonesArray<TLorentzVector>(Build_BranchName(locParticleBranchName, "P4_Measured"), dInitNumComboArraySize);

//MEASURED PID
locBranchRegister.Register_FundamentalArray<Float_t>(Build_BranchName(locParticleBranchName, "Beta_Timing_Measured"), locArraySizeString, dInitNumComboArraySize);
if(locParticleBranchName.substr(0, 6) == "Photon")
locBranchRegister.Register_FundamentalArray<Float_t>(Build_BranchName(locParticleBranchName, "ChiSq_Timing_Measured"), locArraySizeString, dInitNumComboArraySize);

//KINFIT PID
if((locKinFitType != d_NoFit) && (locKinFitType != d_SpacetimeFit) && (locKinFitType != d_P4AndSpacetimeFit))
{
locBranchRegister.Register_FundamentalArray<Float_t>(Build_BranchName(locParticleBranchName, "Beta_Timing_KinFit"), locArraySizeString, dInitNumComboArraySize);
if(locParticleBranchName.substr(0, 6) == "Photon")
	locBranchRegister.Register_FundamentalArray<Float_t>(Build_BranchName(locParticleBranchName, "ChiSq_Timing_KinFit"), locArraySizeString, dInitNumComboArraySize);
}

//KINFIT INFO //at the production vertex
if(locKinFitType != d_NoFit)
{
locBranchRegister.Register_ClonesArray<TLorentzVector>(Build_BranchName(locParticleBranchName, "P4_KinFit"), dInitNumComboArraySize);
if(locKinFitType != d_P4Fit)
	locBranchRegister.Register_ClonesArray<TLorentzVector>(Build_BranchName(locParticleBranchName, "X4_KinFit"), dInitNumComboArraySize);
}
953}

955void DEventWriterROOT::Create_Branches_KinFitData(DTreeBranchRegister& locBranchRegister, JEventLoop* locEventLoop, const DReaction* locReaction, bool locIsMCDataFlag) const
956{

if(!STORE_PULL_INFO && !STORE_ERROR_MATRIX_INFO)
return;
  
  //Get the fit-type (P4, P4 + Vertex, and Vertex-only)
  DKinFitType kfitType = locReaction->Get_KinFitType();
  
  vector<Particle_t> finalchargedPIDs,finalneutralPIDs;
  size_t nReactionSteps = locReaction->Get_NumReactionSteps();
  int nFinalChargedParticles,nFinalNeutralParticles;
  
  vector<string> pidVector = getPIDinReaction(locReaction);
  map<Particle_t,int> assignMap = getNameForParticle(locReaction,pidVector);
   
  //Check, if at least two charged particles are present in the final state: (important for vertex-fit)
  int sumAllChargedParticles = 0;
  int sumAllNeutralParticles = 0;
  //First, check how many neutral/charged particles are in the final stage
  //The pull-information is avaiablable for certain particle combinations (charged + neutral) only:
  //----------------------------------------------------------------------------------------------
  for(size_t loc_i=0;loc_i<nReactionSteps;loc_i++){
      
      //Charged Particles:
      finalchargedPIDs = locReaction->Get_ReactionStep(loc_i)->Get_FinalPIDs(true,d_Charged,true);
      nFinalChargedParticles = finalchargedPIDs.size();
      sumAllChargedParticles += nFinalChargedParticles;
      
       //Neutral Particles:
      finalneutralPIDs = locReaction->Get_ReactionStep(loc_i)->Get_FinalPIDs(true,d_Neutral,true);
      nFinalNeutralParticles = finalneutralPIDs.size();
      sumAllNeutralParticles += nFinalNeutralParticles;
  }
  //----------------------------------------------------------------------------------------------
  
  if(sumAllChargedParticles == 1 && (kfitType == d_VertexFit || kfitType == d_P4AndVertexFit)){
      STORE_PULL_INFO = false;
      jout <<"  "<< endl;
      jout <<">>> WARNING: Only one charged track in the final state! No vertex fit possible! No pulls will be stored... <<<"<< endl;
      jout <<"  "<< endl;
  }
  
  if(sumAllNeutralParticles > 0 && kfitType == d_VertexFit){
      STORE_PULL_INFO = false;
      jout <<"  "<< endl;
      jout <<">>> WARNING: Neutral particles involved in pure vertex fit! No pulls will be stored... <<<"<< endl;
      jout <<"  "<< endl;
  }
  
  if(kfitType == d_NoFit){
STORE_PULL_INFO = false;
      jout <<"  "<< endl;
      jout <<">>> WARNING: No fit specified! No pulls will be stored... <<<"<< endl;
      jout <<"  "<< endl;
}

  setPullFlag(locReaction,STORE_PULL_INFO);
      
  //Get Pulls for the beam:
  //Particle_t beamPID = locReaction->Get_ReactionStep(0)->Get_InitialPID();
  
  string particleCovM = "numEntries_ParticleErrM";
  string showerCovM = "numEntries_ShowerErrM";
  string decayCovM = "numEntries_DecayErrM";
  if(STORE_ERROR_MATRIX_INFO){
     locBranchRegister.Register_Single<Int_t>(particleCovM);
     locBranchRegister.Register_Single<Int_t>(showerCovM); 
  }
  
  //Set the branches for the beam photon:
  if(STORE_PULL_INFO) 
  	setTreePullBranches(locBranchRegister,"ComboBeam",kfitType,dInitNumComboArraySize,true);
  if(STORE_ERROR_MATRIX_INFO) 
  	locBranchRegister.Register_FundamentalArray< Float_t >(Build_BranchName("ComboBeam","ErrMatrix"),particleCovM,nEntriesParticleCov);

  //----------------------------------------------------------------------------------------------
  for(size_t loc_i=0;loc_i<nReactionSteps;loc_i++){
      
      //Charged Particles:
      finalchargedPIDs = locReaction->Get_ReactionStep(loc_i)->Get_FinalPIDs(true,d_Charged,true);
      nFinalChargedParticles = finalchargedPIDs.size();
      //----------------------------------------------------------------------------------------------
      for(int loc_j=0;loc_j<nFinalChargedParticles;loc_j++){
          //Define branches for final charged particles:
          string branchName = assignName(finalchargedPIDs.at(loc_j),assignMap);
          assignMap[finalchargedPIDs.at(loc_j)]--;
          if(branchName != "nada"){
                if(STORE_PULL_INFO) setTreePullBranches(locBranchRegister,branchName,kfitType,dInitNumComboArraySize,false);
                if(STORE_ERROR_MATRIX_INFO) locBranchRegister.Register_FundamentalArray< Float_t >(Build_BranchName(branchName,"ErrMatrix"),particleCovM,nEntriesParticleCov);
          }
      }
      //----------------------------------------------------------------------------------------------
      
      //Neutral Particles:
      finalneutralPIDs = locReaction->Get_ReactionStep(loc_i)->Get_FinalPIDs(true,d_Neutral,true);
      nFinalNeutralParticles = finalneutralPIDs.size();
      //----------------------------------------------------------------------------------------------
      for(int loc_j=0;loc_j<nFinalNeutralParticles;loc_j++){
          //Define branches for final neutral particles:
          string branchName = assignName(finalneutralPIDs.at(loc_j),assignMap);
          assignMap[finalneutralPIDs.at(loc_j)]--;
          if(branchName != "nada"){
              if(STORE_PULL_INFO) setTreePullBranches(locBranchRegister,branchName,kfitType,dInitNumComboArraySize,true);
              if(finalneutralPIDs.at(loc_j) == Gamma && STORE_ERROR_MATRIX_INFO) locBranchRegister.Register_FundamentalArray< Float_t >(Build_BranchName(branchName,"ErrMatrix"),showerCovM,nEntriesShowerCov);
          } 
       }
      //----------------------------------------------------------------------------------------------
   }
  //----------------------------------------------------------------------------------------------
  
  assignMap.clear();
  abundanceMap.clear();
1068}

1070void DEventWriterROOT::Fill_ThrownTree(JEventLoop* locEventLoop) const
1071{
vector<const DMCThrown*> locMCThrowns_FinalState;
locEventLoop->Get(locMCThrowns_FinalState, "FinalState");

vector<const DMCThrown*> locMCThrowns_Decaying;
locEventLoop->Get(locMCThrowns_Decaying, "Decaying");

const DMCReaction* locMCReaction = NULL__null;
locEventLoop->GetSingle(locMCReaction);

ULong64_t locNumPIDThrown_FinalState = 0, locPIDThrown_Decaying = 0;
Compute_ThrownPIDInfo(locMCThrowns_FinalState, locMCThrowns_Decaying, locNumPIDThrown_FinalState, locPIDThrown_Decaying);

vector<const DMCThrown*> locMCThrownsToSave;
map<const DMCThrown*, unsigned int> locThrownIndexMap;
Group_ThrownParticles(locMCThrowns_FinalState, locMCThrowns_Decaying, locMCThrownsToSave, locThrownIndexMap);

vector<const DBeamPhoton*> locTaggedMCGenBeams;
locEventLoop->Get(locTaggedMCGenBeams, "TAGGEDMCGEN");

vector<const DBeamPhoton*> locMCGenBeams;
locEventLoop->Get(locMCGenBeams, "MCGEN");

const DBeamPhoton* locTaggedMCGenBeam = locTaggedMCGenBeams.empty() ? locMCGenBeams[0] : locTaggedMCGenBeams[0]; //if empty: will have to do. 

japp->RootWriteLock();

//primary event info
dThrownTreeFillData.Fill_Single<UInt_t>("RunNumber", locEventLoop->GetJEvent().GetRunNumber());
dThrownTreeFillData.Fill_Single<ULong64_t>("EventNumber", locEventLoop->GetJEvent().GetEventNumber());

//throwns
Fill_ThrownInfo(&dThrownTreeFillData, locMCReaction, locTaggedMCGenBeam, locMCThrownsToSave, locThrownIndexMap, locNumPIDThrown_FinalState, locPIDThrown_Decaying);

//Custom Branches
Fill_CustomBranches_ThrownTree(&dThrownTreeFillData, locEventLoop, locMCReaction, locMCThrownsToSave);

//FILL TTREE
dThrownTreeInterface->Fill(dThrownTreeFillData);


japp->RootUnLock();

1114}

1116void DEventWriterROOT::Fill_DataTrees(JEventLoop* locEventLoop, string locDReactionTag) const
1117{
if(locDReactionTag == "Thrown")
{
cout << "WARNING: CANNOT FILL THROWN TREE WITH THIS FUNCTION." << endl;
return;
}

vector<const DAnalysisResults*> locAnalysisResultsVector;
locEventLoop->Get(locAnalysisResultsVector);

vector<const DReaction*> locReactionsWithTag;
locEventLoop->Get(locReactionsWithTag, locDReactionTag.c_str());

for(size_t loc_i = 0; loc_i < locAnalysisResultsVector.size(); ++loc_i)
{
deque<const DParticleCombo*> locPassedParticleCombos;
locAnalysisResultsVector[loc_i]->Get_PassedParticleCombos(locPassedParticleCombos);
if(locPassedParticleCombos.empty())
	continue;

const DReaction* locReaction = locAnalysisResultsVector[loc_i]->Get_Reaction();
if(!locReaction->Get_EnableTTreeOutputFlag())
	continue;

bool locReactionFoundFlag = false;
for(size_t loc_j = 0; loc_j < locReactionsWithTag.size(); ++loc_j)
{
	if(locReactionsWithTag[loc_j] != locReaction)
		continue;
	locReactionFoundFlag = true;
	break;
}
if(!locReactionFoundFlag)
	continue; //reaction not from this factory, continue

Fill_DataTree(locEventLoop, locReaction, locPassedParticleCombos);
}
1154}

1156void DEventWriterROOT::Fill_DataTree(JEventLoop* locEventLoop, const DReaction* locReaction, deque<const DParticleCombo*>& locParticleCombos) const
1157{
if(locReaction->Get_ReactionName() == "Thrown")
{
cout << "WARNING: CANNOT FILL THROWN TREE WITH THIS FUNCTION." << endl;
return;
}

if(!locReaction->Get_EnableTTreeOutputFlag())
{
cout << "WARNING: ROOT TTREE OUTPUT NOT ENABLED FOR THIS DREACTION (" << locReaction->Get_ReactionName() << ")" << endl;
return;
}

/***************************************************** GET THROWN INFO *****************************************************/

vector<const DMCThrown*> locMCThrowns_FinalState;
locEventLoop->Get(locMCThrowns_FinalState, "FinalState");

vector<const DMCThrown*> locMCThrowns_Decaying;
locEventLoop->Get(locMCThrowns_Decaying, "Decaying");

vector<const DMCThrownMatching*> locMCThrownMatchingVector;
locEventLoop->Get(locMCThrownMatchingVector);
const DMCThrownMatching* locMCThrownMatching = locMCThrownMatchingVector.empty() ? NULL__null : locMCThrownMatchingVector[0];

vector<const DMCReaction*> locMCReactions;
locEventLoop->Get(locMCReactions);
const DMCReaction* locMCReaction = locMCReactions.empty() ? NULL__null : locMCReactions[0];

vector<const DBeamPhoton*> locTaggedMCGenBeams;
locEventLoop->Get(locTaggedMCGenBeams, "TAGGEDMCGEN");

vector<const DBeamPhoton*> locMCGenBeams;
locEventLoop->Get(locMCGenBeams, "MCGEN");

 const DBeamPhoton* locTaggedMCGenBeam = nullptr;

if (locTaggedMCGenBeams.empty()){
    if ( !locMCGenBeams.empty() ) locTaggedMCGenBeam = locMCGenBeams[0];
 }
 else locTaggedMCGenBeam = locTaggedMCGenBeams[0];

//Pre-compute thrown info
ULong64_t locNumPIDThrown_FinalState = 0, locPIDThrown_Decaying = 0;
Compute_ThrownPIDInfo(locMCThrowns_FinalState, locMCThrowns_Decaying, locNumPIDThrown_FinalState, locPIDThrown_Decaying);

//Pre-compute thrown info
vector<const DMCThrown*> locMCThrownsToSave;
map<const DMCThrown*, unsigned int> locThrownIndexMap;
Group_ThrownParticles(locMCThrowns_FinalState, locMCThrowns_Decaying, locMCThrownsToSave, locThrownIndexMap);

/****************************************************** GET PARTICLES ******************************************************/

bool locSaveUnusedFlag = locReaction->Get_SaveUnusedFlag();

//Get PIDs need for reaction
auto locDetectedPIDs = locReaction->Get_FinalPIDs(-1, false, false, d_AllCharges, false);
set<Particle_t> locReactionPIDs;
for(size_t loc_j = 0; loc_j < locDetectedPIDs.size(); ++loc_j)
locReactionPIDs.insert(locDetectedPIDs[loc_j]);

//GET HYPOTHESES
vector<const DChargedTrackHypothesis*> locChargedTrackHypotheses;
vector<const DNeutralParticleHypothesis*> locNeutralParticleHypotheses;
if(locSaveUnusedFlag)
{
locChargedTrackHypotheses = Get_ChargedHypotheses(locEventLoop);
locNeutralParticleHypotheses = Get_NeutralHypotheses(locEventLoop, locReactionPIDs);
}
else
{
locChargedTrackHypotheses = Get_ChargedHypotheses_Used(locEventLoop, locReaction, locParticleCombos);
locNeutralParticleHypotheses = Get_NeutralHypotheses_Used(locEventLoop, locReaction, locReactionPIDs, locParticleCombos);
}

//GET BEAM PHOTONS
bool locBeamUsedFlag = DAnalysis::Get_IsFirstStepBeam(locReaction);
vector<const DBeamPhoton*> locBeamPhotons = Get_BeamPhotons(locParticleCombos);

//create map of particles to array index:
//used for pointing combo particles to the appropriate measured-particle array index
//for hypos, they are the preselect versions if they exist, else the combo versions (e.g. PID not in REST)

//indices: beam
map<pair<oid_t, Particle_t>, size_t> locObjectToArrayIndexMap; //particle_t necessary for neutral showers!
for(size_t loc_i = 0; loc_i < locBeamPhotons.size(); ++loc_i)
{
pair<oid_t, Particle_t> locBeamPair(locBeamPhotons[loc_i]->id, locBeamPhotons[loc_i]->PID());
locObjectToArrayIndexMap[locBeamPair] = loc_i;
}

//indices: charged
for(size_t loc_i = 0; loc_i < locChargedTrackHypotheses.size(); ++loc_i)
{
const DTrackTimeBased* locTrackTimeBased = locChargedTrackHypotheses[loc_i]->Get_TrackTimeBased();
pair<oid_t, Particle_t> locTrackPair(locTrackTimeBased->id, locTrackTimeBased->PID());
locObjectToArrayIndexMap[locTrackPair] = loc_i;
}

//indices: neutral
for(size_t loc_i = 0; loc_i < locNeutralParticleHypotheses.size(); ++loc_i)
{
const DNeutralShower* locNeutralShower = locNeutralParticleHypotheses[loc_i]->Get_NeutralShower();
pair<oid_t, Particle_t> locShowerPair(locNeutralShower->id, locNeutralParticleHypotheses[loc_i]->PID());
locObjectToArrayIndexMap[locShowerPair] = loc_i;
}

/**************************************************** GET MISCELLANEOUS ****************************************************/

//GET DETECTOR MATCHES
const DDetectorMatches* locDetectorMatches = NULL__null;
locEventLoop->GetSingle(locDetectorMatches);

//GET DVERTEX
const DVertex* locVertex = NULL__null;
locEventLoop->GetSingle(locVertex);

//GET TRIGGER
const DTrigger* locTrigger = NULL__null;
locEventLoop->GetSingle(locTrigger);

/************************************************* EXECUTE ANALYSIS ACTIONS ************************************************/
      
japp->RootWriteLock();

Bool_t locIsThrownTopologyFlag = kFALSE;
vector<Bool_t> locIsTrueComboFlags;
vector<Bool_t> locIsBDTSignalComboFlags;
if(locMCReaction != NULL__null)
{
DCutAction_ThrownTopology* locThrownTopologyAction = dCutActionMap_ThrownTopology.find(locReaction)->second;
locIsThrownTopologyFlag = (*locThrownTopologyAction)(locEventLoop, NULL__null); //combo not used/needed
for(size_t loc_i = 0; loc_i < locParticleCombos.size(); ++loc_i)
{
	DCutAction_TrueCombo* locTrueComboAction = dCutActionMap_TrueCombo.find(locReaction)->second;
	locIsTrueComboFlags.push_back((*locTrueComboAction)(locEventLoop, locParticleCombos[loc_i]));

	DCutAction_BDTSignalCombo* locBDTSignalComboAction = dCutActionMap_BDTSignalCombo.find(locReaction)->second;
	locIsBDTSignalComboFlags.push_back((*locBDTSignalComboAction)(locEventLoop, locParticleCombos[loc_i]));
}
}

/***************************************************** FILL TTREE DATA *****************************************************/

//Get tree fill data
DTreeFillData* locTreeFillData = dTreeFillDataMap.find(locReaction)->second;

//PRIMARY EVENT INFO
locTreeFillData->Fill_Single<UInt_t>("RunNumber", locEventLoop->GetJEvent().GetRunNumber());
locTreeFillData->Fill_Single<ULong64_t>("EventNumber", locEventLoop->GetJEvent().GetEventNumber());
locTreeFillData->Fill_Single<UInt_t>("L1TriggerBits", locTrigger->Get_L1TriggerBits());

//PRODUCTION X4
DLorentzVector locProductionX4 = locVertex->dSpacetimeVertex;
TLorentzVector locProductionTX4(locProductionX4.X(), locProductionX4.Y(), locProductionX4.Z(), locProductionX4.T());
locTreeFillData->Fill_Single<TLorentzVector>("X4_Production", locProductionTX4);

//THROWN INFORMATION
if(locMCReaction != NULL__null)
{
Fill_ThrownInfo(locTreeFillData, locMCReaction, locTaggedMCGenBeam, locMCThrownsToSave, locThrownIndexMap, locNumPIDThrown_FinalState, locPIDThrown_Decaying, locMCThrownMatching);
locTreeFillData->Fill_Single<Bool_t>("IsThrownTopology", locIsThrownTopologyFlag);
}

//INDEPENDENT BEAM PARTICLES
if(locBeamUsedFlag)
{
//however, only fill with beam particles that are in the combos
locTreeFillData->Fill_Single<UInt_t>("NumBeam", UInt_t(locBeamPhotons.size()));
for(size_t loc_i = 0; loc_i < locBeamPhotons.size(); ++loc_i)
	Fill_BeamData(locTreeFillData, loc_i, locBeamPhotons[loc_i], locVertex, locMCThrownMatching);
}

//INDEPENDENT CHARGED TRACKS
locTreeFillData->Fill_Single<UInt_t>("NumChargedHypos", UInt_t(locChargedTrackHypotheses.size()));
for(size_t loc_i = 0; loc_i < locChargedTrackHypotheses.size(); ++loc_i)
Fill_ChargedHypo(locTreeFillData, loc_i, locChargedTrackHypotheses[loc_i], locMCThrownMatching, locThrownIndexMap, locDetectorMatches);

//INDEPENDENT NEUTRAL PARTICLES
locTreeFillData->Fill_Single<UInt_t>("NumNeutralHypos", UInt_t(locNeutralParticleHypotheses.size()));
for(size_t loc_i = 0; loc_i < locNeutralParticleHypotheses.size(); ++loc_i)
Fill_NeutralHypo(locTreeFillData, loc_i, locNeutralParticleHypotheses[loc_i], locMCThrownMatching, locThrownIndexMap, locDetectorMatches);

//UNUSED TRACKS
double locSumPMag_UnusedTracks = 0.0;
TVector3 locSumP3_UnusedTracks;
int locNumUnusedTracks = dAnalysisUtilities->Calc_Momentum_UnusedTracks(locEventLoop, locParticleCombos[0], locSumPMag_UnusedTracks, locSumP3_UnusedTracks);
locTreeFillData->Fill_Single<UChar_t>("NumUnusedTracks", locNumUnusedTracks);

//COMBOS
locTreeFillData->Fill_Single<UInt_t>("NumCombos", UInt_t(locParticleCombos.size()));
for(size_t loc_i = 0; loc_i < locParticleCombos.size(); ++loc_i)
{
Fill_ComboData(locTreeFillData, locReaction, locParticleCombos[loc_i], loc_i, locObjectToArrayIndexMap);

//ENERGY OF UNUSED SHOWERS (access to event loop required)
double locEnergy_UnusedShowers = 0.;
double locEnergy_UnusedShowers_Quality = 0.;
int locNumber_UnusedShowers_Quality = 0;
int locNumber_UnusedShowers = dAnalysisUtilities->Calc_Energy_UnusedShowers(locEventLoop, locParticleCombos[loc_i], locEnergy_UnusedShowers, locNumber_UnusedShowers_Quality, locEnergy_UnusedShowers_Quality);
locTreeFillData->Fill_Array<UChar_t>("NumUnusedShowers", locNumber_UnusedShowers, loc_i);
locTreeFillData->Fill_Array<Float_t>("Energy_UnusedShowers", locEnergy_UnusedShowers, loc_i);
locTreeFillData->Fill_Array<UChar_t>("NumUnusedShowers_Quality", locNumber_UnusedShowers_Quality, loc_i);
locTreeFillData->Fill_Array<Float_t>("Energy_UnusedShowers_Quality", locEnergy_UnusedShowers_Quality, loc_i);

//MOMENTUM OF UNUSED TRACKS (access to event loop required)
double locSumPMag_UnusedTracks = 0;
TVector3 locSumP3_UnusedTracks;
dAnalysisUtilities->Calc_Momentum_UnusedTracks(locEventLoop, locParticleCombos[loc_i], locSumPMag_UnusedTracks, locSumP3_UnusedTracks);
locTreeFillData->Fill_Array<Float_t>("SumPMag_UnusedTracks", locSumPMag_UnusedTracks, loc_i);
locTreeFillData->Fill_Array<TVector3>("SumP3_UnusedTracks", locSumP3_UnusedTracks, loc_i);

if(locMCReaction != NULL__null)
{
	locTreeFillData->Fill_Array<Bool_t>("IsTrueCombo", locIsTrueComboFlags[loc_i], loc_i);
	locTreeFillData->Fill_Array<Bool_t>("IsBDTSignalCombo", locIsTrueComboFlags[loc_i], loc_i);
}
}

//Kinematic fit data (pulls and covariance matrices)
Fill_KinFitData(locTreeFillData, locEventLoop, locReaction, locMCReaction, locMCThrownsToSave, locMCThrownMatching, locDetectorMatches, locBeamPhotons, locChargedTrackHypotheses, locNeutralParticleHypotheses, locParticleCombos);

//CUSTOM
Fill_CustomBranches_DataTree(locTreeFillData, locEventLoop, locReaction, locMCReaction, locMCThrownsToSave, locMCThrownMatching, locDetectorMatches, locBeamPhotons, locChargedTrackHypotheses, locNeutralParticleHypotheses, locParticleCombos);

//FILL
DTreeInterface* locTreeInterface = dTreeInterfaceMap.find(locReaction)->second;
locTreeInterface->Fill(*locTreeFillData);	

japp->RootUnLock();

1388}

1390vector<const DBeamPhoton*> DEventWriterROOT::Get_BeamPhotons(const deque<const DParticleCombo*>& locParticleCombos) const
1391{
//however, only fill with beam particles that are in the combos
set<const DBeamPhoton*> locBeamPhotonSet;
vector<const DBeamPhoton*> locBeamPhotons;
for(size_t loc_j = 0; loc_j < locParticleCombos.size(); ++loc_j)
{
const DParticleComboStep* locParticleComboStep = locParticleCombos[loc_j]->Get_ParticleComboStep(0);
const DKinematicData* locKinematicData = locParticleComboStep->Get_InitialParticle_Measured();
if(locKinematicData == NULL__null)
	continue;
const DBeamPhoton* locBeamPhoton = dynamic_cast<const DBeamPhoton*>(locKinematicData);
if(locBeamPhoton == NULL__null)
	continue;
if(locBeamPhotonSet.find(locBeamPhoton) != locBeamPhotonSet.end())
	continue; //already included
locBeamPhotonSet.insert(locBeamPhoton);
locBeamPhotons.push_back(locBeamPhoton);
}

return locBeamPhotons;
1411}

1413vector<const DChargedTrackHypothesis*> DEventWriterROOT::Get_ChargedHypotheses(JEventLoop* locEventLoop) const
1414{
//For default/preselect, save all
//For combo, of new PIDs ONLY, save one of each for each track
//save the one with the same RF bunch as the common bunch

vector<const DChargedTrack*> locChargedTracks;
locEventLoop->Get(locChargedTracks, "Combo");

vector<const DChargedTrackHypothesis*> locChargedHyposToSave;
for(auto& locChargedTrack : locChargedTracks)
locChargedHyposToSave.insert(locChargedHyposToSave.end(), locChargedTrack->dChargedTrackHypotheses.begin(), locChargedTrack->dChargedTrackHypotheses.end());

return locChargedHyposToSave;
1427}

1429vector<const DChargedTrackHypothesis*> DEventWriterROOT::Get_ChargedHypotheses_Used(JEventLoop* locEventLoop, const DReaction* locReaction, const deque<const DParticleCombo*>& locParticleCombos) const
1430{
//get all hypos
vector<const DChargedTrackHypothesis*> locAllHypos = Get_ChargedHypotheses(locEventLoop);

//get used time-based tracks
set<const DTrackTimeBased*> locUsedTimeBasedTracks;
for(auto& locCombo : locParticleCombos)
{
auto locChargedParticles = locCombo->Get_FinalParticles_Measured(locReaction, d_Charged);
for(auto& locParticle : locChargedParticles)
	locUsedTimeBasedTracks.insert(static_cast<const DChargedTrackHypothesis*>(locParticle)->Get_TrackTimeBased());
}

//loop through "all" hypos, removing those that weren't used
for(auto locIterator = locAllHypos.begin(); locIterator != locAllHypos.end();)
{
const DTrackTimeBased* locTrackTimeBased = (*locIterator)->Get_TrackTimeBased();
if(locUsedTimeBasedTracks.find(locTrackTimeBased) != locUsedTimeBasedTracks.end())
	++locIterator;
else
   locIterator = locAllHypos.erase(locIterator); 
}

return locAllHypos;
1454}

1456vector<const DNeutralParticleHypothesis*> DEventWriterROOT::Get_NeutralHypotheses(JEventLoop* locEventLoop, const set<Particle_t>& locReactionPIDs) const
1457{
//For default/preselect, save all
//For combo, of new PIDs ONLY, save one of each for each shower
//save the one with the same RF bunch as the common bunch

vector<const DNeutralParticle*> locNeutralParticles;
locEventLoop->Get(locNeutralParticles, "Combo");

vector<const DNeutralParticleHypothesis*> locNeutralHyposToSave;
for(auto& locNeutralParticle : locNeutralParticles)
locNeutralHyposToSave.insert(locNeutralHyposToSave.end(), locNeutralParticle->dNeutralParticleHypotheses.begin(), locNeutralParticle->dNeutralParticleHypotheses.end());

return locNeutralHyposToSave;
1470}

1472vector<const DNeutralParticleHypothesis*> DEventWriterROOT::Get_NeutralHypotheses_Used(JEventLoop* locEventLoop, const DReaction* locReaction, const set<Particle_t>& locReactionPIDs, const deque<const DParticleCombo*>& locParticleCombos) const
1473{
//get all hypos
vector<const DNeutralParticleHypothesis*> locAllHypos = Get_NeutralHypotheses(locEventLoop, locReactionPIDs);

//get used showers
set<pair<const DNeutralShower*, Particle_t> > locUsedNeutralShowers;
for(auto& locCombo : locParticleCombos)
{
auto locNeutralParticles = locCombo->Get_FinalParticles_Measured(locReaction, d_Neutral);
for(auto& locParticle : locNeutralParticles)
{
	const DNeutralShower* locNeutralShower = static_cast<const DNeutralParticleHypothesis*>(locParticle)->Get_NeutralShower();
	pair<const DNeutralShower*, Particle_t> locShowerPair(locNeutralShower, locParticle->PID());
	locUsedNeutralShowers.insert(locShowerPair);
}
}

//loop through "all" hypos, removing those that weren't used
for(auto locIterator = locAllHypos.begin(); locIterator != locAllHypos.end();)
{
const DNeutralShower* locNeutralShower = (*locIterator)->Get_NeutralShower();
pair<const DNeutralShower*, Particle_t> locShowerPair(locNeutralShower, (*locIterator)->PID());
if(locUsedNeutralShowers.find(locShowerPair) == locUsedNeutralShowers.end())
   locIterator = locAllHypos.erase(locIterator);
else
	++locIterator;
}

return locAllHypos;
1502}

1504ULong64_t DEventWriterROOT::Calc_ParticleMultiplexID(Particle_t locPID) const
1505{
int locPower = ParticleMultiplexPower(locPID);
if(locPower == -1)
return 0;

int locIsFinalStateInt = Is_FinalStateParticle(locPID);
if(locPID == Pi0)
locIsFinalStateInt = 1;

if(locIsFinalStateInt == 1) //decimal
{
ULong64_t locParticleMultiplexID = 1;
for(int loc_i = 0; loc_i < locPower; ++loc_i)
	locParticleMultiplexID *= ULong64_t(10);
return locParticleMultiplexID;
}
//decaying: binary
return (ULong64_t(1) << ULong64_t(locPower)); //bit-shift
1523}

1525void DEventWriterROOT::Compute_ThrownPIDInfo(const vector<const DMCThrown*>& locMCThrowns_FinalState, const vector<const DMCThrown*>& locMCThrowns_Decaying, ULong64_t& locNumPIDThrown_FinalState, ULong64_t& locPIDThrown_Decaying) const
1526{
//THROWN PARTICLES BY PID
locNumPIDThrown_FinalState = 0;
for(size_t loc_i = 0; loc_i < locMCThrowns_FinalState.size(); ++loc_i) //final state
{
Particle_t locPID = locMCThrowns_FinalState[loc_i]->PID();
ULong64_t locPIDMultiplexID = Calc_ParticleMultiplexID(locPID);
if(locPIDMultiplexID == 0)
	continue; //unrecognized PID!!!
unsigned int locCurrentNumParticles = (locNumPIDThrown_FinalState / locPIDMultiplexID) % ULong64_t(10);
if(locCurrentNumParticles != 9)
	locNumPIDThrown_FinalState += locPIDMultiplexID;
}

locPIDThrown_Decaying = 0;
for(size_t loc_i = 0; loc_i < locMCThrowns_Decaying.size(); ++loc_i) //decaying
{
Particle_t locPID = locMCThrowns_Decaying[loc_i]->PID();
ULong64_t locPIDMultiplexID = Calc_ParticleMultiplexID(locPID);
if(locPIDMultiplexID == 0)
	continue; //unrecognized PID!!!
if(locPID != Pi0)
	locPIDThrown_Decaying |= locPIDMultiplexID; //bit-wise or
else //save pi0's as final state instead of decaying
{
	unsigned int locCurrentNumParticles = (locNumPIDThrown_FinalState / locPIDMultiplexID) % ULong64_t(10);
	if(locCurrentNumParticles != 9)
		locNumPIDThrown_FinalState += locPIDMultiplexID;
}
}
1556}

1558void DEventWriterROOT::Group_ThrownParticles(const vector<const DMCThrown*>& locMCThrowns_FinalState, const vector<const DMCThrown*>& locMCThrowns_Decaying, vector<const DMCThrown*>& locMCThrownsToSave, map<const DMCThrown*, unsigned int>& locThrownIndexMap) const
1559{
locMCThrownsToSave.clear();
locMCThrownsToSave.insert(locMCThrownsToSave.end(), locMCThrowns_FinalState.begin(), locMCThrowns_FinalState.end());
locMCThrownsToSave.insert(locMCThrownsToSave.end(), locMCThrowns_Decaying.begin(), locMCThrowns_Decaying.end());

//create map of mcthrown to array index
locThrownIndexMap.clear();
for(size_t loc_i = 0; loc_i < locMCThrownsToSave.size(); ++loc_i)
locThrownIndexMap[locMCThrownsToSave[loc_i]] = loc_i;
1568}

1570void DEventWriterROOT::Fill_ThrownInfo(DTreeFillData* locTreeFillData, const DMCReaction* locMCReaction, const DBeamPhoton* locTaggedMCGenBeam, const vector<const DMCThrown*>& locMCThrowns, const map<const DMCThrown*, unsigned int>& locThrownIndexMap, ULong64_t locNumPIDThrown_FinalState, ULong64_t locPIDThrown_Decaying, const DMCThrownMatching* locMCThrownMatching) const
1571{
//THIS MUST BE CALLED FROM WITHIN A LOCK, SO DO NOT PASS IN JEVENTLOOP! //TOO TEMPTING TO DO SOMETHING BAD

//WEIGHT
locTreeFillData->Fill_Single<Float_t>("MCWeight", locMCReaction->weight);

//THROWN BEAM
locTreeFillData->Fill_Single<Int_t>(Build_BranchName("ThrownBeam", "PID"), PDGtype(locMCReaction->beam.PID()));
locTreeFillData->Fill_Single<Float_t>(Build_BranchName("ThrownBeam", "GeneratedEnergy"), locMCReaction->beam.energy());

DVector3 locThrownBeamX3 = locMCReaction->beam.position();
TLorentzVector locThrownBeamTX4(locThrownBeamX3.X(), locThrownBeamX3.Y(), locThrownBeamX3.Z(), locMCReaction->beam.time());
locTreeFillData->Fill_Single<TLorentzVector>(Build_BranchName("ThrownBeam", "X4"), locThrownBeamTX4);

DLorentzVector locThrownBeamP4 = locTaggedMCGenBeam->lorentzMomentum();
TLorentzVector locThrownBeamTP4(locThrownBeamP4.Px(), locThrownBeamP4.Py(), locThrownBeamP4.Pz(), locThrownBeamP4.E());
locTreeFillData->Fill_Single<TLorentzVector>(Build_BranchName("ThrownBeam", "P4"), locThrownBeamTP4);

//THROWN PRODUCTS
locTreeFillData->Fill_Single<UInt_t>("NumThrown", locMCThrowns.size());
for(size_t loc_i = 0; loc_i < locMCThrowns.size(); ++loc_i)
Fill_ThrownParticleData(locTreeFillData, loc_i, locMCThrowns[loc_i], locThrownIndexMap, locMCThrownMatching);

//PID INFO
locTreeFillData->Fill_Single<ULong64_t>("NumPIDThrown_FinalState", locNumPIDThrown_FinalState);
locTreeFillData->Fill_Single<ULong64_t>("PIDThrown_Decaying", locPIDThrown_Decaying);
1597}

1599void DEventWriterROOT::Fill_ThrownParticleData(DTreeFillData* locTreeFillData, unsigned int locArrayIndex, const DMCThrown* locMCThrown, const map<const DMCThrown*, unsigned int>& locThrownIndexMap, const DMCThrownMatching* locMCThrownMatching) const
1600{
string locParticleBranchName = "Thrown";

//IDENTIFIERS
int locParentIndex = -1; //e.g. photoproduced
map<const DMCThrown*, unsigned int>::const_iterator locIterator;
for(locIterator = locThrownIndexMap.begin(); locIterator != locThrownIndexMap.end(); ++locIterator)
{
if(locIterator->first->myid != locMCThrown->parentid)
	continue;
locParentIndex = locIterator->second;
break;
}
locTreeFillData->Fill_Array<Int_t>(Build_BranchName(locParticleBranchName, "ParentIndex"), locParentIndex, locArrayIndex);
locTreeFillData->Fill_Array<Int_t>(Build_BranchName(locParticleBranchName, "PID"), locMCThrown->pdgtype, locArrayIndex);

//MATCHING
if(locMCThrownMatching != NULL__null)
{
Int_t locMatchID = -1;
double locMatchFOM = -1.0;
if(ParticleCharge(locMCThrown->PID()) != 0)
{
	const DChargedTrack* locChargedTrack = locMCThrownMatching->Get_MatchingChargedTrack(locMCThrown, locMatchFOM);
	if(locChargedTrack != NULL__null)
		locMatchID = locChargedTrack->candidateid;
}
else
{
	//Can't use DNeutralShower JObject::id (must use dShowerID): 
		//Matching done with default-tag showers, but pre-select showers are saved to tree: JObject::id's don't match
	const DNeutralShower* locNeutralShower = locMCThrownMatching->Get_MatchingNeutralShower(locMCThrown, locMatchFOM);
	if(locNeutralShower != NULL__null)
		locMatchID = locNeutralShower->dShowerID;
}
locTreeFillData->Fill_Array<Int_t>(Build_BranchName(locParticleBranchName, "MatchID"), locMatchID, locArrayIndex);
locTreeFillData->Fill_Array<Float_t>(Build_BranchName(locParticleBranchName, "MatchFOM"), locMatchFOM, locArrayIndex);
}

//KINEMATICS: THROWN //at the production vertex
TLorentzVector locX4_Thrown(locMCThrown->position().X(), locMCThrown->position().Y(), locMCThrown->position().Z(), locMCThrown->time());
locTreeFillData->Fill_Array<TLorentzVector>(Build_BranchName(locParticleBranchName, "X4"), locX4_Thrown, locArrayIndex);
TLorentzVector locP4_Thrown(locMCThrown->momentum().X(), locMCThrown->momentum().Y(), locMCThrown->momentum().Z(), locMCThrown->energy());
locTreeFillData->Fill_Array<TLorentzVector>(Build_BranchName(locParticleBranchName, "P4"), locP4_Thrown, locArrayIndex);
1644}

1646void DEventWriterROOT::Fill_BeamData(DTreeFillData* locTreeFillData, unsigned int locArrayIndex, const DBeamPhoton* locBeamPhoton, const DVertex* locVertex, const DMCThrownMatching* locMCThrownMatching) const
1647{
string locParticleBranchName = "Beam";

//IDENTIFIER
locTreeFillData->Fill_Array<Int_t>(Build_BranchName(locParticleBranchName, "PID"), PDGtype(locBeamPhoton->PID()), locArrayIndex);

//MATCHING
if(locMCThrownMatching != NULL__null)
{
       // Tag the thrown beam photon by comparing energy, time and counter
       const DBeamPhoton* locBeamPhotonTaggedMCGEN = locMCThrownMatching->Get_TaggedMCGENBeamPhoton();
Bool_t locIsGeneratorFlag = kFALSE;
if (locBeamPhotonTaggedMCGEN != NULL__null)
  locIsGeneratorFlag = (locBeamPhotonTaggedMCGEN->energy() == locBeamPhoton->energy() && locBeamPhotonTaggedMCGEN->dCounter == locBeamPhoton->dCounter && locBeamPhotonTaggedMCGEN->time() == locBeamPhoton->time()) ? kTRUE : kFALSE;
locTreeFillData->Fill_Array<Bool_t>(Build_BranchName(locParticleBranchName, "IsGenerator"), locIsGeneratorFlag, locArrayIndex);
}

//KINEMATICS: MEASURED

//use production vertex, propagate photon time
DVector3 locTargetCenter = locBeamPhoton->position();
DVector3 locProductionVertex = locVertex->dSpacetimeVertex.Vect();
double locDeltaPath = (locProductionVertex - locTargetCenter).Mag();
bool locDownstreamFlag = ((locProductionVertex.Z() - locTargetCenter.Z()) > 0.0);
double locDeltaT = locDownstreamFlag ? locDeltaPath/29.9792458 : -1.0*locDeltaPath/29.9792458;
double locTime = locBeamPhoton->time() + locDeltaT;

TLorentzVector locX4_Measured(locProductionVertex.X(), locProductionVertex.Y(), locProductionVertex.Z(), locTime);
locTreeFillData->Fill_Array<TLorentzVector>(Build_BranchName(locParticleBranchName, "X4_Measured"), locX4_Measured, locArrayIndex);

DLorentzVector locDP4 = locBeamPhoton->lorentzMomentum();
TLorentzVector locP4_Measured(locDP4.Px(), locDP4.Py(), locDP4.Pz(), locDP4.E());
locTreeFillData->Fill_Array<TLorentzVector>(Build_BranchName(locParticleBranchName, "P4_Measured"), locP4_Measured, locArrayIndex);
1680}

1682void DEventWriterROOT::Fill_ChargedHypo(DTreeFillData* locTreeFillData, unsigned int locArrayIndex, const DChargedTrackHypothesis* locChargedTrackHypothesis, const DMCThrownMatching* locMCThrownMatching, const map<const DMCThrown*, unsigned int>& locThrownIndexMap, const DDetectorMatches* locDetectorMatches) const
1683{
string locParticleBranchName = "ChargedHypo";

//ASSOCIATED OBJECTS
auto locTrackTimeBased = locChargedTrackHypothesis->Get_TrackTimeBased();

const DBCALShower* locBCALShower = NULL__null;
if(locChargedTrackHypothesis->Get_BCALShowerMatchParams() != NULL__null)
locBCALShower = locChargedTrackHypothesis->Get_BCALShowerMatchParams()->dBCALShower;

const DFCALShower* locFCALShower = NULL__null;
if(locChargedTrackHypothesis->Get_FCALShowerMatchParams() != NULL__null)
locFCALShower = locChargedTrackHypothesis->Get_FCALShowerMatchParams()->dFCALShower;

//IDENTIFIERS
locTreeFillData->Fill_Array<Int_t>(Build_BranchName(locParticleBranchName, "TrackID"), locTrackTimeBased->candidateid, locArrayIndex);
locTreeFillData->Fill_Array<Int_t>(Build_BranchName(locParticleBranchName, "PID"), PDGtype(locChargedTrackHypothesis->PID()), locArrayIndex);

//MATCHING
if(locMCThrownMatching != NULL__null)
{
Int_t locThrownIndex = -1;
double locMatchFOM = 0.0;
const DMCThrown* locMCThrown = locMCThrownMatching->Get_MatchingMCThrown(locChargedTrackHypothesis, locMatchFOM);
if(locMCThrown != NULL__null)
	locThrownIndex = locThrownIndexMap.find(locMCThrown)->second;
locTreeFillData->Fill_Array<Int_t>(Build_BranchName(locParticleBranchName, "ThrownIndex"), locThrownIndex, locArrayIndex);
}

//KINEMATICS: MEASURED
DVector3 locPosition = locChargedTrackHypothesis->position();
TLorentzVector locTX4_Measured(locPosition.X(), locPosition.Y(), locPosition.Z(), locChargedTrackHypothesis->time());
locTreeFillData->Fill_Array<TLorentzVector>(Build_BranchName(locParticleBranchName, "X4_Measured"), locTX4_Measured, locArrayIndex);

DLorentzVector locDP4 = locChargedTrackHypothesis->lorentzMomentum();
TLorentzVector locP4_Measured(locDP4.Px(), locDP4.Py(), locDP4.Pz(), locDP4.E());
locTreeFillData->Fill_Array<TLorentzVector>(Build_BranchName(locParticleBranchName, "P4_Measured"), locP4_Measured, locArrayIndex);

// Global PID
locTreeFillData->Fill_Array<Float_t>(Build_BranchName(locParticleBranchName, "PIDFOM"), locChargedTrackHypothesis->Get_FOM(), locArrayIndex);

//TRACKING INFO
locTreeFillData->Fill_Array<UInt_t>(Build_BranchName(locParticleBranchName, "NDF_Tracking"), locTrackTimeBased->Ndof, locArrayIndex);
locTreeFillData->Fill_Array<Float_t>(Build_BranchName(locParticleBranchName, "ChiSq_Tracking"), locTrackTimeBased->chisq, locArrayIndex);
locTreeFillData->Fill_Array<UInt_t>(Build_BranchName(locParticleBranchName, "NDF_DCdEdx"), locChargedTrackHypothesis->Get_NDF_DCdEdx(), locArrayIndex);
locTreeFillData->Fill_Array<Float_t>(Build_BranchName(locParticleBranchName, "ChiSq_DCdEdx"), locChargedTrackHypothesis->Get_ChiSq_DCdEdx(), locArrayIndex);
locTreeFillData->Fill_Array<Float_t>(Build_BranchName(locParticleBranchName, "dEdx_CDC"), locTrackTimeBased->ddEdx_CDC_amp, locArrayIndex);
locTreeFillData->Fill_Array<Float_t>(Build_BranchName(locParticleBranchName, "dEdx_CDC_integral"), locTrackTimeBased->ddEdx_CDC, locArrayIndex);
locTreeFillData->Fill_Array<Float_t>(Build_BranchName(locParticleBranchName, "dEdx_FDC"), locTrackTimeBased->ddEdx_FDC, locArrayIndex);

//HIT ENERGY
double locTOFdEdx = (locChargedTrackHypothesis->Get_TOFHitMatchParams() != NULL__null) ? locChargedTrackHypothesis->Get_TOFHitMatchParams()->dEdx : 0.0;
locTreeFillData->Fill_Array<Float_t>(Build_BranchName(locParticleBranchName, "dEdx_TOF"), locTOFdEdx, locArrayIndex);
double locSCdEdx = (locChargedTrackHypothesis->Get_SCHitMatchParams() != NULL__null) ? locChargedTrackHypothesis->Get_SCHitMatchParams()->dEdx : 0.0;
locTreeFillData->Fill_Array<Float_t>(Build_BranchName(locParticleBranchName, "dEdx_ST"), locSCdEdx, locArrayIndex);
double locBCALEnergy = (locBCALShower != NULL__null) ? locBCALShower->E : 0.0;
locTreeFillData->Fill_Array<Float_t>(Build_BranchName(locParticleBranchName, "Energy_BCAL"), locBCALEnergy, locArrayIndex);
double locBCALPreshowerEnergy = (locBCALShower != NULL__null) ? locBCALShower->E_preshower : 0.0;
locTreeFillData->Fill_Array<Float_t>(Build_BranchName(locParticleBranchName, "Energy_BCALPreshower"), locBCALPreshowerEnergy, locArrayIndex);
if(BCAL_VERBOSE_OUTPUT) {
double locBCALLayer2Energy = (locBCALShower != NULL__null) ? locBCALShower->E_L2 : 0.0;
locTreeFillData->Fill_Array<Float_t>(Build_BranchName(locParticleBranchName, "Energy_BCALLayer2"), locBCALLayer2Energy, locArrayIndex);
double locBCALLayer3Energy = (locBCALShower != NULL__null) ? locBCALShower->E_L3 : 0.0;
locTreeFillData->Fill_Array<Float_t>(Build_BranchName(locParticleBranchName, "Energy_BCALLayer3"), locBCALLayer3Energy, locArrayIndex);
double locBCALLayer4Energy = (locBCALShower != NULL__null) ? locBCALShower->E_L4 : 0.0;
locTreeFillData->Fill_Array<Float_t>(Build_BranchName(locParticleBranchName, "Energy_BCALLayer4"), locBCALLayer4Energy, locArrayIndex);
}

double locFCALEnergy = (locFCALShower != NULL__null) ? locFCALShower->getEnergy() : 0.0;
locTreeFillData->Fill_Array<Float_t>(Build_BranchName(locParticleBranchName, "Energy_FCAL"), locFCALEnergy, locArrayIndex);

//double locCCALEnergy = (locCCALShower != NULL) ? locCCALShower->getEnergy() : 0.0;
//locTreeFillData->Fill_Array<Float_t>(Build_BranchName(locParticleBranchName, "Energy_CCAL"), locCCALEnergy, locArrayIndex);

//SHOWER PROPERTIES
double locSigLongBCAL = (locBCALShower != NULL__null) ? locBCALShower->sigLong : 0.0;
double locSigThetaBCAL = (locBCALShower != NULL__null) ? locBCALShower->sigTheta : 0.0;
double locSigTransBCAL = (locBCALShower != NULL__null) ? locBCALShower->sigTrans : 0.0;
double locRMSTimeBCAL = (locBCALShower != NULL__null) ? locBCALShower->rmsTime : 0.0;
locTreeFillData->Fill_Array<Float_t>(Build_BranchName(locParticleBranchName, "SigLong_BCAL"), locSigLongBCAL, locArrayIndex);
locTreeFillData->Fill_Array<Float_t>(Build_BranchName(locParticleBranchName, "SigTheta_BCAL"), locSigThetaBCAL, locArrayIndex);
locTreeFillData->Fill_Array<Float_t>(Build_BranchName(locParticleBranchName, "SigTrans_BCAL"), locSigTransBCAL, locArrayIndex);
locTreeFillData->Fill_Array<Float_t>(Build_BranchName(locParticleBranchName, "RMSTime_BCAL"), locRMSTimeBCAL, locArrayIndex);

double locE1E9FCAL = (locFCALShower != NULL__null) ? locFCALShower->getE1E9() : 0.0;
double locE9E25FCAL = (locFCALShower != NULL__null) ? locFCALShower->getE9E25() : 0.0;
double locSumUFCAL = (locFCALShower != NULL__null) ? locFCALShower->getSumU() : 0.0;
double locSumVFCAL = (locFCALShower != NULL__null) ? locFCALShower->getSumV() : 0.0;
locTreeFillData->Fill_Array<Float_t>(Build_BranchName(locParticleBranchName, "E1E9_FCAL"), locE1E9FCAL, locArrayIndex);
locTreeFillData->Fill_Array<Float_t>(Build_BranchName(locParticleBranchName, "E9E25_FCAL"), locE9E25FCAL, locArrayIndex);
locTreeFillData->Fill_Array<Float_t>(Build_BranchName(locParticleBranchName, "SumU_FCAL"), locSumUFCAL, locArrayIndex);
locTreeFillData->Fill_Array<Float_t>(Build_BranchName(locParticleBranchName, "SumV_FCAL"), locSumVFCAL, locArrayIndex);

//TIMING INFO
locTreeFillData->Fill_Array<Float_t>(Build_BranchName(locParticleBranchName, "HitTime"), locChargedTrackHypothesis->t1(), locArrayIndex);
double locStartTimeError = locChargedTrackHypothesis->t0_err();
double locRFDeltaTVariance = (*locChargedTrackHypothesis->errorMatrix())(6, 6) + locStartTimeError*locStartTimeError;
locTreeFillData->Fill_Array<Float_t>(Build_BranchName(locParticleBranchName, "RFDeltaTVar"), locRFDeltaTVariance, locArrayIndex);

//MEASURED PID INFO
locTreeFillData->Fill_Array<Float_t>(Build_BranchName(locParticleBranchName, "Beta_Timing"), locChargedTrackHypothesis->measuredBeta(), locArrayIndex);
locTreeFillData->Fill_Array<Float_t>(Build_BranchName(locParticleBranchName, "ChiSq_Timing"), locChargedTrackHypothesis->Get_ChiSq_Timing(), locArrayIndex);
locTreeFillData->Fill_Array<UInt_t>(Build_BranchName(locParticleBranchName, "NDF_Timing"), locChargedTrackHypothesis->Get_NDF_Timing(), locArrayIndex);

//SHOWER MATCHING: BCAL
double locTrackBCAL_DeltaPhi = 999.0, locTrackBCAL_DeltaZ = 999.0;
if(locChargedTrackHypothesis->Get_BCALShowerMatchParams() != NULL__null)
{
locTrackBCAL_DeltaPhi = locChargedTrackHypothesis->Get_BCALShowerMatchParams()->dDeltaPhiToShower;
locTrackBCAL_DeltaZ = locChargedTrackHypothesis->Get_BCALShowerMatchParams()->dDeltaZToShower;
}
locTreeFillData->Fill_Array<Float_t>(Build_BranchName(locParticleBranchName, "TrackBCAL_DeltaPhi"), locTrackBCAL_DeltaPhi, locArrayIndex);
locTreeFillData->Fill_Array<Float_t>(Build_BranchName(locParticleBranchName, "TrackBCAL_DeltaZ"), locTrackBCAL_DeltaZ, locArrayIndex);

//SHOWER MATCHING: FCAL
double locDOCAToShower_FCAL = 999.0;
if(locChargedTrackHypothesis->Get_FCALShowerMatchParams() != NULL__null)
locDOCAToShower_FCAL = locChargedTrackHypothesis->Get_FCALShowerMatchParams()->dDOCAToShower;
locTreeFillData->Fill_Array<Float_t>(Build_BranchName(locParticleBranchName, "TrackFCAL_DOCA"), locDOCAToShower_FCAL, locArrayIndex);

// DIRC
if(DIRC_OUTPUT) {
int locDIRCNumPhotons = 0;
double locDIRCExtrapolatedX = 999;
double locDIRCExtrapolatedY = 999;
double locDIRCThetaC = 999.;
double locDIRCLele = 999.;
double locDIRCLpi = 999.;
double locDIRCLk = 999.;
double locDIRCLp = 999.;
auto locDIRCMatchParams = locChargedTrackHypothesis->Get_DIRCMatchParams();
if(locDIRCMatchParams != NULL__null) {
        locDIRCExtrapolatedX = locDIRCMatchParams->dExtrapolatedPos.X();
	locDIRCExtrapolatedY = locDIRCMatchParams->dExtrapolatedPos.Y();
	locDIRCNumPhotons = locDIRCMatchParams->dNPhotons;
	locDIRCThetaC = locDIRCMatchParams->dThetaC;
	locDIRCLele =  locDIRCMatchParams->dLikelihoodElectron;
	locDIRCLpi =  locDIRCMatchParams->dLikelihoodPion;	
	locDIRCLk =  locDIRCMatchParams->dLikelihoodKaon;
	locDIRCLp =  locDIRCMatchParams->dLikelihoodProton;
}
locTreeFillData->Fill_Array<Int_t>(Build_BranchName(locParticleBranchName, "NumPhotons_DIRC"), locDIRCNumPhotons, locArrayIndex);
locTreeFillData->Fill_Array<Float_t>(Build_BranchName(locParticleBranchName, "ExtrapolatedX_DIRC"), locDIRCExtrapolatedX, locArrayIndex);
locTreeFillData->Fill_Array<Float_t>(Build_BranchName(locParticleBranchName, "ExtrapolatedY_DIRC"), locDIRCExtrapolatedY, locArrayIndex);
locTreeFillData->Fill_Array<Float_t>(Build_BranchName(locParticleBranchName, "ThetaC_DIRC"), locDIRCThetaC, locArrayIndex);
locTreeFillData->Fill_Array<Float_t>(Build_BranchName(locParticleBranchName, "Lele_DIRC"), locDIRCLele, locArrayIndex);
locTreeFillData->Fill_Array<Float_t>(Build_BranchName(locParticleBranchName, "Lpi_DIRC"), locDIRCLpi, locArrayIndex);
locTreeFillData->Fill_Array<Float_t>(Build_BranchName(locParticleBranchName, "Lk_DIRC"), locDIRCLk, locArrayIndex);
locTreeFillData->Fill_Array<Float_t>(Build_BranchName(locParticleBranchName, "Lp_DIRC"), locDIRCLp, locArrayIndex);
}
1833}

1835void DEventWriterROOT::Fill_NeutralHypo(DTreeFillData* locTreeFillData, unsigned int locArrayIndex, const DNeutralParticleHypothesis* locNeutralParticleHypothesis, const DMCThrownMatching* locMCThrownMatching, const map<const DMCThrown*, unsigned int>& locThrownIndexMap, const DDetectorMatches* locDetectorMatches) const
1836{
string locParticleBranchName = "NeutralHypo";
const DNeutralShower* locNeutralShower = locNeutralParticleHypothesis->Get_NeutralShower();

//ASSOCIATED OBJECTS
const DBCALShower* locBCALShower = NULL__null;
locNeutralShower->GetSingle(locBCALShower);
const DFCALShower* locFCALShower = NULL__null;
locNeutralShower->GetSingle(locFCALShower);
const DCCALShower* locCCALShower = NULL__null;
locNeutralShower->GetSingle(locCCALShower);

//IDENTIFIERS
Particle_t locPID = locNeutralParticleHypothesis->PID();
locTreeFillData->Fill_Array<Int_t>(Build_BranchName(locParticleBranchName, "NeutralID"), locNeutralShower->dShowerID, locArrayIndex);
locTreeFillData->Fill_Array<Int_t>(Build_BranchName(locParticleBranchName, "PID"), PDGtype(locPID), locArrayIndex);

//MATCHING
if(locMCThrownMatching != NULL__null)
{
Int_t locThrownIndex = -1;
double locMatchFOM = 0.0;
const DMCThrown* locMCThrown = locMCThrownMatching->Get_MatchingMCThrown(locNeutralParticleHypothesis, locMatchFOM);
if(locMCThrown != NULL__null)
	locThrownIndex = locThrownIndexMap.find(locMCThrown)->second;
locTreeFillData->Fill_Array<Int_t>(Build_BranchName(locParticleBranchName, "ThrownIndex"), locThrownIndex, locArrayIndex);
}

//KINEMATICS: MEASURED
DVector3 locPosition = locNeutralParticleHypothesis->position();
TLorentzVector locX4_Measured(locPosition.X(), locPosition.Y(), locPosition.Z(), locNeutralParticleHypothesis->time());
locTreeFillData->Fill_Array<TLorentzVector>(Build_BranchName(locParticleBranchName, "X4_Measured"), locX4_Measured, locArrayIndex);

DLorentzVector locDP4 = locNeutralParticleHypothesis->lorentzMomentum();
TLorentzVector locP4_Measured(locDP4.Px(), locDP4.Py(), locDP4.Pz(), locDP4.E());
locTreeFillData->Fill_Array<TLorentzVector>(Build_BranchName(locParticleBranchName, "P4_Measured"), locP4_Measured, locArrayIndex);

//MEASURED PID INFO
locTreeFillData->Fill_Array<Float_t>(Build_BranchName(locParticleBranchName, "Beta_Timing"), locNeutralParticleHypothesis->measuredBeta(), locArrayIndex);
locTreeFillData->Fill_Array<Float_t>(Build_BranchName(locParticleBranchName, "ChiSq_Timing"), locNeutralParticleHypothesis->Get_ChiSq(), locArrayIndex);
locTreeFillData->Fill_Array<UInt_t>(Build_BranchName(locParticleBranchName, "NDF_Timing"), locNeutralParticleHypothesis->Get_NDF(), locArrayIndex);
locTreeFillData->Fill_Array<Float_t>(Build_BranchName(locParticleBranchName, "ShowerQuality"), locNeutralShower->dQuality, locArrayIndex);

//SHOWER ENERGY
DetectorSystem_t locDetector = locNeutralShower->dDetectorSystem;
double locBCALEnergy = (locDetector == SYS_BCAL) ? locNeutralShower->dEnergy : 0.0;
locTreeFillData->Fill_Array<Float_t>(Build_BranchName(locParticleBranchName, "Energy_BCAL"), locBCALEnergy, locArrayIndex);
double locBCALPreshowerEnergy = (locDetector == SYS_BCAL) ? static_cast<const DBCALShower*>(locNeutralShower->dBCALFCALShower)->E_preshower : 0.0;
locTreeFillData->Fill_Array<Float_t>(Build_BranchName(locParticleBranchName, "Energy_BCALPreshower"), locBCALPreshowerEnergy, locArrayIndex);
if(BCAL_VERBOSE_OUTPUT) {
double locBCALLayer2Energy = (locBCALShower != NULL__null) ? locBCALShower->E_L2 : 0.0;
locTreeFillData->Fill_Array<Float_t>(Build_BranchName(locParticleBranchName, "Energy_BCALLayer2"), locBCALLayer2Energy, locArrayIndex);
double locBCALLayer3Energy = (locBCALShower != NULL__null) ? locBCALShower->E_L3 : 0.0;
locTreeFillData->Fill_Array<Float_t>(Build_BranchName(locParticleBranchName, "Energy_BCALLayer3"), locBCALLayer3Energy, locArrayIndex);
double locBCALLayer4Energy = (locBCALShower != NULL__null) ? locBCALShower->E_L4 : 0.0;
locTreeFillData->Fill_Array<Float_t>(Build_BranchName(locParticleBranchName, "Energy_BCALLayer4"), locBCALLayer4Energy, locArrayIndex);
}

double locFCALEnergy = (locDetector == SYS_FCAL) ? locNeutralShower->dEnergy : 0.0;
locTreeFillData->Fill_Array<Float_t>(Build_BranchName(locParticleBranchName, "Energy_FCAL"), locFCALEnergy, locArrayIndex);

double locCCALEnergy = (locDetector == SYS_CCAL) ? locNeutralShower->dEnergy : 0.0;
locTreeFillData->Fill_Array<Float_t>(Build_BranchName(locParticleBranchName, "Energy_CCAL"), locCCALEnergy, locArrayIndex);

//SHOWER POSITION
DLorentzVector locHitDX4 = locNeutralShower->dSpacetimeVertex;
TLorentzVector locTX4_Shower(locHitDX4.X(), locHitDX4.Y(), locHitDX4.Z(), locHitDX4.T());
locTreeFillData->Fill_Array<TLorentzVector>(Build_BranchName(locParticleBranchName, "X4_Shower"), locTX4_Shower, locArrayIndex);

//SHOWER PROPERTIES
double locSigLongBCAL = (locBCALShower != NULL__null) ? locBCALShower->sigLong : 0.0;
double locSigThetaBCAL = (locBCALShower != NULL__null) ? locBCALShower->sigTheta : 0.0;
double locSigTransBCAL = (locBCALShower != NULL__null) ? locBCALShower->sigTrans : 0.0;
double locRMSTimeBCAL = (locBCALShower != NULL__null) ? locBCALShower->rmsTime : 0.0;
locTreeFillData->Fill_Array<Float_t>(Build_BranchName(locParticleBranchName, "SigLong_BCAL"), locSigLongBCAL, locArrayIndex);
locTreeFillData->Fill_Array<Float_t>(Build_BranchName(locParticleBranchName, "SigTheta_BCAL"), locSigThetaBCAL, locArrayIndex);
locTreeFillData->Fill_Array<Float_t>(Build_BranchName(locParticleBranchName, "SigTrans_BCAL"), locSigTransBCAL, locArrayIndex);
locTreeFillData->Fill_Array<Float_t>(Build_BranchName(locParticleBranchName, "RMSTime_BCAL"), locRMSTimeBCAL, locArrayIndex);

if(FCAL_VERBOSE_OUTPUT) {
double locE1E9FCAL = (locFCALShower != NULL__null) ? locFCALShower->getE1E9() : 0.0;
double locE9E25FCAL = (locFCALShower != NULL__null) ? locFCALShower->getE9E25() : 0.0;
double locSumUFCAL = (locFCALShower != NULL__null) ? locFCALShower->getSumU() : 0.0;
double locSumVFCAL = (locFCALShower != NULL__null) ? locFCALShower->getSumV() : 0.0;
locTreeFillData->Fill_Array<Float_t>(Build_BranchName(locParticleBranchName, "E1E9_FCAL"), locE1E9FCAL, locArrayIndex);
locTreeFillData->Fill_Array<Float_t>(Build_BranchName(locParticleBranchName, "E9E25_FCAL"), locE9E25FCAL, locArrayIndex);
locTreeFillData->Fill_Array<Float_t>(Build_BranchName(locParticleBranchName, "SumU_FCAL"), locSumUFCAL, locArrayIndex);
locTreeFillData->Fill_Array<Float_t>(Build_BranchName(locParticleBranchName, "SumV_FCAL"), locSumVFCAL, locArrayIndex);
}

//Track DOCA to Shower - BCAL
double locNearestTrackBCALDeltaPhi = 999.0, locNearestTrackBCALDeltaZ = 999.0;
if(locBCALShower != NULL__null)
{
if(!locDetectorMatches->Get_DistanceToNearestTrack(locBCALShower, locNearestTrackBCALDeltaPhi, locNearestTrackBCALDeltaZ))
{
	locNearestTrackBCALDeltaPhi = 999.0;
	locNearestTrackBCALDeltaZ = 999.0;
}
else if((locNearestTrackBCALDeltaPhi > 999.0) || (locNearestTrackBCALDeltaZ > 999.0))
{
	locNearestTrackBCALDeltaPhi = 999.0;
	locNearestTrackBCALDeltaZ = 999.0;
}
}
locTreeFillData->Fill_Array<Float_t>(Build_BranchName(locParticleBranchName, "TrackBCAL_DeltaPhi"), locNearestTrackBCALDeltaPhi, locArrayIndex);
locTreeFillData->Fill_Array<Float_t>(Build_BranchName(locParticleBranchName, "TrackBCAL_DeltaZ"), locNearestTrackBCALDeltaZ, locArrayIndex);

//Track DOCA to Shower - FCAL
double locDistanceToNearestTrack_FCAL = 999.0;
if(locFCALShower != NULL__null)
{
if(!locDetectorMatches->Get_DistanceToNearestTrack(locFCALShower, locDistanceToNearestTrack_FCAL))
	locDistanceToNearestTrack_FCAL = 999.0;
if(locDistanceToNearestTrack_FCAL > 999.0)
	locDistanceToNearestTrack_FCAL = 999.0;
}
locTreeFillData->Fill_Array<Float_t>(Build_BranchName(locParticleBranchName, "TrackFCAL_DOCA"), locDistanceToNearestTrack_FCAL, locArrayIndex);

//PHOTON PID INFO
double locStartTimeError = locNeutralParticleHypothesis->t0_err();
double locPhotonRFDeltaTVar = (*locNeutralParticleHypothesis->errorMatrix())(6, 6) + locStartTimeError*locStartTimeError;
if(locPID != Gamma)
locPhotonRFDeltaTVar = 0.0;
locTreeFillData->Fill_Array<Float_t>(Build_BranchName(locParticleBranchName, "PhotonRFDeltaTVar"), locPhotonRFDeltaTVar, locArrayIndex);
1961}

1963void DEventWriterROOT::Fill_ComboData(DTreeFillData* locTreeFillData, const DReaction* locReaction, const DParticleCombo* locParticleCombo, unsigned int locComboIndex, const map<pair<oid_t, Particle_t>, size_t>& locObjectToArrayIndexMap) const
1964{
//MAIN CLASSES
const DKinFitResults* locKinFitResults = locParticleCombo->Get_KinFitResults();
const DEventRFBunch* locEventRFBunch = locParticleCombo->Get_EventRFBunch();

//IS COMBO CUT
locTreeFillData->Fill_Array<Bool_t>("IsComboCut", kFALSE, locComboIndex);

//RF INFO
double locRFTime = (locEventRFBunch != NULL__null) ? locEventRFBunch->dTime : numeric_limits<double>::quiet_NaN();
locTreeFillData->Fill_Array<Float_t>("RFTime_Measured", locRFTime, locComboIndex);

//KINFIT INFO
DKinFitType locKinFitType = locReaction->Get_KinFitType();
bool locKinFitFlag = (locKinFitType != d_NoFit);
if(locKinFitFlag)
{
if(locKinFitResults != NULL__null)
{
	locTreeFillData->Fill_Array<Float_t>("ChiSq_KinFit", locKinFitResults->Get_ChiSq(), locComboIndex);
	locTreeFillData->Fill_Array<UInt_t>("NDF_KinFit", locKinFitResults->Get_NDF(), locComboIndex);
	if((locKinFitType == d_SpacetimeFit) || (locKinFitType == d_P4AndSpacetimeFit))
	{
		double locRFTime_KinFit = -9.9E9; //NOT IMPLEMENTED YET
		locTreeFillData->Fill_Array<Float_t>("RFTime_KinFit", locRFTime_KinFit, locComboIndex);
	}
}
else
{
	locTreeFillData->Fill_Array<Float_t>("ChiSq_KinFit", 0.0, locComboIndex);
	locTreeFillData->Fill_Array<UInt_t>("NDF_KinFit", 0, locComboIndex);
	if((locKinFitType == d_SpacetimeFit) || (locKinFitType == d_P4AndSpacetimeFit))
		locTreeFillData->Fill_Array<Float_t>("RFTime_KinFit", -9.9E9, locComboIndex);
}
}

//STEP DATA
for(size_t loc_i = 0; loc_i < locParticleCombo->Get_NumParticleComboSteps(); ++loc_i)
Fill_ComboStepData(locTreeFillData, locReaction, locParticleCombo, loc_i, locComboIndex, locKinFitType, locObjectToArrayIndexMap);
2003}

2005void DEventWriterROOT::Fill_ComboStepData(DTreeFillData* locTreeFillData, const DReaction* locReaction, const DParticleCombo* locParticleCombo, unsigned int locStepIndex, unsigned int locComboIndex, DKinFitType locKinFitType, const map<pair<oid_t, Particle_t>, size_t>& locObjectToArrayIndexMap) const
2006{
auto locReactionVertexInfo = dVertexInfoMap.find(locReaction)->second;
auto locReactionStep = locReaction->Get_ReactionStep(locStepIndex);
const TList* locUserInfo = dTreeInterfaceMap.find(locReaction)->second->Get_UserInfo(); //No Lock!  But this should be unchanging at this point anyway
const TMap* locPositionToNameMap = (TMap*)locUserInfo->FindObject("PositionToNameMap");

auto locParticleComboStep = locParticleCombo->Get_ParticleComboStep(locStepIndex);
DLorentzVector locStepX4 = locParticleComboStep->Get_SpacetimeVertex();
TLorentzVector locStepTX4(locStepX4.X(), locStepX4.Y(), locStepX4.Z(), locStepX4.T());

//beam & production vertex
Particle_t locInitialPID = locReactionStep->Get_InitialPID();
const DKinematicData* locInitialParticle = locParticleComboStep->Get_InitialParticle();
const DBeamPhoton* locBeamPhoton = dynamic_cast<const DBeamPhoton*>(locInitialParticle);
if(locBeamPhoton != NULL__null)
{
const DKinematicData* locInitParticleMeasured = locParticleComboStep->Get_InitialParticle_Measured();
const DBeamPhoton* locMeasuredBeamPhoton = dynamic_cast<const DBeamPhoton*>(locInitParticleMeasured);

//get array index
pair<oid_t, Particle_t> locBeamPair(locMeasuredBeamPhoton->id, locMeasuredBeamPhoton->PID());
size_t locBeamIndex = locObjectToArrayIndexMap.find(locBeamPair)->second;

Fill_ComboBeamData(locTreeFillData, locComboIndex, locBeamPhoton, locBeamIndex, locKinFitType);
}
else //decaying
{
//get the branch name
ostringstream locPositionStream;
locPositionStream << locStepIndex << "_-1";
TObjString* locObjString = (TObjString*)locPositionToNameMap->GetValue(locPositionStream.str().c_str());
string locParticleBranchName = (const char*)(locObjString->GetString());

auto locP4FitFlag = ((locKinFitType == d_P4Fit) || (locKinFitType == d_P4AndVertexFit) || (locKinFitType == d_P4AndSpacetimeFit));
if(IsFixedMass(locInitialPID) && locReactionStep->Get_KinFitConstrainInitMassFlag() && locP4FitFlag)
{
	TLorentzVector locDecayP4;
	if(locInitialParticle == NULL__null)
	{
		//fit failed to converge, calc from other particles
		DLorentzVector locDecayDP4 = dAnalysisUtilities->Calc_FinalStateP4(locReaction, locParticleCombo, locStepIndex, false);
		locDecayDP4.SetE(sqrt(locDecayDP4.Vect().Mag2() + ParticleMass(locInitialPID)*ParticleMass(locInitialPID)));
		locDecayP4.SetPxPyPzE(locDecayDP4.Px(), locDecayDP4.Py(), locDecayDP4.Pz(), locDecayDP4.E());
	}
	else
		locDecayP4.SetPxPyPzE(locInitialParticle->momentum().X(), locInitialParticle->momentum().Y(), locInitialParticle->momentum().Z(), locInitialParticle->energy());
	locTreeFillData->Fill_Array<TLorentzVector>(Build_BranchName(locParticleBranchName, "P4_KinFit"), locDecayP4, locComboIndex);
}

if((locStepIndex == 0) || IsDetachedVertex(locInitialPID))
	locTreeFillData->Fill_Array<TLorentzVector>(Build_BranchName(locParticleBranchName, "X4"), locStepTX4, locComboIndex);

auto locStepVertexInfo = locReactionVertexInfo->Get_StepVertexInfo(locStepIndex);
auto locParentVertexInfo = locStepVertexInfo->Get_ParentVertexInfo();
auto locVertexKinFitFlag = ((locKinFitType != d_P4Fit) && (locKinFitType != d_NoFit));
if(IsDetachedVertex(locInitialPID) && locVertexKinFitFlag && (locParentVertexInfo != nullptr) && locStepVertexInfo->Get_FittableVertexFlag() && locParentVertexInfo->Get_FittableVertexFlag())
{
	auto locKinFitParticle = locParticleComboStep->Get_InitialKinFitParticle();
	auto locPathLengthSigma = locKinFitParticle->Get_PathLengthUncertainty();
	locTreeFillData->Fill_Array<Float_t>(Build_BranchName(locParticleBranchName, "PathLengthSigma"), locPathLengthSigma, locComboIndex);
}
}

//final state particles
for(size_t loc_i = 0; loc_i < locParticleComboStep->Get_NumFinalParticles(); ++loc_i)
{
Particle_t locPID = locReactionStep->Get_FinalPID(loc_i);
const DKinematicData* locKinematicData = locParticleComboStep->Get_FinalParticle(loc_i);
const DKinematicData* locKinematicData_Measured = locParticleComboStep->Get_FinalParticle_Measured(loc_i);

//decaying particle
if(DAnalysis::Get_DecayStepIndex(locReaction, locStepIndex, loc_i) >= 0)
	continue;

//get the branch name
ostringstream locPositionStream;
locPositionStream << locStepIndex << "_" << loc_i;
TObjString* locObjString = (TObjString*)locPositionToNameMap->GetValue(locPositionStream.str().c_str());
string locParticleBranchName = (const char*)(locObjString->GetString());

//missing particle
if(locReactionStep->Get_MissingParticleIndex() == int(loc_i))
{
	if((locKinFitType == d_P4Fit) || (locKinFitType == d_P4AndVertexFit) || (locKinFitType == d_P4AndSpacetimeFit))
	{
		TLorentzVector locMissingP4;
		if(locKinematicData == NULL__null)
		{
			//fit failed to converge, calc from other particles
			DLorentzVector locMissingDP4 = dAnalysisUtilities->Calc_MissingP4(locReaction, locParticleCombo, false);
			locMissingDP4.SetE(sqrt(locMissingDP4.Vect().Mag2() + ParticleMass(locPID)*ParticleMass(locPID)));
			locMissingP4.SetPxPyPzE(locMissingDP4.Px(), locMissingDP4.Py(), locMissingDP4.Pz(), locMissingDP4.E());
		}
		else
			locMissingP4.SetPxPyPzE(locKinematicData->momentum().X(), locKinematicData->momentum().Y(), locKinematicData->momentum().Z(), locKinematicData->energy());

		locTreeFillData->Fill_Array<TLorentzVector>(Build_BranchName(locParticleBranchName, "P4_KinFit"), locMissingP4, locComboIndex);
	}
	continue;
}

//fill the data
if(ParticleCharge(locPID) == 0)
{
	const DNeutralParticleHypothesis* locNeutralHypo = dynamic_cast<const DNeutralParticleHypothesis*>(locKinematicData);
	const DNeutralParticleHypothesis* locMeasuredNeutralHypo = dynamic_cast<const DNeutralParticleHypothesis*>(locKinematicData_Measured);

	//get array index
	const DNeutralShower* locNeutralShower = locNeutralHypo->Get_NeutralShower();
	pair<oid_t, Particle_t> locNeutralPair(locNeutralShower->id, locNeutralHypo->PID());
	size_t locNeutralIndex = locObjectToArrayIndexMap.find(locNeutralPair)->second;

	Fill_ComboNeutralData(locTreeFillData, locComboIndex, locParticleBranchName, locMeasuredNeutralHypo, locNeutralHypo, locNeutralIndex, locKinFitType);
}
else
{
	const DChargedTrackHypothesis* locChargedHypo = dynamic_cast<const DChargedTrackHypothesis*>(locKinematicData);
	const DChargedTrackHypothesis* locMeasuredChargedHypo = dynamic_cast<const DChargedTrackHypothesis*>(locKinematicData_Measured);

	//get array index
	const DTrackTimeBased* locTrackTimeBased = locChargedHypo->Get_TrackTimeBased();
	pair<oid_t, Particle_t> locTrackPair(locTrackTimeBased->id, locChargedHypo->PID());
	size_t locChargedIndex = locObjectToArrayIndexMap.find(locTrackPair)->second;

	Fill_ComboChargedData(locTreeFillData, locComboIndex, locParticleBranchName, locMeasuredChargedHypo, locChargedHypo, locChargedIndex, locKinFitType);
}
}
2133}

2135void DEventWriterROOT::Fill_ComboBeamData(DTreeFillData* locTreeFillData, unsigned int locComboIndex, const DBeamPhoton* locBeamPhoton, size_t locBeamIndex, DKinFitType locKinFitType) const
2136{
string locParticleBranchName = "ComboBeam";

//IDENTIFIER
locTreeFillData->Fill_Array<Int_t>(Build_BranchName(locParticleBranchName, "BeamIndex"), locBeamIndex, locComboIndex);

//KINEMATICS: KINFIT
if(locKinFitType != d_NoFit)
{
if(locKinFitType != d_P4Fit)
{
	DVector3 locPosition = locBeamPhoton->position();
	TLorentzVector locX4_KinFit(locPosition.X(), locPosition.Y(), locPosition.Z(), locBeamPhoton->time());
	locTreeFillData->Fill_Array<TLorentzVector>(Build_BranchName(locParticleBranchName, "X4_KinFit"), locX4_KinFit, locComboIndex);
}

//if charged, bends in b-field, update p4 when vertex changes
if(((locKinFitType != d_VertexFit) && (locKinFitType != d_SpacetimeFit)) || (ParticleCharge(locBeamPhoton->PID()) != 0))
{
	DLorentzVector locDP4 = locBeamPhoton->lorentzMomentum();
	TLorentzVector locP4_KinFit(locDP4.Px(), locDP4.Py(), locDP4.Pz(), locDP4.E());
	locTreeFillData->Fill_Array<TLorentzVector>(Build_BranchName(locParticleBranchName, "P4_KinFit"), locP4_KinFit, locComboIndex);
}
}
2160}

2162void DEventWriterROOT::Fill_ComboChargedData(DTreeFillData* locTreeFillData, unsigned int locComboIndex, string locParticleBranchName, const DChargedTrackHypothesis* locMeasuredChargedHypo, const DChargedTrackHypothesis* locChargedHypo, size_t locChargedIndex, DKinFitType locKinFitType) const
2163{
//IDENTIFIER
locTreeFillData->Fill_Array<Int_t>(Build_BranchName(locParticleBranchName, "ChargedIndex"), locChargedIndex, locComboIndex);

//MEASURED PID
locTreeFillData->Fill_Array<Float_t>(Build_BranchName(locParticleBranchName, "Beta_Timing_Measured"), locMeasuredChargedHypo->measuredBeta(), locComboIndex);
locTreeFillData->Fill_Array<Float_t>(Build_BranchName(locParticleBranchName, "ChiSq_Timing_Measured"), locMeasuredChargedHypo->Get_ChiSq_Timing(), locComboIndex);

//KINFIT PID
if((locKinFitType != d_NoFit) && (locKinFitType != d_SpacetimeFit) && (locKinFitType != d_P4AndSpacetimeFit))
{
locTreeFillData->Fill_Array<Float_t>(Build_BranchName(locParticleBranchName, "Beta_Timing_KinFit"), locChargedHypo->measuredBeta(), locComboIndex);
locTreeFillData->Fill_Array<Float_t>(Build_BranchName(locParticleBranchName, "ChiSq_Timing_KinFit"), locChargedHypo->Get_ChiSq_Timing(), locComboIndex);
}

//KINFIT
if(locKinFitType != d_NoFit)
{
//KINEMATICS
if(locKinFitType != d_P4Fit)
{
	DVector3 locPosition = locChargedHypo->position();
	TLorentzVector locX4_KinFit(locPosition.X(), locPosition.Y(), locPosition.Z(), locChargedHypo->time());
	locTreeFillData->Fill_Array<TLorentzVector>(Build_BranchName(locParticleBranchName, "X4_KinFit"), locX4_KinFit, locComboIndex);
}

//update even if vertex-only fit, because charged momentum propagated through b-field
DLorentzVector locDP4 = locChargedHypo->lorentzMomentum();
TLorentzVector locP4_KinFit(locDP4.Px(), locDP4.Py(), locDP4.Pz(), locDP4.E());
locTreeFillData->Fill_Array<TLorentzVector>(Build_BranchName(locParticleBranchName, "P4_KinFit"), locP4_KinFit, locComboIndex);
}
2194}

2196void DEventWriterROOT::Fill_ComboNeutralData(DTreeFillData* locTreeFillData, unsigned int locComboIndex, string locParticleBranchName, const DNeutralParticleHypothesis* locMeasuredNeutralHypo, const DNeutralParticleHypothesis* locNeutralHypo, size_t locNeutralIndex, DKinFitType locKinFitType) const
2197{
//IDENTIFIER
locTreeFillData->Fill_Array<Int_t>(Build_BranchName(locParticleBranchName, "NeutralIndex"), locNeutralIndex, locComboIndex);

//KINEMATICS: MEASURED
DVector3 locPosition = locMeasuredNeutralHypo->position();
TLorentzVector locX4_Measured(locPosition.X(), locPosition.Y(), locPosition.Z(), locMeasuredNeutralHypo->time());
locTreeFillData->Fill_Array<TLorentzVector>(Build_BranchName(locParticleBranchName, "X4_Measured"), locX4_Measured, locComboIndex);

DLorentzVector locDP4 = locMeasuredNeutralHypo->lorentzMomentum();
TLorentzVector locP4_Measured(locDP4.Px(), locDP4.Py(), locDP4.Pz(), locDP4.E());
locTreeFillData->Fill_Array<TLorentzVector>(Build_BranchName(locParticleBranchName, "P4_Measured"), locP4_Measured, locComboIndex);

//MEASURED PID INFO
locTreeFillData->Fill_Array<Float_t>(Build_BranchName(locParticleBranchName, "Beta_Timing_Measured"), locMeasuredNeutralHypo->measuredBeta(), locComboIndex);
if(locParticleBranchName.substr(0, 6) == "Photon")
locTreeFillData->Fill_Array<Float_t>(Build_BranchName(locParticleBranchName, "ChiSq_Timing_Measured"), locMeasuredNeutralHypo->Get_ChiSq(), locComboIndex);

//KINFIT PID INFO
if((locKinFitType != d_NoFit) && (locKinFitType != d_SpacetimeFit) && (locKinFitType != d_P4AndSpacetimeFit))
{
locTreeFillData->Fill_Array<Float_t>(Build_BranchName(locParticleBranchName, "Beta_Timing_KinFit"), locNeutralHypo->measuredBeta(), locComboIndex);
if(locParticleBranchName.substr(0, 6) == "Photon")
	locTreeFillData->Fill_Array<Float_t>(Build_BranchName(locParticleBranchName, "ChiSq_Timing_KinFit"), locNeutralHypo->Get_ChiSq(), locComboIndex);
}

//KINFIT
if(locKinFitType != d_NoFit)
{
//KINEMATICS
if(locKinFitType != d_P4Fit)
{
	DVector3 locPosition = locNeutralHypo->position();
	TLorentzVector locX4_KinFit(locPosition.X(), locPosition.Y(), locPosition.Z(), locNeutralHypo->time());
	locTreeFillData->Fill_Array<TLorentzVector>(Build_BranchName(locParticleBranchName, "X4_KinFit"), locX4_KinFit, locComboIndex);
}

//update even if vertex-only fit, because neutral momentum defined by vertex
DLorentzVector locDP4 = locNeutralHypo->lorentzMomentum();
TLorentzVector locP4_KinFit(locDP4.Px(), locDP4.Py(), locDP4.Pz(), locDP4.E());
locTreeFillData->Fill_Array<TLorentzVector>(Build_BranchName(locParticleBranchName, "P4_KinFit"), locP4_KinFit, locComboIndex);
}
2239}

2241void DEventWriterROOT::Fill_KinFitData(DTreeFillData* locTreeFillData, JEventLoop* locEventLoop, const DReaction* locReaction, const DMCReaction* locMCReaction, const vector<const DMCThrown*>& locMCThrowns,
const DMCThrownMatching* locMCThrownMatching, const DDetectorMatches* locDetectorMatches,
const vector<const DBeamPhoton*>& locBeamPhotons, const vector<const DChargedTrackHypothesis*>& locChargedHypos,
const vector<const DNeutralParticleHypothesis*>& locNeutralHypos, const deque<const DParticleCombo*>& locParticleCombos) const
2245{
if(!STORE_PULL_INFO && !STORE_ERROR_MATRIX_INFO)
return;

if(STORE_ERROR_MATRIX_INFO){
 locTreeFillData->Fill_Single<Int_t>("numEntries_ParticleErrM", nEntriesParticleCov);
 locTreeFillData->Fill_Single<Int_t>("numEntries_ShowerErrM", nEntriesShowerCov);
}
  bool writeOutPulls = getPullFlag( locReaction );
      
  //Retreive and fill pull-information:
  //######################################################################################################################
  DKinFitType kfitType = locReaction->Get_KinFitType();
  
  //Now fill the pulls for all particles:
  int nParticleCombos = locParticleCombos.size();
  int nMeasuredFinalParticles;
  double currentCharge;
  Particle_t currentPID;
      
  //----------------------------------------------------------------------------
  for(int iPC=0;iPC<nParticleCombos;iPC++){
      Get_PullsFromFit(locParticleCombos.at(iPC));
      
      vector<string> pidVector = getPIDinReaction(locReaction);
      map<Particle_t,int> assignMap = getNameForParticle(locReaction,pidVector);
      
      
      //Look at beam photons:
      //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      const DParticleComboStep* locParticleComboStep = locParticleCombos[iPC]->Get_ParticleComboStep(0);
      const DKinematicData* beamPhoton = locParticleComboStep->Get_InitialParticle_Measured();
      
      if(kfitType == d_P4Fit || kfitType == d_P4AndVertexFit){
  if(writeOutPulls){	
           map<DKinFitPullType, double> someBeamMap = getPulls(beamPhoton,NULL__null,kfitType);
           fillTreePullBranches(locTreeFillData,"ComboBeam",kfitType,someBeamMap,iPC,true);
           someBeamMap.clear();
     }
   }
      //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      
      //Look at the decay products:
      //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      vector<const DKinematicData*> finalParticlesMeasured = locParticleCombos.at(iPC)->Get_FinalParticles_Measured(locReaction, d_AllCharges);
      vector<const DKinematicData*> finalParticles = locParticleCombos.at(iPC)->Get_FinalParticles(locReaction, d_AllCharges);
      
      vector<const JObject*> finalParticleObjects = locParticleCombos.at(iPC)->Get_FinalParticle_SourceObjects(d_AllCharges);
      nMeasuredFinalParticles = finalParticlesMeasured.size();
      
      //----------------------------------------------------------------------------
      for(int iFP=0;iFP<nMeasuredFinalParticles;iFP++){
          const DKinematicData* part = finalParticlesMeasured.at(iFP);
          const DKinematicData* partFit = finalParticles.at(iFP);
          const  DNeutralShower* sh = dynamic_cast<const DNeutralShower*>(finalParticleObjects.at(iFP));
          bool isNeutral = false;
          
          map<DKinFitPullType, double>  someMap;
          if(writeOutPulls) someMap = getPulls(part,sh,kfitType);
          //-----------------------------------------------------
          if(sh == NULL__null){
		currentCharge = part->charge();
              currentPID = part->PID();
	}else{
		currentCharge = 0.0;
              currentPID = Gamma;
    }
          
          if(currentCharge == 0) isNeutral = true;
          //-----------------------------------------------------
          
          string branchName = assignName(currentPID,assignMap);
          assignMap[currentPID]--;
          //-----------------------------------------------------
          if(branchName != "nada"){
              if(STORE_PULL_INFO){
                  fillTreePullBranches(locTreeFillData,branchName,kfitType,someMap,iPC,isNeutral);
                  fillTreeTrackPullBranches(locTreeFillData,branchName,kfitType,someMap,iPC,isNeutral,part,partFit);
              }
              if(STORE_ERROR_MATRIX_INFO)
              	fillTreeErrMBranches(locTreeFillData,branchName,kfitType,part,sh,isNeutral);
          }
          //-----------------------------------------------------
          if(STORE_PULL_INFO) someMap.clear();
  }
      //----------------------------------------------------------------------------
      //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      
      myPullsMap.clear();
      assignMap.clear();
      abundanceMap.clear();
      
  }

2339}


2342//I guess this is the most clumsiest way of assigning a few names to a tree-branch, but it is doing what it is supposed to do...
2343//*************************************************************************************************************************************
2344vector<string> DEventWriterROOT::getPIDinReaction(const DReaction* locReaction)const
2345{
  vector <string> myVector;
  vector<Particle_t> particle = locReaction->Get_FinalPIDs(-1,false,false,d_AllCharges,true);
  int nParticleIDs = particle.size();
  //----------------------------------------------------
  for(int k=0;k<nParticleIDs;k++){
      myVector.push_back(EnumString(particle.at(k)));
  }
  //----------------------------------------------------
  return myVector;
2355}

2357map<Particle_t, int> DEventWriterROOT::getNameForParticle(const DReaction* someReaction, vector<string> someVector)const
2358{
  map <Particle_t,int> myMap;
  vector<Particle_t> particle = someReaction->Get_FinalPIDs(-1,false,false,d_AllCharges,true);
  int nParticleIDs = particle.size();
  int currentCounter = 0;
  int nEl = someVector.size();
  //----------------------------------------------------
  for(int k=0;k<nParticleIDs;k++){
      currentCounter = 0;
      
      //----------------------------------------------------
      for(int j=0;j<nEl;j++){
          if(EnumString(particle.at(k)) == someVector.at(j)){
              currentCounter++;
              someVector.at(j) = "";
          }
      }
      //----------------------------------------------------
      if(currentCounter > 0){
          abundanceMap[particle.at(k)] = currentCounter;
          myMap[particle.at(k)] = currentCounter;
      }
  }
  //----------------------------------------------------
  return myMap;
2383}

2385string DEventWriterROOT::assignName(Particle_t someParticle, map<Particle_t,int> someMap)const
2386{
  string myName = "nada";
  int abundance = abundanceMap.find(someParticle)->second;
  if(abundance == 1) {
      myName = EnumString(someParticle);
  } else if(abundance > 1) {
      int index = abundance - someMap.find(someParticle)->second +1;
      string addName;
      ostringstream convert;
      convert << index;
      addName = convert.str();
      if(strcmp(EnumString(someParticle),"Gamma")==0) {
          myName = "Photon" + addName;
      } else {
      	myName = EnumString(someParticle) + addName;
}
  }
  
  return myName;
2405}

2407void DEventWriterROOT::Get_PullsFromFit(const DParticleCombo* particleCombos) const
2408{
  const DKinFitResults* fitResults = particleCombos->Get_KinFitResults();
  if(fitResults != NULL__null){
      fitResults->Get_Pulls(myPullsMap);
  }
2413}

2415map<DKinFitPullType, double> DEventWriterROOT::getPulls(const JObject* particle, const JObject* shower, DKinFitType yourFitType) const{
map<DKinFitPullType, double> myMap;

if(particle != NULL__null || shower != NULL__null) {

   if(yourFitType == d_P4Fit || yourFitType == d_VertexFit) {
   myMap = myPullsMap.find(particle)->second;
  }else if(yourFitType == d_P4AndVertexFit ){ //Option noFit is not considered
      if(shower == NULL__null) {
	   myMap = myPullsMap.find(particle)->second;
   } else myMap = myPullsMap.find(shower)->second; 	
  }
  }

return myMap;
2430}


2433//Get the pull-features:
2434void DEventWriterROOT::setTreePullBranches(DTreeBranchRegister& locBranchRegister,string yourBranchName,DKinFitType yourFitType, int yourNCombos, bool isNeutral) const
2435{
  string locArraySizeString = "NumCombos";
  
  if(yourFitType == d_P4Fit){
      //Regular Pulls:
      locBranchRegister.Register_FundamentalArray<Double_t>(Build_BranchName(yourBranchName,"Px_Pull"),locArraySizeString, yourNCombos);
      locBranchRegister.Register_FundamentalArray<Double_t>(Build_BranchName(yourBranchName,"Py_Pull"),locArraySizeString, yourNCombos);
      locBranchRegister.Register_FundamentalArray<Double_t>(Build_BranchName(yourBranchName,"Pz_Pull"),locArraySizeString, yourNCombos);
      
      //Pulls used for tracking:
      locBranchRegister.Register_FundamentalArray<Double_t>(Build_BranchName(yourBranchName,"QPt_Pull"),locArraySizeString, yourNCombos);
      locBranchRegister.Register_FundamentalArray<Double_t>(Build_BranchName(yourBranchName,"Phi_Pull"),locArraySizeString, yourNCombos);
      locBranchRegister.Register_FundamentalArray<Double_t>(Build_BranchName(yourBranchName,"tanLambda_Pull"),locArraySizeString, yourNCombos);
locBranchRegister.Register_FundamentalArray<Double_t>(Build_BranchName(yourBranchName,"D_Pull"),locArraySizeString, yourNCombos);
    
  }else if(yourFitType == d_P4AndVertexFit){
      //------------------------------------------------------------
      if(isNeutral){
          locBranchRegister.Register_FundamentalArray<Double_t>(Build_BranchName(yourBranchName,"E_Pull"),locArraySizeString, yourNCombos);
          locBranchRegister.Register_FundamentalArray<Double_t>(Build_BranchName(yourBranchName,"Xx_Pull"),locArraySizeString, yourNCombos);
          locBranchRegister.Register_FundamentalArray<Double_t>(Build_BranchName(yourBranchName,"Xy_Pull"),locArraySizeString, yourNCombos);
          locBranchRegister.Register_FundamentalArray<Double_t>(Build_BranchName(yourBranchName,"Xz_Pull"),locArraySizeString, yourNCombos);
      }else{
          //Regular Pulls:
          locBranchRegister.Register_FundamentalArray<Double_t>(Build_BranchName(yourBranchName,"Px_Pull"),locArraySizeString, yourNCombos);
          locBranchRegister.Register_FundamentalArray<Double_t>(Build_BranchName(yourBranchName,"Py_Pull"),locArraySizeString, yourNCombos);
          locBranchRegister.Register_FundamentalArray<Double_t>(Build_BranchName(yourBranchName,"Pz_Pull"),locArraySizeString, yourNCombos);
          locBranchRegister.Register_FundamentalArray<Double_t>(Build_BranchName(yourBranchName,"Xx_Pull"),locArraySizeString, yourNCombos);
          locBranchRegister.Register_FundamentalArray<Double_t>(Build_BranchName(yourBranchName,"Xy_Pull"),locArraySizeString, yourNCombos);
          locBranchRegister.Register_FundamentalArray<Double_t>(Build_BranchName(yourBranchName,"Xz_Pull"),locArraySizeString, yourNCombos);
          
          //Pulls used for tracking:
          locBranchRegister.Register_FundamentalArray<Double_t>(Build_BranchName(yourBranchName,"QPt_Pull"),locArraySizeString, yourNCombos);
          locBranchRegister.Register_FundamentalArray<Double_t>(Build_BranchName(yourBranchName,"Phi_Pull"),locArraySizeString, yourNCombos);
          locBranchRegister.Register_FundamentalArray<Double_t>(Build_BranchName(yourBranchName,"tanLambda_Pull"),locArraySizeString, yourNCombos);
   locBranchRegister.Register_FundamentalArray<Double_t>(Build_BranchName(yourBranchName,"D_Pull"),locArraySizeString, yourNCombos);
      }
      //------------------------------------------------------------
  }else if(yourFitType == d_VertexFit && !isNeutral){
          locBranchRegister.Register_FundamentalArray<Double_t>(Build_BranchName(yourBranchName,"Xx_Pull"),locArraySizeString, yourNCombos);
          locBranchRegister.Register_FundamentalArray<Double_t>(Build_BranchName(yourBranchName,"Xy_Pull"),locArraySizeString, yourNCombos);
          locBranchRegister.Register_FundamentalArray<Double_t>(Build_BranchName(yourBranchName,"Xz_Pull"),locArraySizeString, yourNCombos);
  }
2478}

2480//Fill the tree branches:
2481void DEventWriterROOT::fillTreePullBranches(DTreeFillData* locTreeFillData,string yourBranchName,DKinFitType yourFitType,map<DKinFitPullType, double> yourPullsMap, int yourIndex, bool isNeutral) const
2482{
  
      if(yourFitType == d_P4Fit){
          locTreeFillData->Fill_Array<Double_t>(Build_BranchName(yourBranchName, "Px_Pull"), yourPullsMap.find(d_PxPull)->second,yourIndex);
          locTreeFillData->Fill_Array<Double_t>(Build_BranchName(yourBranchName, "Py_Pull"), yourPullsMap.find(d_PyPull)->second,yourIndex);
          locTreeFillData->Fill_Array<Double_t>(Build_BranchName(yourBranchName, "Pz_Pull"), yourPullsMap.find(d_PzPull)->second,yourIndex);
  
      }else if(yourFitType == d_P4AndVertexFit){
         //------------------------------------------------------------
         if(isNeutral){
              locTreeFillData->Fill_Array<Double_t>(Build_BranchName(yourBranchName, "E_Pull"), yourPullsMap.find(d_EPull)->second,yourIndex);
              locTreeFillData->Fill_Array<Double_t>(Build_BranchName(yourBranchName, "Xx_Pull"), yourPullsMap.find(d_XxPull)->second,yourIndex);
              locTreeFillData->Fill_Array<Double_t>(Build_BranchName(yourBranchName, "Xy_Pull"), yourPullsMap.find(d_XyPull)->second,yourIndex);
              locTreeFillData->Fill_Array<Double_t>(Build_BranchName(yourBranchName, "Xz_Pull"), yourPullsMap.find(d_XzPull)->second,yourIndex);
         }else{
              locTreeFillData->Fill_Array<Double_t>(Build_BranchName(yourBranchName, "Px_Pull"), yourPullsMap.find(d_PxPull)->second,yourIndex);
              locTreeFillData->Fill_Array<Double_t>(Build_BranchName(yourBranchName, "Py_Pull"), yourPullsMap.find(d_PyPull)->second,yourIndex);
              locTreeFillData->Fill_Array<Double_t>(Build_BranchName(yourBranchName, "Pz_Pull"), yourPullsMap.find(d_PzPull)->second,yourIndex);
              locTreeFillData->Fill_Array<Double_t>(Build_BranchName(yourBranchName, "Xx_Pull"), yourPullsMap.find(d_XxPull)->second,yourIndex);
              locTreeFillData->Fill_Array<Double_t>(Build_BranchName(yourBranchName, "Xy_Pull"), yourPullsMap.find(d_XyPull)->second,yourIndex);
              locTreeFillData->Fill_Array<Double_t>(Build_BranchName(yourBranchName, "Xz_Pull"), yourPullsMap.find(d_XzPull)->second,yourIndex);
         }
         //------------------------------------------------------------
      }else if(yourFitType == d_VertexFit && !isNeutral){
             locTreeFillData->Fill_Array<Double_t>(Build_BranchName(yourBranchName, "Xx_Pull"), yourPullsMap.find(d_XxPull)->second,yourIndex);
             locTreeFillData->Fill_Array<Double_t>(Build_BranchName(yourBranchName, "Xy_Pull"), yourPullsMap.find(d_XyPull)->second,yourIndex);
             locTreeFillData->Fill_Array<Double_t>(Build_BranchName(yourBranchName, "Xz_Pull"), yourPullsMap.find(d_XzPull)->second,yourIndex);
      }
2510}

2512//Make sure, the pull information is written out, when available
2513void DEventWriterROOT::setPullFlag(const DReaction* currentReaction, bool myFlag) const
2514{
  writePulls[currentReaction] = myFlag;
2516}

2518bool DEventWriterROOT::getPullFlag(const DReaction* currentReaction) const
2519{
  bool outFlag = writePulls.find(currentReaction)->second;
  return outFlag;
2522}

2524//Now fill the individual elements of the error matrix to the tree:
2525void DEventWriterROOT::fillTreeErrMBranches(DTreeFillData* locTreeFillData,string yourBranchName,DKinFitType yourFitType, const DKinematicData* particle, const DNeutralShower* shower, bool isNeutral) const
2526{
  if(yourFitType == d_P4Fit) {
fillTreeParticleErrMBranches(locTreeFillData,yourBranchName,particle);
} else if(yourFitType == d_P4AndVertexFit) {
if(isNeutral) {
	fillTreeShowerErrMBranches(locTreeFillData,yourBranchName,shower);
} else {
	fillTreeParticleErrMBranches(locTreeFillData,yourBranchName,particle);
}
} else if(yourFitType == d_VertexFit && !isNeutral) {
fillTreeParticleErrMBranches(locTreeFillData,yourBranchName,particle);
}     	
2538}

2540//Fill the 7x7 error matrix elements:
2541void DEventWriterROOT::fillTreeParticleErrMBranches(DTreeFillData* locTreeFillData,string yourBranchName, const DKinematicData* particle) const
2542{
if(particle != NULL__null){	
  const TMatrixTSym<float> errMatrix = *( particle->errorMatrix() );

  int mIndex = 0;
  float matrixEl = 0.0;
  for(int row=0;row<7;row++){
   for(int col=0;col<7;col++){
	    if(col >= row){
			matrixEl = errMatrix[row][col];
			locTreeFillData->Fill_Array<Float_t>(Build_BranchName(yourBranchName, "ErrMatrix"),matrixEl,mIndex);
			mIndex++;
		}
   }
  }

  }
2559}

2561//Fill the 5x5 error matrix elements:
2562void DEventWriterROOT::fillTreeShowerErrMBranches(DTreeFillData* locTreeFillData,string yourBranchName, const DNeutralShower* shower) const
2563{
if(shower != NULL__null){
  const TMatrixTSym <float> errMatrix = *( shower->dCovarianceMatrix );
  float matrixEl = 0.0;
  
  int mIndex = 0;
  for(int row=0;row<5;row++){
   for(int col=0;col<5;col++){
	    if(col >= row){
		   matrixEl = errMatrix[row][col];
	       locTreeFillData->Fill_Array<Float_t>(Build_BranchName(yourBranchName, "ErrMatrix"),matrixEl,mIndex);
		   mIndex++;
		}
   }
  }

  }
2580}

2582//GET THE TRACKING PULLS:

2584//Get the difference in pull errors for one momentum component:
2585double DEventWriterROOT::getDiffSquaredErrP_Component(double yourRecComponent,double yourFitComponent,double yourPull) const
2586{
  double diff = yourRecComponent - yourFitComponent;
  double error_component = 0.0;
  
  if(diff != 0.0 && yourPull != 0.0){
      error_component = diff*diff / (yourPull*yourPull);
  }
  
  return error_component;
2595}

2597// Get the difference in pull errors for the position
2598vector< vector<double> > DEventWriterROOT::getSquaredErrX(const DKinematicData* particle, const DKinematicData* particleFit, map<DKinFitPullType, double> yourPullsMap) const{
  
  double x_pull = yourPullsMap.find(d_XxPull)->second;
  double y_pull = yourPullsMap.find(d_XyPull)->second;
  
  double x_rec = particle->x();
  double x_fit = particleFit->x();
  
  double y_rec = particle->y();
  double y_fit = particleFit->y();
  
  double error_x = getDiffSquaredErrP_Component(x_rec,x_fit,x_pull);
  double error_y = getDiffSquaredErrP_Component(y_rec,y_fit,y_pull);
  
  const TMatrixTSym<float> errMatrix = *( particle->errorMatrix() ); //--> Reconstructed errors!
  double error_x_rec = errMatrix[3][3];
  double error_y_rec = errMatrix[4][4];
  
  double error_x_fit = error_x_rec - error_x;
  double error_y_fit = error_y_rec - error_y;
  
  vector <double> xComponent;
  vector <double> yComponent;
  
  xComponent.push_back(error_x_rec);
  xComponent.push_back(error_x_fit);
  
  yComponent.push_back(error_y_rec);
  yComponent.push_back(error_y_fit);
  
  vector < vector<double> > myErrors;
  myErrors.push_back(xComponent);
  myErrors.push_back(yComponent);
 
  return myErrors;

2634}

2636//Get the difference in pull errors for the momentum
2637vector< vector<double> > DEventWriterROOT::getSquaredErrP(const DKinematicData* particle, const DKinematicData* particleFit, map<DKinFitPullType, double> yourPullsMap) const
2638{
  
  double px_pull = yourPullsMap.find(d_PxPull)->second;
  double py_pull = yourPullsMap.find(d_PyPull)->second;
  double pz_pull = yourPullsMap.find(d_PzPull)->second;
  
  double px_rec = particle->px();
  double px_fit = particleFit->px();
  
  double py_rec = particle->py();
  double py_fit = particleFit->py();
  
  double pz_rec = particle->pz();
  double pz_fit = particleFit->pz();
  
  double error_x = getDiffSquaredErrP_Component(px_rec,px_fit,px_pull);
  double error_y = getDiffSquaredErrP_Component(py_rec,py_fit,py_pull);
  double error_z = getDiffSquaredErrP_Component(pz_rec,pz_fit,pz_pull);
  
  const TMatrixTSym<float> errMatrix = *( particle->errorMatrix() ); //--> Reconstructed errors!
  double error_x_rec = errMatrix[0][0];
  double error_y_rec = errMatrix[1][1];
  double error_z_rec = errMatrix[2][2];
  
  double error_x_fit = error_x_rec - error_x;
  double error_y_fit = error_y_rec - error_y;
  double error_z_fit = error_z_rec - error_z;
  
  vector <double> xComponent;
  vector <double> yComponent;
  vector <double> zComponent;
  
  xComponent.push_back(error_x_rec);
  xComponent.push_back(error_x_fit);
  
  yComponent.push_back(error_y_rec);
  yComponent.push_back(error_y_fit);
  
  zComponent.push_back(error_z_rec);
  zComponent.push_back(error_z_fit);
  
  vector < vector<double> > myErrors;
  myErrors.push_back(xComponent);
  myErrors.push_back(yComponent);
  myErrors.push_back(zComponent);
 
  return myErrors;
2685}

2687//1.) Start with phi:
2688//Calculate the error for one state (fitted,reconstructed):
2689double DEventWriterROOT::getPhiError(const DKinematicData* particle, vector< vector<double> > yourErrorMatrix, int isFitted) const
2690{
  double myPhiError = 0.0;
  double myPx = particle->px();
  double myPy = particle->py();
  
  //Phi = artan(Px/Py) --> caclulate derivatives from this...
  if(myPy != 0.0){
      
      double myDPx = yourErrorMatrix[0][isFitted];
      double myDPy = yourErrorMatrix[1][isFitted];
      
      double arg = myPx / myPy;
      double factor1 = 1.0 / pow( myPy*(1.0 + arg*arg), 2);
      double factor2 = (myDPx + arg*arg*myDPy);
      
      myPhiError = factor1*factor2;
  }
  
  return myPhiError;
2709}


2712//Now calculate the pull:
2713double DEventWriterROOT::getPhiPull(const DKinematicData* particle, const DKinematicData* particleFit, vector< vector<double> > yourErrorMatrix) const
2714{
  double myPhiPull = -100.0;
  
  double phi_rec = ( particle->momentum() ).Phi();
  double phi_fit = ( particleFit->momentum() ).Phi();
  
  double dPhi_rec = getPhiError(particle,yourErrorMatrix,0);
  double dPhi_fit = getPhiError(particleFit,yourErrorMatrix,1);
  
  double pull_norm = dPhi_rec - dPhi_fit;
  
  if(pull_norm > 0.0){
      myPhiPull = (phi_rec - phi_fit) / sqrt(pull_norm);
  }
  
  return myPhiPull;
2730}

2732//2.) Tan(lambda):
2733//Get the error of lambda:
2734double DEventWriterROOT::getLambdaError(const DKinematicData* particle, vector< vector<double> > yourErrorMatrix, int isFitted) const
2735{
  double myTanLambdaError = 0.0;
  //lambda = pi/2 - theta --> dlambda = dtheta...
  
  double myPx = particle->px();
  double myPy = particle->py();
  double myPz = particle->pz();
  
  double myDPx = yourErrorMatrix[0][isFitted];
  double myDPy = yourErrorMatrix[1][isFitted];
  double myDPz = yourErrorMatrix[2][isFitted];
  
  double r2 = myPx*myPx + myPy*myPy + myPz*myPz;
  double lambda = 0.5*M_PI3.14159265358979323846 - ( particle->momentum() ).Theta();
  //---------------------------
  if(r2 > 0.0 && cos(lambda) != 0.0){
      double r = sqrt(r2);
      double factor1 = 1.0 / ( pow(r,6) * pow(cos(lambda),4) );
      double sum1 = myPx*myPz;
      double sum2 = myPy*myPz;
      double sum3 = (r - myPz*myPz);
      double factor2 = sum1*sum1*myDPx + sum2*sum2*myDPy + sum3*sum3*myDPz;
      
      myTanLambdaError = factor1*factor2;
  }
  //---------------------------
  
  
  return myTanLambdaError;
2764}


2767//Go for the pull itself:
2768double DEventWriterROOT::getTanLambdaPull(const DKinematicData* particle, const DKinematicData* particleFit, vector< vector<double> > yourErrorMatrix) const
2769{
  double myTanLambdaPull = -100.0;
  
  //Reconstructed values:
  double lambda_rec = 0.5*M_PI3.14159265358979323846 - ( particle->momentum() ).Theta();
  double tanLambda_rec = tan(lambda_rec);
  double dtanLambda_rec = getLambdaError(particle,yourErrorMatrix,0);
  
  //Fitted values:
  double lambda_fit = 0.5*M_PI3.14159265358979323846 - ( particleFit->momentum() ).Theta();
  double tanLambda_fit = tan(lambda_fit);
  double dtanLambda_fit = getLambdaError(particleFit,yourErrorMatrix,1);
  
  double pull_norm = dtanLambda_rec-dtanLambda_fit;
  //---------------------------
  if(pull_norm > 0.0){
      myTanLambdaPull = (tanLambda_rec-tanLambda_fit) / sqrt(pull_norm);
  }
  //---------------------------
  
  return myTanLambdaPull;
2790}

2792//3.) q/pt:
2793//Calculate error for 1/pt:
2794double DEventWriterROOT::getQPTError(const DKinematicData* particle, vector< vector<double> > yourErrorMatrix, int isFitted) const
2795{
  double myQPTError = 0.0;
  double myPx = particle->px();
  double myPy = particle->py();
  double r2 = myPx*myPx + myPy*myPy;
  
  double myDPx = yourErrorMatrix[0][isFitted];
  double myDPy = yourErrorMatrix[1][isFitted];
  //---------------------------
  if(r2 > 0.0){
      double q = particle->charge();
      double factor1 = q*q / pow(r2,3);
      double factor2 = (myPx*myPx*myDPx + myPy*myPy*myDPy);
      
      myQPTError = factor1*factor2;
  }
  //---------------------------
  
  return myQPTError;
2814}


2817//Get the pull:
2818double DEventWriterROOT::getQPTPull(const DKinematicData* particle, const DKinematicData* particleFit, vector< vector<double> > yourErrorMatrix) const
2819{
  double myPTPull = -100.0;
  
  double particle_charge = particle->charge();
  
  double myPx_rec = particle->px();
  double myPy_rec = particle->py();
  double r_rec = myPx_rec*myPx_rec + myPy_rec*myPy_rec;
  
  double myPx_fit = particleFit->px();
  double myPy_fit = particleFit->py();
  double r_fit = myPx_fit*myPx_fit + myPy_fit*myPy_fit;
  
  //---------------------------
  if(r_rec > 0.0 && r_fit > 0.0){
      double qpt_rec = particle_charge/sqrt(r_rec);
      double dqpt_rec = getQPTError(particle,yourErrorMatrix,0);
      
      double qpt_fit = particle_charge/sqrt(r_fit);
      double dqpt_fit = getQPTError(particleFit,yourErrorMatrix,1);
      
      double pull_norm = dqpt_rec-dqpt_fit;
      //---------------------------
      if(pull_norm > 0.0){
 myPTPull = (qpt_rec-qpt_fit)/sqrt(pull_norm);
      }
      //---------------------------
  }
  //---------------------------
  
  
  return myPTPull;
2851}

2853//4.) D:
2854//Get the error of D:
2855double DEventWriterROOT::getDError(const DKinematicData* particle, vector< vector<double> > yourErrorMatrix, int isFitted) const
2856{

double myDsq=particle->position().Perp2()+1e-8;
double myx=particle->x();
double myy=particle->y();

double myVarX=yourErrorMatrix[0][isFitted];
double myVarY=yourErrorMatrix[1][isFitted];
double myDError=(myx*myx*myVarX+myy*myy*myVarY)/myDsq;

return myDError;
2867}

2869//Go for the pull itself:
2870double DEventWriterROOT::getDPull(const DKinematicData* particle, const DKinematicData* particleFit, vector< vector<double> > yourErrorMatrix) const
2871{
double myDPull=-100.0;

//Reconstructed values:
double D_rec=particle->position().Perp();
double dx_rec=particle->x();
double dy_rec=particle->y();
double phi_rec=particle->momentum().Phi();
double cosphi_rec=cos(phi_rec);
double sinphi_rec=sin(phi_rec);
if ((dx_rec>0.0 && sinphi_rec>0.0) || (dy_rec<0.0 && cosphi_rec>0.0)
    || (dy_rec>0.0 && cosphi_rec<0.0) || (dx_rec<0.0 && sinphi_rec<0.0))
  D_rec*=-1.;
double dD_rec=getDError(particle,yourErrorMatrix,0);

//Fitted values:
double D_fit=particleFit->position().Perp();
double dx_fit=particleFit->x();
double dy_fit=particleFit->y();
double phi_fit=particleFit->momentum().Phi();
double cosphi_fit=cos(phi_fit);
double sinphi_fit=sin(phi_fit);
if ((dx_fit>0.0 && sinphi_fit>0.0) || (dy_fit<0.0 && cosphi_fit>0.0)
    || (dy_fit>0.0 && cosphi_fit<0.0) || (dx_fit<0.0 && sinphi_fit<0.0))
  D_fit*=-1.;
double dD_fit=getDError(particle,yourErrorMatrix,1);

double pull_norm = dD_rec-dD_fit;
//---------------------------
if(pull_norm > 0.0){
  myDPull = (D_rec-D_fit) / sqrt(pull_norm);
}
//---------------------------
  
return myDPull;
2906}

2908//Create a vector containing all pulls:
2909vector<double> DEventWriterROOT::collectTrackingPulls(const DKinematicData* particle, const DKinematicData* particleFit, map<DKinFitPullType, double> yourPullsMap) const
2910{
  //1.) Get the pre-and post-fit errors:
  vector< vector<double> > myErrors = getSquaredErrP(particle,particleFit,yourPullsMap);
  vector< vector<double> > myXErrors = getSquaredErrX(particle,particleFit,yourPullsMap);
  
  //2.) Calculate the pulls:
  //2.1.) q/pt:
  double qpt_pull = getQPTPull(particle,particleFit,myErrors);
  //2.2.) phi:
  double phi_pull = getPhiPull(particle,particleFit,myErrors);
  //2.3.) tan(lambda):
  double tanLambda_pull = getTanLambdaPull(particle,particleFit,myErrors);
  //2.4.) D
  double D_pull = getDPull(particle,particleFit,myXErrors);
  
  //3.) Store everything into a vector:
  vector<double> myTrackingPulls;
  myTrackingPulls.push_back(qpt_pull);
  myTrackingPulls.push_back(phi_pull);
  myTrackingPulls.push_back(tanLambda_pull);
  myTrackingPulls.push_back(D_pull);
  
  return myTrackingPulls;
2933}

2935//Fill the tracking pulls:
2936void DEventWriterROOT::fillTreeTrackPullBranches(DTreeFillData* locTreeFillData,string yourBranchName,DKinFitType yourFitType,map<DKinFitPullType, double> yourPullsMap, int yourIndex, bool isNeutral, const DKinematicData* particle, const DKinematicData* particleFit) const
2937{
  if( !isNeutral && (yourFitType == d_P4Fit || yourFitType == d_P4AndVertexFit) ){//Get the tracking pulls for charged tracks only...
      //1.) Get the vector containing the tracking pulls:
      vector<double> myTrackingPulls = collectTrackingPulls(particle,particleFit,yourPullsMap);
      
      //2.) Fill the pulls to the tree:
      locTreeFillData->Fill_Array<Double_t>(Build_BranchName(yourBranchName, "QPt_Pull"),myTrackingPulls[0],yourIndex);
      locTreeFillData->Fill_Array<Double_t>(Build_BranchName(yourBranchName, "Phi_Pull"),myTrackingPulls[1],yourIndex);
      locTreeFillData->Fill_Array<Double_t>(Build_BranchName(yourBranchName, "tanLambda_Pull"),myTrackingPulls[2],yourIndex);
locTreeFillData->Fill_Array<Double_t>(Build_BranchName(yourBranchName, "D_Pull"),myTrackingPulls[3],yourIndex);
  }
2948}


←

/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//
7 
8#ifndef _JEventLoop_
9#define _JEventLoop_
10 
11#include <sys/time.h>
12 
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;
25 
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>
36 
37// The following is here just so we can use ROOT's THtml class to generate documentation.
38#include "cint.h"
39 
40 
41// Place everything in JANA namespace
42namespace jana{
43 
44 
45template<class T> class JFactory;
46class JApplication;
47class JEventProcessor;
48 
49 
50class JEventLoop{
51	public:
52	
53		friend class JApplication;
54	
55		enum data_source_t{
56			DATA_NOT_AVAILABLE = 1,
57			DATA_FROM_CACHE,
58			DATA_FROM_SOURCE,
59			DATA_FROM_FACTORY
60		};
61		
62		typedef struct{
63			string caller_name;
64			string caller_tag;
65			string callee_name;
66			string callee_tag;
67			double start_time;
68			double end_time;
69			data_source_t data_source;
70		}call_stack_t;
71		
72		typedef struct{
73			const char* factory_name;
74			string tag;
75			const char* filename;
76			int line;
77		}error_call_stack_t;
78 
79	                                   JEventLoop(JApplication *app); ///< Constructor
80							   virtual ~JEventLoop(); ////< Destructor
81				   virtual const char* className(void){return static_className();}
82					static const char* static_className(void){return "JEventLoop";}
83		
84						 JApplication* GetJApplication(void) const {return app;} ///< Get pointer to the JApplication object
85                                  void RefreshProcessorListFromJApplication(void); ///< Re-copy the list of JEventProcessors from JApplication
86                      virtual jerror_t AddFactory(JFactory_base* factory); ///< Add a factory
87							  jerror_t RemoveFactory(JFactory_base* factory); ///< Remove a factory
88                        JFactory_base* GetFactory(const string data_name, const char *tag="", bool allow_deftag=true); ///< Get a specific factory pointer
89                vector<JFactory_base*> GetFactories(void) const {return factories;} ///< Get all factory pointers
90                                  void GetFactoryNames(vector<string> &factorynames); ///< Get names of all factories in name:tag format
91                                  void GetFactoryNames(map<string,string> &factorynames); ///< Get names of all factories in map with key=name, value=tag
92                    map<string,string> GetDefaultTags(void) const {return default_tags;}
93                              jerror_t ClearFactories(void); ///< Reset all factories in preparation for next event.
94							  jerror_t PrintFactories(int sparsify=0); ///< Print a list of all factories.
95                              jerror_t Print(const string data_name, const char *tag=""); ///< Print the data of the given type
96 
97						 JCalibration* GetJCalibration();
98                template<class T> bool GetCalib(string namepath, map<string,T> &vals);
99                template<class T> bool GetCalib(string namepath, vector<T> &vals);
100                template<class T> bool GetCalib(string namepath, T &val);
101 
102                            JGeometry* GetJGeometry();
103                template<class T> bool GetGeom(string namepath, map<string,T> &vals);
104			    template<class T> bool GetGeom(string namepath, T &val);
105 
106                     JResourceManager* GetJResourceManager(void);
107                                string GetResource(string namepath);
108                template<class T> bool GetResource(string namepath, T vals, int event_number=0);
109 
110                                  void Initialize(void); ///< Do initializations just before event processing starts
111                              jerror_t Loop(void); ///< Loop over events
112                              jerror_t OneEvent(uint64_t event_number); ///< Process a specific single event (if source supports it)
113                              jerror_t OneEvent(void); ///< Process a single event
114                           inline void Pause(void){pause = 1;} ///< Pause event processing
115                           inline void Resume(void){pause = 0;} ///< Resume event processing
116                           inline void Quit(void){quit = 1;} ///< Clean up and exit the event loop
117                           inline bool GetQuit(void) const {return quit;}
118                                  void QuitProgram(void);
119 
120		                               // Support for random access of events
121		                          bool HasRandomAccess(void);
122		                          void AddToEventQueue(uint64_t event_number){ next_events_to_process.push_back(event_number); }
123		                          void AddToEventQueue(list<uint64_t> &event_numbers) { next_events_to_process.insert(next_events_to_process.end(), event_numbers.begin(), event_numbers.end()); }
124		                list<uint64_t> GetEventQueue(void){ return next_events_to_process; }
125		                          void ClearEventQueue(void){ next_events_to_process.clear(); }
126 
127        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).
128        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)
129        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
130            template<class T> jerror_t GetFromSource(vector<const T*> &t, JFactory_base *factory=NULL__null); ///< Get data object pointers from source.
131                        inline JEvent& GetJEvent(void){return event;} ///< Get pointer to the current JEvent object.
132                           inline void SetJEvent(JEvent *event){this->event = *event;} ///< Set the JEvent pointer.
133                           inline void SetAutoFree(int auto_free){this->auto_free = auto_free;} ///< Set the Auto-Free flag on/off
134                      inline pthread_t GetPThreadID(void) const {return pthread_id;} ///< Get the pthread of the thread to which this JEventLoop belongs
135                                double GetInstantaneousRate(void) const {return rate_instantaneous;} ///< Get the current event processing rate
136                                double GetIntegratedRate(void) const {return rate_integrated;} ///< Get the current event processing rate
137                                double GetLastEventProcessingTime(void) const {return delta_time_single;}
138                          unsigned int GetNevents(void) const {return Nevents;}
139 
140                           inline bool CheckEventBoundary(uint64_t event_numberA, uint64_t event_numberB);
141 
142                           inline bool GetCallStackRecordingStatus(void){ return record_call_stack; }
143                           inline void DisableCallStackRecording(void){ record_call_stack = false; }
144                           inline void EnableCallStackRecording(void){ record_call_stack = true; }
145		                     inline void CallStackStart(JEventLoop::call_stack_t &cs, const string &caller_name, const string &caller_tag, const string callee_name, const string callee_tag);
146		                     inline void CallStackEnd(JEventLoop::call_stack_t &cs);
147           inline vector<call_stack_t> GetCallStack(void){return call_stack;} ///< Get the current factory call stack
148                           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
149                           inline void AddToErrorCallStack(error_call_stack_t &cs){error_call_stack.push_back(cs);} ///< Add layer to the factory call stack
150     inline vector<error_call_stack_t> GetErrorCallStack(void){return error_call_stack;} ///< Get the current factory error call stack
151                                  void PrintErrorCallStack(void); ///< Print the current factory call stack
152 
153                        const JObject* FindByID(JObject::oid_t id); ///< Find a data object by its identifier.
154            template<class T> const T* FindByID(JObject::oid_t id); ///< Find a data object by its type and identifier
155                        JFactory_base* FindOwner(const JObject *t); ///< Find the factory that owns a data object by pointer
156                        JFactory_base* FindOwner(JObject::oid_t id); ///< Find a factory that owns a data object by identifier
157 
158                                       // User defined references
159                template<class T> void SetRef(T *t);    ///< Add a user reference to this JEventLoop (must be a pointer)
160                  template<class T> T* GetRef(void);   ///< Get a user-defined reference of a specific type
161          template<class T> vector<T*> GetRefsT(void); ///< Get all user-defined refrences of a specicif type
162     vector<pair<const char*, void*> > GetRefs(void){ return user_refs; }  ///< Get copy of full list of user-defined references
163                template<class T> void RemoveRef(T *t); ///< Remove user reference from list
164 
165									   // Convenience methods wrapping JEvent methods of same name
166		                      uint64_t GetStatus(void){return event.GetStatus();}
167		                          bool GetStatusBit(uint32_t bit){return event.GetStatusBit(bit);}
168		                          bool SetStatusBit(uint32_t bit, bool val=true){return event.SetStatusBit(bit, val);}
169		                          bool ClearStatusBit(uint32_t bit){return event.ClearStatusBit(bit);}
170		                          void ClearStatus(void){event.ClearStatus();}
171		                          void SetStatusBitDescription(uint32_t bit, string description){event.SetStatusBitDescription(bit, description);}
172		                        string GetStatusBitDescription(uint32_t bit){return event.GetStatusBitDescription(bit);}
173		                          void GetStatusBitDescriptions(map<uint32_t, string> &status_bit_descriptions){return event.GetStatusBitDescriptions(status_bit_descriptions);}
174                                string StatusWordToString(void);
175 
176	private:
177		JEvent event;
178		vector<JFactory_base*> factories;
179		vector<JEventProcessor*> processors;
180		vector<error_call_stack_t> error_call_stack;
181		vector<call_stack_t> call_stack;
182		JApplication *app;
183		JThread *jthread;
184		bool initialized;
185		bool print_parameters_called;
186		int pause;
187		int quit;
188		int auto_free;
189		pthread_t pthread_id;
190		map<string, string> default_tags;
191		vector<pair<string,string> > auto_activated_factories;
192		bool record_call_stack;
193		string caller_name;
194		string caller_tag;
195		vector<uint64_t> event_boundaries;
196		int32_t event_boundaries_run; ///< Run number boundaries were retrieved from (possbily 0)
197		list<uint64_t> next_events_to_process;
198		
199		uint64_t Nevents;			      ///< Total events processed (this thread)
200		uint64_t Nevents_rate;		   ///< Num. events accumulated for "instantaneous" rate
201		double delta_time_single;		///< Time spent processing last event
202		double delta_time_rate;			///< Integrated time accumulated "instantaneous" rate (partial number of events)
203		double delta_time;				///< Total time spent processing events (this thread)
204		double rate_instantaneous;		///< Latest instantaneous rate
205		double rate_integrated;			///< Rate integrated over all events
206   
207		static data_source_t null_data_source;
208 
209		vector<pair<const char*, void*> > user_refs;
210};
211 
212 
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
218 
219#if !defined(__CINT__) && !defined(__CLING__)
220 
221//-------------
222// GetSingle
223//-------------
224template<class T>
225JFactory<T>* JEventLoop::GetSingle(const T* &t, const char *tag, bool exception_if_not_one)
226{
227	/// This is a convenience method that can be used to get a pointer to the single
228	/// object of type T from the specified factory. It simply calls the Get(vector<...>) method
229	/// and copies the first pointer into "t" (or NULL if something other than 1 object is returned).
230	/// 
231	/// This is intended to address the common situation in which there is an interest
232	/// in the event if and only if there is exactly 1 object of type T. If the event
233	/// has no objects of that type or more than 1 object of that type (for the specified
234	/// factory) then an exception of type "unsigned long" is thrown with the value
235	/// being the number of objects of type T. You can supress the exception by setting
236	/// exception_if_not_one to false. In that case, you will have to check if t==NULL to
237	/// know if the call succeeded.
238	vector<const T*> v;
239	JFactory<T> *fac = Get(v, tag);
2
←
Calling 'JEventLoop::Get'→
240 
241	if(v.size()!=1){
242		t = NULL__null;
243		if(exception_if_not_one) throw v.size();
244	}
245	
246	t = v[0];
247	
248	return fac;
249}
250 
251//-------------
252// Get
253//-------------
254template<class T> 
255JFactory<T>* JEventLoop::Get(vector<const T*> &t, const char *tag, bool allow_deftag)
256{
257	/// Retrieve or generate the array of objects of
258	/// type T for the curent event being processed
259	/// by this thread.
260	///
261	/// By default, preference is given to reading the
262	/// objects from the data source(e.g. file) before generating
263	/// them in the factory. A flag exists in the factory
264	/// however to change this so that the factory is
265	/// given preference.
266	///
267	/// Note that regardless of the setting of this flag,
268	/// the data are only either read in or generated once.
269	/// Ownership of the objects will always be with the
270	/// factory so subsequent calls will always return pointers to
271	/// the same data.
272	///
273	/// If the factory is called on to generate the data,
274	/// it is done by calling the factory's Get() method
275	/// which, in turn, calls the evnt() method. 
276	/// 
277	/// First, we just call the GetFromFactory() method.
278	/// It will make the initial decision as to whether
279	/// it should look in the source first or not. If
280	/// it returns NULL, then the factory couldn't be
281	/// found so we automatically try the file.
282	///
283	/// Note that if no factory exists to hold the objects
284	/// from the file, one can be created automatically
285	/// providing the <i>JANA:AUTOFACTORYCREATE</i>
286	/// configuration parameter is set.
287 
288	// Check if a tag was specified for this data type to use for the
289	// default.
290	const char *mytag = tag2.1
'tag' is not equal to NULL
2.1
'tag' is not equal to NULL
2.1
'tag' is not equal to NULL
==NULL__null ? "":tag; // protection against NULL tags
3
←
'?' condition is false→
291	if(strlen(mytag)==0 && allow_deftag3.1
'allow_deftag' is true
3.1
'allow_deftag' is true
3.1
'allow_deftag' is true
){
4
←
Taking true branch→
292		map<string, string>::const_iterator iter = default_tags.find(T::static_className());
293		if(iter!=default_tags.end())tag = iter->second.c_str();
5
←
Assuming the condition is true→
6
←
Taking true branch→
7
←
Value assigned to 'tag'→
294	}
295	
296	
297	// If we are trying to keep track of the call stack then we
298	// need to add a new call_stack_t object to the the list
299	// and initialize it with the start time and caller/callee
300	// info.
301	call_stack_t cs;
302 
303	// Optionally record starting info of call stack entry
304	if(record_call_stack) CallStackStart(cs, caller_name, caller_tag, T::static_className(), tag);
8
←
Assuming field 'record_call_stack' is false→
9
←
Taking false branch→
305 
306	// Get the data (or at least try to)
307	JFactory<T>* factory=NULL__null;
308	try{
309		factory = GetFromFactory(t, tag, cs.data_source, allow_deftag);
10
←
Passing value via 2nd parameter 'tag'→
11
←
Calling 'JEventLoop::GetFromFactory'→
310		if(!factory){
311			// No factory exists for this type and tag. It's possible
312			// that the source may be able to provide the objects
313			// but it will need a place to put them. We can create a
314			// dumb JFactory just to hold the data in case the source
315			// can provide the objects. Before we do though, make sure
316			// the user condones this via the presence of the
317			// "JANA:AUTOFACTORYCREATE" config parameter.
318			string p;
319			try{
320				gPARMS->GetParameter("JANA:AUTOFACTORYCREATE", p);
321			}catch(...){}
322			if(p.size()==0){
323				jout<<std::endl;
324				_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;
325				jout<<"No factory of type \""<<T::static_className()<<"\" with tag \""<<tag<<"\" exists."<<std::endl;
326				jout<<"If you are reading objects from a file, I can auto-create a factory"<<std::endl;
327				jout<<"of the appropriate type to hold the objects, but this feature is turned"<<std::endl;
328				jout<<"off by default. To turn it on, set the \"JANA:AUTOFACTORYCREATE\""<<std::endl;
329				jout<<"configuration parameter. This can usually be done by passing the"<<std::endl;
330				jout<<"following argument to the program from the command line:"<<std::endl;
331				jout<<std::endl;
332				jout<<"   -PJANA:AUTOFACTORYCREATE=1"<<std::endl;
333				jout<<std::endl;
334				jout<<"Note that since the most commonly expected occurance of this situation."<<std::endl;
335				jout<<"is an error, the program will now throw an exception so that the factory."<<std::endl;
336				jout<<"call stack can be printed."<<std::endl;
337				jout<<std::endl;
338				throw exception();
339			}else{
340				AddFactory(new JFactory<T>(tag));
341				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;
342			
343				// Now try once more. The GetFromFactory method will call
344				// GetFromSource since it's empty.
345				factory = GetFromFactory(t, tag, cs.data_source, allow_deftag);
346			}
347		}
348	}catch(exception &e){
349		// Uh-oh, an exception was thrown. Add us to the call stack
350		// and re-throw the exception
351		error_call_stack_t ecs;
352		ecs.factory_name = T::static_className();
353		ecs.tag = tag;
354		ecs.filename = NULL__null;
355		error_call_stack.push_back(ecs);
356		throw e;
357	}
358	
359	// If recording the call stack, update the end_time field and add to stack
360	if(record_call_stack) CallStackEnd(cs);
361	
362	return factory;
363}
364 
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{
371	// We need to find the factory providing data type T with
372	// tag given by "tag". 
373	vector<JFactory_base*>::iterator iter=factories.begin();
374	JFactory<T> *factory = NULL__null;
375	string className(T::static_className());
376	
377	// Check if a tag was specified for this data type to use for the
378	// default.
379	const char *mytag = tag==NULL__null ? "":tag; // protection against NULL tags
12
←
Assuming 'tag' is equal to NULL→
13
←
Assuming pointer value is null→
14
←
'?' condition is true→
380	if(strlen(mytag)==0 && allow_deftag14.1
'allow_deftag' is true
14.1
'allow_deftag' is true
14.1
'allow_deftag' is true
){
15
←
Taking true branch→
381		map<string, string>::const_iterator iter = default_tags.find(className);
382		if(iter!=default_tags.end())tag = iter->second.c_str();
16
←
Assuming the condition is false→
17
←
Taking false branch→
383	}
384	
385	for(; iter!=factories.end(); iter++){
18
←
Calling 'operator!=<jana::JFactory_base **, std::vector<jana::JFactory_base *>>'→
21
←
Returning from 'operator!=<jana::JFactory_base **, std::vector<jana::JFactory_base *>>'→
22
←
Loop condition is true.  Entering loop body→
386		// It turns out a long standing bug in g++ makes dynamic_cast return
387		// zero improperly when used on objects created on one side of
388		// a dynamically shared object (DSO) and the cast occurs on the 
389		// other side. I saw bug reports ranging from 2001 to 2004. I saw
390		// saw it first-hand on LinuxEL4 using g++ 3.4.5. This is too bad
391		// since it is much more elegant (and safe) to use dynamic_cast.
392		// To avoid this problem which can occur with plugins, we check
393		// the name of the data classes are the same. (sigh)
394		//factory = dynamic_cast<JFactory<T> *>(*iter);
395		if(className == (*iter)->GetDataClassName())factory = (JFactory<T>*)(*iter);
23
←
Taking true branch→
396		if(factory == NULL__null)continue;
24
←
Assuming 'factory' is not equal to NULL→
25
←
Taking false branch→
397		const char *factag = factory->Tag()==NULL__null ? "":factory->Tag();
26
←
Assuming the condition is true→
27
←
'?' condition is true→
398		if(!strcmp(factag, tag)){
28
←
Null pointer passed to 2nd parameter expecting 'nonnull'
399			break;
400		}else{
401			factory=NULL__null;
402		}
403	}
404	
405	// If factory not found, just return now
406	if(!factory){
407		data_source = DATA_NOT_AVAILABLE;
408		return NULL__null;
409	}
410	
411	// OK, we found the factory. If the evnt() routine has already
412	// been called, then just call the factory's Get() routine
413	// to return a copy of the existing data
414	if(factory->evnt_was_called()){
415		factory->CopyFrom(t);
416		data_source = DATA_FROM_CACHE;
417		return factory;
418	}
419	
420	// Next option is to get the objects from the data source
421	if(factory->GetCheckSourceFirst()){
422		// If the object type/tag is found in the source, it
423		// will return NOERROR, even if there are zero instances
424		// of it. If it is not available in the source then it
425		// will return OBJECT_NOT_AVAILABLE.
426		
427		jerror_t err = GetFromSource(t, factory);
428		if(err == NOERROR){
429			// A return value of NOERROR means the source had the objects
430			// even if there were zero of them.(If the source had no
431			// information about the objects OBJECT_NOT_AVAILABLE would 
432			// have been returned.)
433			// The GetFromSource() call will eventually lead to a call to
434			// the GetObjects() method of the concrete class derived
435			// from JEventSource. That routine should copy the object
436			// pointers into the factory using the factory's CopyTo()
437			// method which also sets the evnt_called flag for the factory.
438			// Note also that the "t" vector is then filled with a call
439			// to the factory's CopyFrom() method in JEvent::GetObjects().
440			// All we need to do now is just set the factory pointers in
441			// the newly generated JObjects and return the factory pointer.
442 
443			factory->SetFactoryPointers();
444			data_source = DATA_FROM_SOURCE;
445 
446			return factory;
447		}
448	}
449		
450	// OK. It looks like we have to have the factory make this.
451	// Get pointers to data from the factory.
452	factory->Get(t);
453	factory->SetFactoryPointers();
454	data_source = DATA_FROM_FACTORY;
455	
456	return factory;
457}
458 
459//-------------
460// GetFromSource
461//-------------
462template<class T> 
463jerror_t JEventLoop::GetFromSource(vector<const T*> &t, JFactory_base *factory)
464{
465	/// This tries to get objects from the event source.
466	/// "factory" must be a valid pointer to a JFactory
467	/// object since that will take ownership of the objects
468	/// created by the source.
469	/// This should usually be called from JEventLoop::GetFromFactory
470	/// which is called from JEventLoop::Get. The latter will
471	/// create a dummy JFactory of the proper flavor and tag if
472	/// one does not already exist so if objects exist in the
473	/// file without a corresponding factory to create them, they
474	/// can still be used.
475	if(!factory)throw OBJECT_NOT_AVAILABLE;
476	
477	return event.GetObjects(t, factory);
478}
479 
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{
485	/// This is used to fill initial info into a call_stack_t stucture
486	/// for recording the call stack. It should be matched with a call
487	/// to CallStackEnd. It is normally called from the Get() method
488	/// above, but may also be used by external actors to manipulate
489	/// the call stack (presumably for good and not evil).
490 
491	struct itimerval tmr;
492	getitimer(ITIMER_PROFITIMER_PROF, &tmr);
493 
494	cs.caller_name    = this->caller_name;
495	cs.caller_tag     = this->caller_tag;
496	this->caller_name = cs.callee_name = callee_name;
497	this->caller_tag  = cs.callee_tag  = callee_tag;
498	cs.start_time = tmr.it_value.tv_sec + tmr.it_value.tv_usec/1.0E6;
499}
500 
501//-------------
502// CallStackEnd
503//-------------
504inline void JEventLoop::CallStackEnd(JEventLoop::call_stack_t &cs)
505{
506	/// Complete a call stack entry. This should be matched
507	/// with a previous call to CallStackStart which was
508	/// used to fill the cs structure.
509 
510	struct itimerval tmr;
511	getitimer(ITIMER_PROFITIMER_PROF, &tmr);
512	cs.end_time = tmr.it_value.tv_sec + tmr.it_value.tv_usec/1.0E6;
513	caller_name = cs.caller_name;
514	caller_tag  = cs.caller_tag;
515	call_stack.push_back(cs);
516}
517 
518//-------------
519// CheckEventBoundary
520//-------------
521inline bool JEventLoop::CheckEventBoundary(uint64_t event_numberA, uint64_t event_numberB)
522{
523	/// Check whether the two event numbers span one or more boundaries
524	/// in the calibration/conditions database for the current run number.
525	/// Return true if they do and false if they don't. The first parameter
526	/// "event_numberA" is also checked if it lands on a boundary in which
527	/// case true is also returned. If event_numberB lands on a boundary,
528	/// but event_numberA does not, then false is returned.
529	///
530	/// This method is not expected to be called by a user. It is, however called,
531	/// everytime a JFactory's Get() method is called.
532 
533	// Make sure our copy of the boundaries is up to date
534	if(event.GetRunNumber()!=event_boundaries_run){
535		event_boundaries.clear(); // in case we can't get the JCalibration pointer
536		JCalibration *jcalib = GetJCalibration();
537		if(jcalib)jcalib->GetEventBoundaries(event_boundaries);
538		event_boundaries_run = event.GetRunNumber();
539	}
540	
541	// Loop over boundaries
542	for(unsigned int i=0; i<event_boundaries.size(); i++){
543		uint64_t eb = event_boundaries[i];
544		if((eb - event_numberA)*(eb - event_numberB) < 0.0 || eb==event_numberA){ // think about it ....
545			// events span a boundary or is on a boundary. Return true
546			return true;
547		}
548	}
549 
550	return false;
551}
552 
553//-------------
554// FindByID
555//-------------
556template<class T> 
557const T* JEventLoop::FindByID(JObject::oid_t id)
558{
559	/// This is a templated method that can be used in place
560	/// of the non-templated FindByID(oid_t) method if one knows
561	/// the class of the object with the specified id.
562	/// This method is faster than calling the non-templated
563	/// FindByID and dynamic_cast-ing the JObject since
564	/// this will only search the objects of factories that
565	/// produce the desired data type.
566	/// This method will cast the JObject pointer to one
567	/// of the specified type. To use this method,
568	/// a type is specified in the call as follows:
569	///
570	/// const DMyType *t = loop->FindByID<DMyType>(id);
571	
572	// Loop over factories looking for ones that provide
573	// specified data type.
574	for(unsigned int i=0; i<factories.size(); i++){
575		if(factories[i]->GetDataClassName() != T::static_className())continue;
576 
577		// This factory provides data of type T. Search it for
578		// the object with the specified id.
579		const JObject *my_obj = factories[i]->GetByID(id);
580		if(my_obj)return dynamic_cast<const T*>(my_obj);
581	}
582 
583	return NULL__null;
584}
585 
586//-------------
587// GetCalib (map)
588//-------------
589template<class T>
590bool JEventLoop::GetCalib(string namepath, map<string,T> &vals)
591{
592	/// Get the JCalibration object from JApplication for the run number of
593	/// the current event and call its Get() method to get the constants.
594 
595	// Note that we could do this by making "vals" a generic type T thus, combining
596	// this with the vector version below. However, doing this explicitly will make
597	// it easier for the user to understand how to call us.
598 
599	vals.clear();
600 
601	JCalibration *calib = GetJCalibration();
602	if(!calib){
603		_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;
604		return true;
605	}
606	
607	return calib->Get(namepath, vals, event.GetEventNumber());
608}
609 
610//-------------
611// GetCalib (vector)
612//-------------
613template<class T> bool JEventLoop::GetCalib(string namepath, vector<T> &vals)
614{
615	/// Get the JCalibration object from JApplication for the run number of
616	/// the current event and call its Get() method to get the constants.
617 
618	vals.clear();
619 
620	JCalibration *calib = GetJCalibration();
621	if(!calib){
622		_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;
623		return true;
624	}
625	
626	return calib->Get(namepath, vals, event.GetEventNumber());
627}
628 
629//-------------
630// GetCalib (single)
631//-------------
632template<class T> bool JEventLoop::GetCalib(string namepath, T &val)
633{
634	/// This is a convenience method for getting a single entry. It
635	/// simply calls the vector version and returns the first entry.
636	/// It returns true if the vector version returns true AND there
637	/// is at least one entry in the vector. No check is made for there
638	/// there being more than one entry in the vector.
639 
640	vector<T> vals;
641	bool ret = GetCalib(namepath, vals);
642	if(vals.empty()) return true;
643	val = vals[0];
644	
645	return ret;
646}
647 
648//-------------
649// GetGeom (map)
650//-------------
651template<class T>
652bool JEventLoop::GetGeom(string namepath, map<string,T> &vals)
653{
654	/// Get the JGeometry object from JApplication for the run number of
655	/// the current event and call its Get() method to get the constants.
656 
657	// Note that we could do this by making "vals" a generic type T thus, combining
658	// this with the vector version below. However, doing this explicitly will make
659	// it easier for the user to understand how to call us.
660 
661	vals.clear();
662 
663	JGeometry *geom = GetJGeometry();
664	if(!geom){
665		_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;
666		return true;
667	}
668	
669	return geom->Get(namepath, vals);
670}
671 
672//-------------
673// GetGeom (atomic)
674//-------------
675template<class T> bool JEventLoop::GetGeom(string namepath, T &val)
676{
677	/// Get the JCalibration object from JApplication for the run number of
678	/// the current event and call its Get() method to get the constants.
679 
680	JGeometry *geom = GetJGeometry();
681	if(!geom){
682		_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;
683		return true;
684	}
685	
686	return geom->Get(namepath, val);
687}
688 
689//-------------
690// SetRef
691//-------------
692template<class T>
693void JEventLoop::SetRef(T *t)
694{
695	pair<const char*, void*> p(typeid(T).name(), (void*)t);
696	user_refs.push_back(p);
697}
698 
699//-------------
700// GetResource
701//-------------
702template<class T> bool JEventLoop::GetResource(string namepath, T vals, int event_number)
703{
704	JResourceManager *resource_manager = GetJResourceManager();
705	if(!resource_manager){
706		string mess = string("Unable to get the JResourceManager object (namepath=\"")+namepath+"\")";
707		throw JException(mess);
708	}
709 
710	return resource_manager->Get(namepath, vals, event_number);
711}
712 
713//-------------
714// GetRef
715//-------------
716template<class T>
717T* JEventLoop::GetRef(void)
718{
719	/// Get a user-defined reference (a pointer)
720	for(unsigned int i=0; i<user_refs.size(); i++){
721		if(user_refs[i].first == typeid(T).name()) return (T*)user_refs[i].second;
722	}
723 
724	return NULL__null;
725}
726 
727//-------------
728// GetRefsT
729//-------------
730template<class T>
731vector<T*> JEventLoop::GetRefsT(void)
732{
733	vector<T*> refs;
734	for(unsigned int i=0; i<user_refs.size(); i++){
735		if(user_refs[i].first == typeid(T).name()){
736			refs.push_back((T*)user_refs[i].second);
737		}
738	}
739 
740	return refs;
741}
742 
743//-------------
744// RemoveRef
745//-------------
746template<class T>
747void JEventLoop::RemoveRef(T *t)
748{
749	vector<pair<const char*, void*> >::iterator iter;
750	for(iter=user_refs.begin(); iter!= user_refs.end(); iter++){
751		if(iter->second == (void*)t){
752			user_refs.erase(iter);
753			return;
754		}
755	}
756	_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}
758 
759 
760#endif //__CINT__  __CLING__
761 
762} // Close JANA namespace
763 
764 
765 
766#endif // _JEventLoop_
767 

←

/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(); }
19
←
Assuming the condition is true→
20
←
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