Updated changes from main

.github/pull_request_template.md

@@ -0,0 +1,12 @@
+## Describe your changes
+
+## Checklist before submitting your PR
+- [ ] This PR implements a single change (one new/modified Tool, or a set of changes to implement one new/modified feature)
+- [ ] This PR alters the minimum number of files to affect this change
+- [ ] If this PR includes a new Tool, a README and minimal demonstration ToolChain is provided
+- [ ] If a new Tool/ToolChain requires model or configuration files, their paths are not hard-coded, and means of generating those files is described in the readme, with examples provided on /pnfs/annie/persistent
+- [ ] For every `new` usage, there is a reason the data must be on the heap
+- [ ] For every `new` there is a `delete`, unless I explicitly know why (e.g. ROOT or a BoostStore takes ownership)
+
+## Additional Material
+Attach any validation or demonstration files here. You may also link to relavant docdb articles.

DataModel/ANNIEGeometry.cpp

@@ -190,8 +190,8 @@ void ANNIEGeometry::SetGeometry(){
   //WCSim
   fCylRadius = 152.0; //cm
   fCylLength = 396.0; //cm
-  fCylFiducialRadius = 0.0;
-  fCylFiducialLength = 0.0;
+  fCylFiducialRadius = 100.0; //cm
+  fCylFiducialLength = 200.0; //cm
   fXoffset = 0.0;
   fYoffset = 0.0; 
   fZoffset = 0.0;
@@ -337,11 +337,11 @@ bool ANNIEGeometry::InsideDetector(double x, double y, double z)
 bool ANNIEGeometry::InsideFiducialVolume(double x, double y, double z)
 {
   if( fGeoType==ANNIEGeometry::kCylinder ){
-    if( z>=-0.5*fCylFiducialLength && z<=+0.5*fCylFiducialLength
-     && x*x+y*y<=fCylFiducialRadius*fCylFiducialRadius ){
+    if( y>=-0.5*fCylFiducialLength && y<=+0.5*fCylFiducialLength
+     && x*x+z*z<=fCylFiducialRadius*fCylFiducialRadius ){
       return 1;
     }
-  } else std::cout<<"WTF ANNIE Should Be A Cylinder!!!!"<<std::endl;
+  } else std::cout<<"WTH ANNIE Should Be A Cylinder!!!!"<<std::endl;
 
   return 0;
 }
@@ -377,42 +377,43 @@ double ANNIEGeometry::DistanceToEdge(double x, double y, double z)
       double dr = 0.0;
       if( fCylRadius>dr ) dr = fCylRadius;
       if( 0.5*fCylLength>dr ) dr = 0.5*fCylLength;
-      if( -sqrt(x*x+y*y)+fCylRadius<dr ) dr = -sqrt(x*x+y*y)+fCylRadius;
-      if( -z+0.5*fCylLength<dr ) dr = -z+0.5*fCylLength;
-      if( +z+0.5*fCylLength<dr ) dr = +z+0.5*fCylLength;
+      if( -sqrt(x*x+z*z)+fCylRadius<dr ) dr = -sqrt(x*x+z*z)+fCylRadius;
+      if( -y+0.5*fCylLength<dr ) dr = -y+0.5*fCylLength;
+      if( y+0.5*fCylLength<dr ) dr = y+0.5*fCylLength;
       return dr;
     }
 
     // outside detector (convention: -ve dr)
     else{
 
       // side region
-      if( z>=-0.5*fCylLength && z<=+0.5*fCylLength ){
-        return -sqrt(x*x+y*y)+fCylRadius;
+      if( y>=-0.5*fCylLength && y<=+0.5*fCylLength ){
+        return -sqrt(x*x+z*z)+fCylRadius;
       }
 
-      // top region
-      if( z<=-0.5*fCylLength
-       && x*x+y*y<fCylRadius*fCylRadius ){
-        return +z+0.5*fCylLength;
+      // below tank
+      if( y<=-0.5*fCylLength
+       && x*x+z*z<fCylRadius*fCylRadius ){
+        return y+0.5*fCylLength;
       }
-      if( z>=+0.5*fCylLength
-       && x*x+y*y<fCylRadius*fCylRadius ){
-        return -z+0.5*fCylLength;
+      // above tank
+      if( y>=+0.5*fCylLength
+       && x*x+z*z<fCylRadius*fCylRadius ){
+        return -y+0.5*fCylLength;
       }
 
       // corner regions
-      if( z>=+0.5*fCylLength
-       && x*x+y*y>=fCylRadius ){
-        double dr = sqrt(x*x+y*y)-fCylRadius;
-        double dz = -z+0.5*fCylLength;
-        return -sqrt(dr*dr+dz*dz);
+      if( y>=+0.5*fCylLength
+       && x*x+z*z>=fCylRadius ){
+        double dr = sqrt(x*x+z*z)-fCylRadius;
+        double dy = -y+0.5*fCylLength;
+        return -sqrt(dr*dr+dy*dy);
       }
-      if( z<=-0.5*fCylLength
-       && x*x+y*y>=fCylRadius ){
-        double dr = sqrt(x*x+y*y)-fCylRadius;
-        double dz = +z+0.5*fCylLength;
-        return -sqrt(dr*dr+dz*dz);
+      if( y<=-0.5*fCylLength
+       && x*x+z*z>=fCylRadius ){
+        double dr = sqrt(x*x+z*z)-fCylRadius;
+        double dy = y+0.5*fCylLength;
+        return -sqrt(dr*dr+dy*dy);
       }
     }
   }

GetToolDAQ.sh

@@ -497,6 +497,7 @@ then
     pip3 install uproot==4.3.7
     pip3 install xgboost==1.6.2
     pip3 install tensorflow==2.10.0
+    pip3 install torch==1.12.1
     pip3 install PyQt5
     # set tensorflow verbosity to suppress info messages about not having a GPU or maximal acceleration
     # https://stackoverflow.com/questions/35911252/disable-tensorflow-debugging-information/42121886#42121886

Makefile

@@ -94,7 +94,7 @@ clean:
 	rm -f include/*.h
 	rm -f lib/*.so
 	rm -f Analyse
-	rm -f UserTools/*/*.o
+	find UserTools/* -type f -name '*.o' -follow -writable -delete
 	rm -f DataModel/*.o
 	rm -f DataModel/DataModel_Linkdef.hh
 	rm -f DataModel/DataModel_RootDict*

Setup.sh

@@ -4,6 +4,8 @@
 
 ToolDAQapp=`pwd`
 
+ulimit -n 4096
+
 export LIBGL_ALWAYS_INDIRECT=1
 
 source ${ToolDAQapp}/ToolDAQ/root-6.24.06/install/bin/thisroot.sh
@@ -24,7 +26,7 @@ export PATH=${ToolDAQapp}/ToolDAQ/LHAPDF-6.3.0/install/bin:$GENIE/bin:$GENIE_REW
 export PATH=/ToolAnalysis/ToolDAQ/fsplit:$PATH
 export TF_CPP_MIN_LOG_LEVEL=2
 
-export FW_SEARCH_PATH=${ToolDAQapp}/UserTools/ReweightEventsGenie
+export FW_SEARCH_PATH=${ToolDAQapp}/UserTools/ReweightFlux
 
 for folder in `ls -d ${PWD}/UserTools/*/ `
 do

UserTools/AssignBunchTimingMC/AssignBunchTimingMC.cpp

@@ -0,0 +1,245 @@
+#include "AssignBunchTimingMC.h"
+
+AssignBunchTimingMC::AssignBunchTimingMC():Tool(){}
+
+//------------------------------------------------------------------------------
+
+bool AssignBunchTimingMC::Initialise(std::string configfile, DataModel &data){
+
+    // Get configuration variables and set default values if necessary
+
+    if ( !configfile.empty() ) m_variables.Initialise(configfile);
+    m_data = &data;
+
+    bool got_verbosity = m_variables.Get("verbosity", verbosity);
+    bool got_width     = m_variables.Get("bunchwidth", fbunchwidth);
+    bool got_interval  = m_variables.Get("bunchinterval", fbunchinterval);
+    bool got_count     = m_variables.Get("bunchcount", fbunchcount);
+    bool got_sample    = m_variables.Get("sampletype", fsample);
+    bool got_trigger   = m_variables.Get("prompttriggertime", ftriggertime);
+
+    if (!got_verbosity) {
+        verbosity = 0;
+        logmessage = "Warning (AssignBunchTimingMC): \"verbosity\" not set in the config, defaulting to 0";
+        Log(logmessage, v_warning, verbosity);
+    }
+
+
+    // default bunch parameters from MicroBooNE paper: https://doi.org/10.1103/PhysRevD.108.052010
+
+    if (!got_width) {
+        fbunchwidth = 1.308;  // ns
+        logmessage = ("Warning (AssignBunchTimingMC): \"bunchwidth\" not "
+            "set in the config file. Using default: 1.308ns");
+        Log(logmessage, v_warning, verbosity);
+    }
+
+    if (!got_interval) {
+        fbunchinterval = 18.936;  // ns
+        logmessage = ("Warning (AssignBunchTimingMC): \"bunchinterval\" not "
+            "set in the config file. Using default: 18.936ns");
+        Log(logmessage, v_warning, verbosity);
+    }
+
+    if (!got_count) {
+        fbunchwidth = 81;
+        logmessage = ("Warning (AssignBunchTimingMC): \"bunchcount\" not "
+            "set in the config file. Using default: 81");
+        Log(logmessage, v_warning, verbosity);
+    }
+
+    if (!got_sample) {
+        fsample = 0;    // assume they are using the Tank samples
+        logmessage = ("Warning (AssignBunchTimingMC): \"sampletype\" not "
+            "set in the config file. Using default: 0 (GENIE tank samples)");
+        Log(logmessage, v_warning, verbosity);
+    }
+
+    if (!got_trigger) {
+        ftriggertime = 0;    // assume they are using the old WCSim prompt trigger time (t = 0, first particle arrival)
+        logmessage = ("Warning (AssignBunchTimingMC): \"prompttriggertime\" not "
+            "set in the config file. Using default: 0 (WCSim prompt trigger time = 0)");
+        Log(logmessage, v_warning, verbosity);
+    }
+
+
+    if (verbosity >= v_message) {
+        std::cout<<"------------------------------------"<<"\n";
+        std::cout<<"AssignBunchTimingMC: Config settings"<<"\n";
+        std::cout<<"------------------------------------"<<"\n";
+        std::cout<<"bunch width       = "<<fbunchwidth<<" ns"<<"\n";
+        std::cout<<"bunch interval    = "<<fbunchinterval<<" ns"<<"\n";
+        std::cout<<"number of bunches = "<<fbunchcount<<"\n";
+        std::cout<<"sample type       = "<<(fsample == 0 ? "(0) Tank" : "(1) World")<<"\n";
+        std::cout<<"trigger time      = "<<(ftriggertime == 0 ? "(0) prompt trigger starts when first particle arrives (default WCSim)" 
+                                                   : "(1) prompt trigger starts at beam dump (modified WCSim)")<<"\n";
+    }
+
+    fbunchTimes = new std::map<double, double>;   // set up pointer; Store will handle deletion in Finalize
+
+    return true;
+}
+
+//------------------------------------------------------------------------------
+
+
+bool AssignBunchTimingMC::Execute()
+{
+
+    if (verbosity >= v_debug) {
+        std::cout << "AssignBunchTimingMC: Executing tool..." << std::endl;
+    }
+
+    if (!LoadStores())      // Load info from store
+        return false;
+    if (verbosity >= v_debug) {
+        std::cout << "AssignBunchTimingMC: Store info loading successful" << std::endl;
+    }
+
+    fbunchTimes->clear();   // clear map
+
+    BNBtiming();            // grab BNB structure
+    if (verbosity >= v_debug) {
+        std::cout << "AssignBunchTimingMC: BNB timing successful" << std::endl;
+    }
+
+    for (std::pair<double, std::vector<MCHit>>&& apair : *fClusterMapMC) {
+        double totalHitTime = 0;
+        int hitCount = 0;
+        int totalHits = apair.second.size();
+
+        CalculateClusterAndBunchTimes(apair.second, totalHitTime, hitCount, totalHits);
+
+        // store the cluster time in a map (e.g., keyed by cluster identifier)
+        fbunchTimes->emplace(apair.first, bunchTime);
+    }
+
+    if (verbosity >= v_debug) {
+        std::cout << "AssignBunchTimingMC: Assigned bunch time successfully. Writing to Store..." << std::endl;
+    }
+
+    // send to store
+    m_data->Stores.at("ANNIEEvent")->Set("bunchTimes", fbunchTimes);
+
+    return true;
+
+}
+
+
+//------------------------------------------------------------------------------
+
+bool AssignBunchTimingMC::Finalise()
+{
+  return true;
+}
+
+//------------------------------------------------------------------------------
+
+
+
+
+bool AssignBunchTimingMC::LoadStores()
+{
+    // grab necessary information from Stores
+
+    bool get_MCClusters = m_data->CStore.Get("ClusterMapMC", fClusterMapMC);
+    if (!get_MCClusters) {
+        Log("AssignBunchTimingMC: no ClusterMapMC in the CStore! Are you sure you ran the ClusterFinder tool?", v_error, verbosity);
+        return false;
+    }
+
+    bool get_AnnieEvent = m_data->Stores.count("ANNIEEvent");
+    if (!get_AnnieEvent) {
+        Log("AssignBunchTimingMC: no ANNIEEvent store!", v_error, verbosity);
+        return false;
+    }
+
+    bool get_MCHits = m_data->Stores.at("ANNIEEvent")->Get("MCHits", fMCHitsMap);
+    if (!get_MCHits) {
+        Log("AssignBunchTimingMC: no MCHits in the ANNIEEvent!", v_error, verbosity);
+        return false;
+    }
+
+    bool get_neutrino_vtxt = m_data->Stores["GenieInfo"]->Get("NuIntxVtx_T",TrueNuIntxVtx_T);   // grab neutrino interaction time
+    if (!get_neutrino_vtxt) {
+        Log("AssignBunchTimingMC: no GENIE neutrino interaction time info in the ANNIEEvent! Are you sure you ran the LoadGENIEEvent tool?", v_error, verbosity);
+        return false;
+    }
+
+    return true;
+}
+
+
+//------------------------------------------------------------------------------
+
+
+void AssignBunchTimingMC::BNBtiming()
+{
+    // Determined from GENIE samples (as of Dec 2024)
+    const double tank_time = 33.0;    // Tank neutrino arrival time: 33ns
+    const double world_time = 33.0;   // WORLD neutrino arrival time: 33ns  (As of Dec 2024, World samples have not yet been re-produced fully)
+
+    if (ftriggertime == 0) {
+        new_nu_time = (fsample == 0) ? (TrueNuIntxVtx_T - tank_time) : (TrueNuIntxVtx_T - world_time);
+    } else if (ftriggertime == 1) {
+        new_nu_time = (fsample == 0) ? -tank_time : -world_time;
+    }
+
+    // seed random number generator with the current time
+    unsigned seed = std::chrono::system_clock::now().time_since_epoch().count();
+    generator.seed(seed);
+    if (verbosity >= v_debug) {
+        std::cout << "AssignBunchTimingMC: Random seed selected: " << seed << std::endl;
+    }
+
+
+    // per event, assign BNB structure (random bunch number + jitter based on intrinsic bunch width + new_nu_time)
+    bunchNumber = rand() % fbunchcount;
+    std::normal_distribution<double> distribution(0, fbunchwidth); 
+    jitter = distribution(generator);
+
+    if (verbosity >= v_message) {
+        std::string logmessage = "AssignBunchTimingMC: bunchNumber = " + std::to_string(bunchNumber) + " | t0 = " + std::to_string(new_nu_time);
+        Log(logmessage, v_debug, verbosity);
+    }
+
+}
+
+
+//------------------------------------------------------------------------------
+
+
+void AssignBunchTimingMC::CalculateClusterAndBunchTimes(std::vector<MCHit> const &mchits, double &totalHitTime, int &hitCount, int &totalHits)
+{
+
+    // loop over the hits to get their times
+    for (auto mchit : mchits) { 
+        double hitTime = mchit.GetTime();
+        totalHitTime += hitTime;
+        hitCount++;
+        if (verbosity >= v_debug) {
+            std::string logmessage = "AssignBunchTimingMC: (" + std::to_string(hitCount) + "/" + std::to_string(totalHits) + ") MCHit time = " + std::to_string(hitTime);
+            Log(logmessage, v_debug, verbosity);
+        }
+    }
+
+    // find nominal cluster time (average hit time)
+    double clusterTime = (hitCount > 0) ? totalHitTime / hitCount : -9999;
+    if (verbosity >= v_debug) {
+        std::string logmessage = "AssignBunchTimingMC: ClusterTime = " + std::to_string(clusterTime);
+        Log(logmessage, v_debug, verbosity);
+    }
+
+    // calculate BunchTime
+    bunchTime = fbunchinterval * bunchNumber + clusterTime + jitter + new_nu_time ;
+    if (verbosity >= v_debug) {
+        std::string logmessage = "AssignBunchTimingMC: bunchTime = " + std::to_string(bunchTime);
+        Log(logmessage, v_debug, verbosity);
+    }
+
+}
+
+
+//------------------------------------------------------------------------------
+
+// done