Getting started in SBND

Viewing charged particles in SBND using a particle gun generator and analyzer module

(Requires some knowledge of reconstruction chains, c++, python and a recent version of sbndcode. I ran on v09\_43_00)

If you still need to setup your sbndcode area, follow steps 1-3 on the SBND Commissioning Page - Get Started Guide

1 - Visual Studio Code

Although there are many text editor options, I’d recommend Visual Studio Code for many reasons:

  1. It has syntax highlighting for almost every language including Fermilab Hierarchical Configuration Language or .fcl files.
  2. You can view .root files using the Explorer page.
  3. It supports ssh connections with X11 forwarding.
  4. There are many keyboard shortcuts such as multiline commenting and search and replace.

To get started, let’s first add an ssh host on VSCode. Whichever gpvm machine you use should be fine. Next, we should add extensions. For this tutorial, you’ll need the following extensions:

-Fermilab Hierarchical Configuration Language support
-Python and if you use Jupyter Notebooks
-Root File Viewer

Now we can get started on doing physics in SBND! First, open the terminal in VSCode to make a directory and your .fcl file:

mkdir /sbnd/app/users/$USER/tutorial
mkdir /sbnd/app/users/$USER/tutorial/data 
cd /sbnd/app/users/$USER/tutorial/data
touch muon_gun.fcl

We will edit this .fcl file in the next section.

2 - Particle Gun

Using a particle gun is a simple way to simulate a single particle, or multiple particles at a time in the SBND detector (This can be done in ICARUS as well but I’ve done it here in SBND). Let’s look at a particle gun .fcl file and see what’s going on! Add the following lines to the file muon_gun.fcl

#Include .fcl for producing particles
#include "prodsingle_sbnd_proj.fcl" 

physics.producers.generator.PadOutVectors: true  #Duplicates single element vectors to match length of longest vector
physics.producers.generator.PDG: [13]            #Generate 1 muon
physics.producers.generator.P0: [3]              #p = 3 GeV/c 
physics.producers.generator.SigmaP: [0]          #No variance
physics.producers.generator.SigmaX: [0]          #
physics.producers.generator.SigmaY: [0]          #    
physics.producers.generator.SigmaZ: [0]          #    
physics.producers.generator.PosDist: 0           #0:Uniform 1:Gaussian         
physics.producers.generator.X0: [-150]           #Start track in TPC0 (x<0)
physics.producers.generator.Y0: [-100]           #Start in lower half of detector
physics.producers.generator.Z0: [-50]            #Start 50 cm upstream
physics.producers.generator.Theta0XZ: [25]       #Muon trajectory in X-Z plane
physics.producers.generator.Theta0YZ: [30]       #Muon trajectory in Y-Z plane
physics.producers.generator.SigmaThetaXZ: [0]    #No variance
physics.producers.generator.SigmaThetaYZ: [0]    #
physics.producers.generator.AngleDist: 0         #0:Uniform 1:Gaussian

(More information on all the settable parameters can be found here) </br></br>To get started using this particle gun, copy the above lines into muon_gun.fcl, save the file, and run the following lines in the command line.

lar -c muon_gun.fcl -n 10 #Generates 10 particles
lar -c standard_g4_sbnd.fcl -s prodsingle_sbnd_SinglesGen-*.root #Runs GEANT4
lar -c standard_detsim_sbnd.fcl -s prodsingle_sbnd_SinglesGen-*_G4-*.root #Runs detector simulation

The output file is an ART file containing LArSoft data products. We can see exactly what data products are in this file by executing

lar -c eventdump.fcl prodsingle_sbnd_SinglesGen-*_G4-*_DetSim-*.root -n 1

We will be looking primarily at MCParticle data products for this tutorial, and we can check that our DetSim file has them by executing

lar -c eventdump.fcl prodsingle_sbnd_SinglesGen-*_G4-*_DetSim-*.root -n 1 | grep MCParticle

where the output shows a vector with the simb::MCParticle class. We will be using this class to get the G4 level information about our muons. First, however, we want to see what our detector sees using the TITUS Event Display

3 - TITUS Event Display

To view our event in TITUS, we should first setup the display by opening a fresh terminal, logging into a gpvm machine, and executing source /sbnd/app/users/mdeltutt/static_evd/ To open the event display for our muon gun, run -s /sbnd/app/users/$USER/tutorial/data/prodsingle_sbnd_SinglesGen-*_G4-*_DetSim-*.root. When it opens, select Raw Digit on the right pane, wait a few seconds, and VOILÀ! You are now seeing what SBND sees when a muon passes through the detector! Clicking next in the top left pane allows you to see all the other muon events and they should all look the same!

