Muon matching code

#Setup.sh#

@@ -0,0 +1,19 @@
+#!/bin/bash
+
+#Application path location of applicaiton
+
+Dependencies=`pwd`/Dependencies
+
+export LD_LIBRARY_PATH=`pwd`/lib:${Dependencies}/TMVA/lib:$LD_LIBRARY_PATH
+
+for folder in `ls -d ${PWD}/UserTools/*/ `
+do
+    export PYTHONPATH=$folder:${PYTHONPATH}
+done
+
+export SEGFAULT_SIGNALS="all"
+export LD_LIBRARY_PATH=$PWD/Dependencies/relic_sk4_ana/data_reduc/third/lib:$LD_LIBRARY_PATH
+source /usr/local/sklib_gcc8/bashenv_gcc8_skofl_21a+atmpd_21a
+source /usr/local/sklib_gcc8/skofl_21a/env.sh
+
+export NEUTRON=1

.gitignore

@@ -16,3 +16,8 @@ DataModel/NTagDataModelDict.cxx
 BStore_RootDict.cxx
 BStore_RootDict.h
 BStore_Linkdef_1.hh
+.cache
+.dir-locals.el
+*.json
+*.png
+*.pdf

DataModel/Constants.cpp

@@ -135,3 +135,11 @@ double PdgToMass(int code){
 	return particle->Mass()*1000.;      // converted to MeV
 }
 
+std::string TriggerIDToTrigger(int code){
+	if(constants::Trigger_ID_To_Trigger.count(code)){
+		return constants::Trigger_ID_To_Trigger.at(code);
+	} else {
+		return std::to_string(code);
+	}
+}
+

DataModel/Constants.h

@@ -24,6 +24,7 @@ int G3ParticleCodeToPdg(int code);
 int StringToG3ParticleCode(std::string name);
 std::string G3ParticleCodeToString(int code);
 double PdgToMass(int code);
+std::string TriggerIDToTrigger(int code);
 
 enum class SKROOTMODE : int { NONE = 4, ZEBRA = 3, READ = 2, WRITE = 1, COPY = 0 };
 
@@ -498,6 +499,43 @@ namespace constants{
 		{"also_multiple_mu",muboy_classes::also_multiple_mu}
 	};
 