4 - Analyzer Module

The SBND Commissioning Page - Get Started Guide covers the hitdumper module, which stores most relevant detector information such as hit and crt readouts. Our focus will be the truth level information, and how to write a simple analyzer module. First, make sure you have setup the LArSoft code space, then execute the following

cd $MRB_SOURCE/sbndcode/sbndcode
mkdir PGUNtutorial
cd PGUNtutorial
cetskelgen analyzer MyAnalyzer

cetskelgen produces a skeleton of a module with the class name being MyAnalyzer. Specifying analyzer after cetskelgen determines the type of art module being used, namely an analyzer module. There are also producer and filter modules which won’t be covered here, but you can familiarize yourself with them here. Before we get too far ahead, we need to add the sub directory to the CMakeLists.txt so that our module will be compiled. To do this execute

echo 'add_subdirectory(PGUNtutorial)' >> ../CMakeLists.txt

and it will add this line to the end of the CMakeLists.txt file. We also need a CMakeLists.txt file in our current directory, so go ahead and make this file and copy the following lines into it

              MODULE_LIBRARIES larcorealg_Geometry
                               larsim_Simulation lardataobj_Simulation
                               ${ART_ROOT_IO_TFILE_SUPPORT} ${ROOT_CORE}
                               art_Persistency_Common canvas
                               art_Persistency_Provenance canvas
                               art_Utilities canvas

You should now have two files, and CMakeLists.txt in the folder $MRB_SOURCE/sbndcode/sbndcode/PGUNtutorial. For reference your module should look something like this in VSCode: Module We will now step through the analyzer module to build a module that can read truth level information and output a TTree root file. For help on making any analyzer module, you can refer to the ultimate analyzer module, the Analysis Tree (which I will be referencing myself to make this module). However, this tutorial will cover the basics so we won’t need to invoke the power of the Analysis Tree yet. In the following sections, I’ll be adding headers, functions, and pointers to data products.

4.1 - MyAnalyzer Module

We are only going to be looking at the MCParticle class, so we’ll need to add the following headers to the top of

//Include these for MCParticle
#include "nusimdata/SimulationBase/MCTruth.h"
//Include these for TTree building
#include "art_root_io/TFileService.h"
#include "art_root_io/TFileDirectory.h"
#include "TTree.h"

After the required functions sections, we’ll add two functions:

//My functions
  void beginJob();
  void reset();

The beginJob() function initializes the variables and TTree, void analyze(art::Event const& e) override pushes data to each TBranch, and void reset() resets the variables after they’ve been written for a single event. Under private:, we’ll initialize our variables and module labels.

  //Simulation information
  TTree*   fEventTree;
  Int_t    run;                  
  Int_t    subrun;               
  Int_t    event;
  //Geant info
  Int_t     no_primaries;                 
  std::vector<Int_t>    pdg;              
  std::vector<Int_t>    status;           
  std::vector<Float_t>  Eng;              
  std::vector<Float_t>  EndE;
  std::vector<Float_t>  Mass;
  std::vector<Float_t>  Px;
  std::vector<Float_t>  Py;
  std::vector<Float_t>  Pz;
  std::vector<Float_t>  P;
  std::vector<Float_t>  StartPointx;
  std::vector<Float_t>  StartPointy;
  std::vector<Float_t>  StartPointz;
  std::vector<Float_t>  StartT;  
  std::vector<Float_t>  EndT;          
  std::vector<Float_t>  EndPointx;
  std::vector<Float_t>  EndPointy;
  std::vector<Float_t>  EndPointz;
  std::vector<Float_t>  theta_xz;    
  std::vector<Float_t>  theta_yz;    
  std::vector<Int_t>    NumberDaughters;
  std::vector<Int_t>    TrackId;
  std::vector<Int_t>    Mother;

  //Module labels
  std::string fMCShowerModuleLabel;
  std::string fMCTrackModuleLabel;

For now, we’ll ignore MyAnalyzer::MyAnalyzer(fhicl::ParameterSet const& p), which passes .fcl parameters into the module. Add the beginJob() function under the fhicl parameter set function.

void MyAnalyzer::beginJob(){
  // Implementation of required member function here.
  std::cout<<"job begin..."<<std::endl;
  art::ServiceHandle<art::TFileService> tfs;
  //Make TTree
  fEventTree = tfs->make<TTree>("Event", "Neutrino interaction info.");

  //Simulation branches
  fEventTree->Branch("event", &event,"event/I");
  fEventTree->Branch("run", &run,"run/I");
  fEventTree->Branch("subrun", &subrun,"subrun/I");

  //Geant info

We should now populate our analyze function, this will be the main piece that pushes the necessary data to the correct branches.

void MyAnalyzer::analyze(art::Event const& e)
  reset(); //Initialize parameters

  //Gets particle information
  art::ServiceHandle<cheat::ParticleInventoryService> pi_serv;

  //Geant info
  const sim::ParticleList& plist = pi_serv->ParticleList();
  sim::ParticleList::const_iterator itPart = plist.begin(), pend = plist.end(); // iterator to pairs (track id, particle)
  for(size_t iPart = 0; (iPart < plist.size()) && (itPart != pend); ++iPart){
    const simb::MCParticle* pPart = (itPart++)->second;
    if (!pPart) {
      throw art::Exception(art::errors::LogicError)
      << "GEANT particle #" << iPart << " returned a null pointer";
    }//endif pPart
  //Geant info
  status.push_back( pPart->StatusCode());
  theta_xz.push_back( std::atan2(pPart->Px(), pPart->Pz()));
  theta_yz.push_back( std::atan2(pPart->Py(), pPart->Pz()));
  }//endfor iPart
}//end analyze

Finally, we’ll complete our module by building the reset function, which clears all of the vectors.

void MyAnalyzer::reset(){
  //Geant info
  StartT.clear()  ;
  EndT.clear()     ;     
  theta_xz.clear()  ;  
  theta_yz.clear()   ; 
}//end reset

4.2 - MyAnalyzer FHiCL —————————————————————————– Great! Now you have an analyzer module capable of creating a TTree for analysis. Let’s make the .fcl file that will allow us to run the module on an art root file. Create and empty file and name it MyAnalyzer.fcl. First we’ll need to attach the service .fcl files, necessary to reference the data products. Add the following lines at the top of MyAnalyzer.fcl

#include "simulationservices_sbnd.fcl"
#include "particleinventoryservice.fcl"
#include "backtrackerservice.fcl"
#include "rootoutput_sbnd.fcl"

process_name: MyAnalyzer

  #Load the service that manages root files for histograms.
  TFileService: { fileName: "MyAnalyzer.root" }
  RandomNumberGenerator: {}                     #ART native random number generator
  FileCatalogMetadata: @local::sbnd_file_catalog_mc
  ParticleInventoryService: @local::standard_particleinventoryservice

The last thing we need is our custom physics list. We only have a single analyzer, so we can simply specify the analyzer module name and be done with writing our .fcl file.

 filters:  {} 
   MyAnalyzer:{module_type: "MyAnalyzer"}

  #define the producer and filter modules for this path, order matters, 
  #filters reject all following items.  see lines starting physics.producers below
  ana:  [ MyAnalyzer]

  #define the output stream, there could be more than one if using filters 
  #stream1:  [ out1 ]

  #trigger_paths is a keyword and contains the paths that modify the art::event, 
  #ie filters and producers
  #trigger_paths: [reco] 

  #end_paths is a keyword and contains the paths that do not modify the art::Event, 
  #ie analyzers and output streams.  these all run simultaneously
  end_paths:     [  ana ] 

Now you have a working analyzer module, and a .fcl file to go along with it! We just have a few more steps before we can start looking at our simulation data. To learn more about .fcl files, consider referring to The ART Workbook section 24.

4.3 - Running MyAnalyzer

Now we should take a moment to check that we have everything before compiling our code. Give a quick ls command and your output should be three files: CMakeLists.txt, MyAnalyzer.fcl and If you have these three, you’re good to go. Open a fresh terminal and run ssh $ to log into the build machine. Setup your working area (same as gpvm) and execute the following:

mrb i -j64

If the stage install is not a success, check that your CMakeLists.txt are properly filled in, and that you have all of the correct header files for both the module and .fcl file. For reference, the complete tutorial for building the analyzer module lives at /sbnd/app/users/brindenc/sbndcode_v09_43_00/srcs/sbndcode/sbndcode/PGUNtutorial/. If the build is successful, head over to the tutorial space you made earlier and run the analyzer module:

cd /sbnd/app/users/$USER/tutorial/data
lar -c $MRB_SOURCE/sbndcode/sbndcode/PGUNtutorial/MyAnalyzer.fcl -s prodsingle_sbnd_SinglesGen-*_G4-*_DetSim-*.root

You should now have MyAnalyzer.root in your data directory. Nice work! Now, click on MyAnalyzer.root by opening your current folder in the explorer pane and selecting the file. It should look something like this when you select the Eng TBranch VSCode2 VSCode3