+	static const std::map<int, std::string> Trigger_ID_To_Trigger{
+		// from skhead.h
+		// "BLANK #" entries are empty ID numbers in skhead.h
+		{0, "LE (sftw. trig.)"},
+		{1, "HE (sftw. trig.)"},
+		{2, "SLE (sftw. trig.)"},
+		{3, "OD (sftw. trig.)"},
+		{4, "Periodic (SKI-III)"},
+		{5, "After/CAL (SKI-III)"},
+		{6, "VETO START"},
+		{7, "VETO STOP"},
+		{8, "BLANK 8"},
+		{9, "BLANK 9"},
+		{10, "BLANK 10"},
+		{11, "Random Wide Trigger"},
+		{12, "Laser (ID, Usho Laser)"},
+		{13, "LED"},
+		{14, "Ni"},
+		{15, "Laser (OD, AutoTQlaser)"},
+		{16, "LE (hitsum)"},
+		{17, "HE (histum)"},
+		{18, "SLE (hitsum)"},
+		{19, "OD (hitsum)"},
+		{20, "BLANK 20"},
+		{21, "BLANK 21"},
+		{22, "SN Burst"},
+		{23, "mu->e Decay"},
+		{24, "LINAC"},
+		{25, "LINAC Microwave"},
+		{26, "BLANK 26"},
+		{27, "Periodic (simple)"},
+		{28, "SHE"},
+		{29, "AFT"},
+		{30, "Pedestal"},
+		{31, "T2K"}
+	};
+
 	// maps of regular expression error codes to descriptive strings
 	// -------------------------------------------------------------
 	/*

DataModel/DataModel.cpp

@@ -14,7 +14,7 @@ DataModel::DataModel() : eventVariables_p(new BStore(false,true)), eventVariable
 }
 
 DataModel::~DataModel(){
-	if(rootTApp) delete rootTApp;
+	//if(rootTApp) delete rootTApp; // segfaults on application termination. J. Fannon 
 	//if(connectionTable) delete connectionTable;   // segfaults on application termination, maybe?
 }
 
@@ -57,15 +57,29 @@ TApplication* DataModel::GetTApp(){
 	return rootTApp;
 }
 
-bool DataModel::RegisterReader(std::string readerName, MTreeReader* reader, std::function<bool()> hasAFT, std::function<bool()> loadSHE, std::function<bool()> loadAFT, std::function<bool(int)> loadCommon){
+bool DataModel::RegisterReader(std::string readerName, MTreeReader* reader, std::function<bool()> hasAFT, std::function<bool()> loadSHE, std::function<bool()> loadAFT, std::function<bool(int)> loadCommon, std::function<int(long, bool)> getTreeEntry){
 	Trees.emplace(readerName, reader);
 	hasAFTs.emplace(readerName, hasAFT);
 	loadSHEs.emplace(readerName, loadSHE);
 	loadAFTs.emplace(readerName, loadAFT);
 	loadCommons.emplace(readerName, loadCommon);
+	getEntries.emplace(readerName, getTreeEntry);
 	return true;
 }
 
+int DataModel::getTreeEntry(std::string ReaderName, long entrynum){
+	if(ReaderName==""){
+		// if no name given but we have only one TreeReader Tool, use that
+		if(hasAFTs.size()) ReaderName = hasAFTs.begin()->first;
+	}
+	if(getEntries.count(ReaderName)){
+		return getEntries.at(ReaderName)(entrynum, true);
+	} else {
+		std::cerr << "getTreeEntry requested for Unknown reader "<<ReaderName<<std::endl;
+	}
+	return 0;
+}
+
 bool DataModel::HasAFT(std::string ReaderName){
 	if(ReaderName==""){
 		// if no name given but we have only one TreeReader Tool, use that

DataModel/DataModel.h

@@ -30,6 +30,10 @@
 #include "EventTrueCaptures.h"
 #include "PMTHitCluster.h"
 
+#include "ParticleCand.h"
+
+#include "skroot_loweC.h"
+
 class MTreeReader;
 class MTreeSelection;
 class TreeReader;
@@ -67,14 +71,16 @@ class DataModel {
   std::unordered_map<std::string, std::function<bool()>> loadSHEs;
   std::unordered_map<std::string, std::function<bool()>> loadAFTs;
   std::unordered_map<std::string, std::function<bool(int)>> loadCommons;
+  std::unordered_map<std::string, std::function<int(long, bool)>> getEntries;
 
   Store vars; ///< This Store can be used for any variables. It is an inefficent ascii based storage and command line arguments will be placed in here along with ToolChain variables
   BStore CStore; ///< This is a more efficent binary Store that can be used to store a dynamic set of inter Tool variables, very useful for constants and and flags hence the name CStore
   std::map<std::string,BStore*> Stores;  ///< This is a map of named BStore pointers which can be deffined to hold a nammed collection of any type of BStore. It is usefull to store data collections that needs subdividing into differnt stores.
 
   // This function is used to register a TreeReader tool's member functions with the DataModel,
   // which provides access from other Tools
-  bool RegisterReader(std::string readerName, MTreeReader* reader, std::function<bool()> hasAFT, std::function<bool()> loadSHE, std::function<bool()> loadAFT, std::function<bool(int)> loadCommon);
+  bool RegisterReader(std::string readerName, MTreeReader* reader, std::function<bool()> hasAFT, std::function<bool()> loadSHE, std::function<bool()> loadAFT, std::function<bool(int)> loadCommon, std::function<int(long,bool)> getTreeEntry);
+  int getTreeEntry(std::string ReaderName="", long entrynum=0);
   // These retain function pointers to call the corresponding TreeReader functions.
   // The TreeReader instance is obtained from the name specified in their config file.
   bool HasAFT(std::string ReaderName="");
@@ -88,8 +94,8 @@ class DataModel {
   // wrapper to check if we've called this yet, since we should probably only call it the once?
   void KZInit();
   bool kz_initialized=false;
-
-
+  
+  
   // Event vars
   BStore* eventVariables_p; // TODO replace with a pointer and update tools to use -> instead of .
   BStore &eventVariables;   // use references to preserve current behaviour...
@@ -107,7 +113,34 @@ class DataModel {
 
   std::map<std::string, Store*> tool_configs;
   StoreToTTree StoreConverter;
-
+
+  //flag for a new muon - used for the spallation reduction
+  bool newMuon = false;
+
+  //flag for a new relic candidate - used for the spallation reduction
+  bool newRelic = false;
+
+  //deques of the struct ParticleCand to store event info of ALL muon candidates and relic candidates
+  std::deque<ParticleCand> muonCandDeque;
+  std::deque<ParticleCand> relicCandDeque;
+
+  //deque for muons that need to be reconstructed (if they have matched to a relic they need to be reconstructed)
+  std::vector<ParticleCand> muonsToRec;
+
+  //vector for relic candidate events that needs to be written out - used for the spallation reduction
+  std::vector<ParticleCand> writeOutRelics;
+
+  //map for lowe information to store when reconstructing events mid-chain and wanting to retrieve the information later
+  std::map<long, skroot_lowe_common> loweCommonBufferMap;
+
+  //vector to store weightings in from various SRN models at energy intervals of 0.5 MeV.
+  std::vector<std::vector<float>> SRNWeights;
+
+  //vector to store modelNames for the SRN Models in
+  std::vector<std::string> SRNModelNames;
+
+  bool applyReweight = false;
+
  private:
 
 

DataModel/MTreeReader.cpp

@@ -187,7 +187,7 @@ int MTreeReader::ParseBranches(){
 			TClass* ac = TClass::GetClass(lf->GetTypeName());
 			if(ac!=nullptr){
 				branch_istobject.emplace(branchname,ac->InheritsFrom("TObject"));
-			} else {
+			}else {
 				std::cerr<<"Unknown class for branch "<<lf->GetBranch()->GetTitle()
 						 <<"! Please make a dictionary."<<std::endl;
 				// what are the consequences of this? Do we need to remove it from

DataModel/ParticleCand.h

@@ -0,0 +1,17 @@
+#ifndef PARTICLECAND_H
+#define PARTICLECAND_H
+
+#include "skroot_loweC.h"
+
+struct ParticleCand {
+	int EventNumber;
+	float EventTime;
+	int EntryNumber;
+	int PID = 0; //0 = muon 1 = LowE
+	std::vector<int> matchedParticleEvNum;
+	std::vector<float> matchedParticleTimeDiff;
+	std::vector<float> matchedParticleBSEnergy;
+	skroot_lowe_common LowECommon;
+};
+
+#endif

DataModel/fortran_routines.h

@@ -46,13 +46,33 @@ extern "C" void skcread_(int*, int*);
 extern "C" void skroot_set_tree_(int*);
 extern "C" void skroot_get_entry_(int*);
 
+extern "C" int softtrg_inittrgtbl_(int*, int*, int*,int*);
+extern "C" int softtrg_inittrgtbl(int , int, int, int);
+
+extern "C" void fix_maxqisk_();
+extern "C" void lfmufit_sk4_();
+
+extern "C" void newmufit_(float (*)[3], float (*)[3], float*);
+
+//extern "C" void makededx_(float (*)[4], float (*)[3], int (*)[11146], float (*)[11146], float (*)[11146], float (*)[11146][3], int*, float (*)[200]);
+
+extern "C" void makededx_intg_(float (*)[4], float (*)[3], float*, int (*)[11146], float (*)[11146], float (*)[11146], float (*)[11146][3], int*, int*, float (*)[200], int (*)[334380], int*);
+
+extern "C" void mfmuselect_(float (*)[3], float (*)[3], float*, int*);
+
+extern "C" void mffastfast_(float (*)[3], float (*)[3], int*);
+
+extern "C" void muboy_zbs_(int*, int*, float (*)[4], float (*)[3], float*, float*, int*, float (*)[36], int*);
+
 // from $ATMPD_ROOT/src/programs/TreeBuilder/examples/fort_fopen.F
 extern "C" void fort_fopen_(int*, const char*, char*, int* ,int);
 
 // read runinf (y/m/d, start time, end time, etc)
 extern "C" void runinfsk_();
 
 //
+
+extern "C" void softtrg_get_cond_(int*, int*, int*, int*, int*);
 extern "C" void softtrg_set_cond_(int*, int*, int*, int*, int*);
 extern "C" void get_sub_triggers_(int*, int*, int*, int*);
 extern "C" void set_timing_gate_(int*);
@@ -75,6 +95,7 @@ extern "C" void lfallfit_sk4_final_qe43_(float*, int*, int*, int*, int*);
 extern "C" void lfallfit_sk4_data_(float*, int*, int*);
 extern "C" void lfallfit_sk4_gain_corr_(float*, int*, int*, int*, int*);
 extern "C" void lfallfit_sk4_mc_(float*, int*, int*);
+extern "C" void lfallfit_sk6_mc_(float*, int*, int*, int*, int*);
 // skroot_lowe_ common block
 
 extern "C" void rluxgo_(int*, int*, int*, int*);

FluxModels/ando.txt

@@ -0,0 +1,33 @@
+3 1.05239
+4 1.26815
+5 1.28199
+6 1.18259
+7 1.0366
+8 0.88051
+9 0.732413
+10 0.601257
+11 0.489383
+12 0.396725
+13 0.320967
+14 0.259602
+15 0.210367
+16 0.17068
+17 0.138815
+18 0.113143
+19 0.0923727
+20 0.075644
+21 0.0620584
+22 0.0510444
+23 0.0420982
+24 0.0347902
+25 0.0288424
+26 0.023958
+27 0.0199445
+28 0.0166415
+29 0.0139067
+30 0.0116483
+31 0.00977277
+32 0.008211
+33 0.00691467
+34 0.00582942
+35 0.00492365