Basic rat macros and code

README.md

.cc Files

To run a file called Bob.cc on the rat created root file jim.root do the following

root
.L Bob.cc+
Bob( jim.root )

.py Files

To run a file called frank.py on the rat created root file jim.root do the following

python frank.py jim.root

AirFill.mac

# To use an air filled detector add this command to a macro before the /run/initialize
/rat/db/set DETECTOR geo_file "geo/snoplus_air.geo"

# To just use an air filled inner volume of the AV (scint volume) add this instead
/rat/db/set GEO[scint] material "air" 

autowms.conf

[
OutputStorage  =  "/tmp/jobOutput";
JdlDefaultAttributes =  [
    RetryCount  =  3;
    rank  = - other.GlueCEStateEstimatedResponseTime;
    PerusalFileEnable  =  false;
    AllowZippedISB  =  true;
    requirements  =  other.GlueCEStateStatus == "Production";
    ShallowRetryCount  =  10;
    SignificantAttributes  =  {"Requirements", "Rank", "FuzzyRank"};
    MyProxyServer  =  "lcgrbp01.gridpp.rl.ac.uk";
    ];
ErrorStorage  =  "/tmp";
VirtualOrganisation  =  "snoplus.snolab.ca";
ListenerStorage  =  "/tmp";
WMProxyEndpoints  =  {"https://lcgwms03.gridpp.rl.ac.uk:7443/glite_wms_wmproxy_server","https://lcgwms02.gridpp.rl.ac.uk:7443/glite_wms_wmproxy_server","https://wms01.grid.hep.ph.ic.ac.uk:7443/glite_wms_wmproxy_server","https://wms02.grid.hep.ph.ic.ac.uk:7443/glite_wms_wmproxy_server"};
WorkloadManagerProxy =  [
	maxInputSandboxFiles = 20;
	];

]

Basic.mac

# To simulate electrons in a basic macro
/run/initialize
/rat/proc frontend
/rat/proc trigger
/rat/proc eventbuilder
/rat/proc calibratePMT
# Save the output
/rat/proc outroot
/rat/procset file "Speedy.root"

# Generate events
/generator/add combo gun:fill
/generator/vtx/set e- 0 0 0 1.0
/generator/pos/set 0.0 0.0 0.0
/generator/rate/set 1
/run/beamOn 10
exit

DAQ.mac

# The DAQ Processor loop, should be added to all macros that you wish to analyses the hit information
/rat/proc frontend # Front end electrons simulation
/rat/proc trigger # Trigger simulation
/rat/proc eventbuilder # Event builder (combine hits into an event)
/rat/proc calibratePMT # Calibrate the PMT hit times and charges

DiscModel.mac

# To change to the disc optical model add these commands before /run/intialize
/rat/db/set GEO[innerPMT] grey_disc [1]
/rat/db/set GEO[outerPMT] grey_disc [1]

DoubleBeta.mac

# To generate neutrinoless double beta decay events, use this generator
/generator/add combo doublebeta:fill
/generator/vtx/set Nd150 0
/generator/pos/set 0 0 0
/generator/rate/set 1

# To generate neutrino present double beta decay events, use this generator
/generator/add combo doublebeta:fill
/generator/vtx/set Nd150 2 
/generator/pos/set 0 0 0
/generator/rate/set 1

gangaJob.py

j=Job()
j.backend='LCG'
j.application='RATUser'
j.application.outputDir='RATUser/tutorialTest'
j.application.ratMacro='run_electrons_snoplus.mac'
j.application.ratVersion='1818a6b5853a930c56f1b9acd73107b6bcb7a796'
j.application.ratBaseVersion='4'
j.submit()

print 'YOUR GANGA JOBID: ',j.fqid
print 'CHECK WITH j.status and j.peek()'

hello.jdl

Executable = "hello.scr";
Arguments = "";
StdOutput = "hw.out";
StdError = "hw.err";
InputSandbox = {"hello.scr"};
OutputSandbox = {"hw.out", "hw.err"};
VirtualOrganisation = "snoplus.snolab.ca";

hello.scr

#!/bin/bash
echo "Hello World!  SNO+ is on the grid."

La138.mac

# To simulate 138La decays use this generator
/generator/add decaychain 138La:fill:uniform
/generator/pos/set 0 0 0
/generator/rate/set 1

# Alternatively to add screening corrections (nuclear model screening) use this line instead of the first above
/generator/add decaychain 138La:fill:exponential::screen

Load.cc

////////////////////////////////////////////////////////
/// My preferred method to load a root file in root.
/// Must load in root as .L Load.cc+
///
/// 11/10/2012 : New File
////////////////////////////////////////////////////////

#include <RAT/DS/Root.hh>
#include <RAT/DS/EV.hh>
#include <RAT/DS/Run.hh>

#include <TH1D.h>
#include <TTree.h>
#include <TFile.h>

#include <time.h>
using namespace std;

void
LoadRootFile( const char* lpFile,
              TTree **tree,
              RAT::DS::Root **rDS,
              RAT::DS::Run **rRun );

void 
DoSomething( const char* inFile )
{
  // Load the root file first
  RAT::DS::Root* rDS;
  RAT::DS::Run* rRun;
  TTree *tree;
  LoadRootFile( inFile, &tree, &rDS, &rRun );
  
  time_t codeStart = time( NULL );

  for( int iEvent = 0; iEvent < tree->GetEntries(); iEvent++ )
    {
      if( iEvent % 100 == 0 )
        cout << iEvent << " finished at " << time( NULL ) - codeStart << endl;
      tree->GetEntry( iEvent );
      // Do something here
    }
}

void
LoadRootFile( const char* lpFile,
              TTree **tree,
              RAT::DS::Root **rDS,
              RAT::DS::Run **rRun )
{
  TFile *file = new TFile( lpFile );
  (*tree) = (TTree*)file->Get( "T" );
  TTree *runTree = (TTree*)file->Get( "runT" );

  *rDS = new RAT::DS::Root();

  (*tree)->SetBranchAddress( "ds", &(*rDS) );

  *rRun = new RAT::DS::Run();

  runTree->SetBranchAddress( "run", &(*rRun) );
  runTree->GetEntry();
}

NoMuHadronic.mac

# Turn off muon and hadronic physics add these before /run/initialize
/glg4debug/glg4param omit_muon_processes 1.0
/glg4debug/glg4param omit_hadronic_processes 1.0

plot_nhit.py

#!/usr/bin/env python
import ROOT
import rat
import sys

nhit_plot = ROOT.TH1D("nhit", "nhit", 3000, 0, 3000)
for ds in rat.dsreader(sys.argv[1]):
    for iEV in range(0, ds.GetEVCount()):
        nhit_plot.Fill(ds.GetEV(iEV).GetNhits())

nhit_plot.Draw()
raw_input("RET to exit")

plot_pmt_hit_times.py

#!/usr/bin/env python
import ROOT
import rat
import sys

# PMT Cal is always between 0 and 500 ns
hit_times = ROOT.TH1D("hitTimes", "hitTimes", 500, 0.0, 500.0);
for ds in rat.dsreader(sys.argv[1]):
    for iev in range(ds.GetEVCount()):
        ev = ds.GetEV(iev)
        for ipmt in range(ev.GetPMTCalCount()):
            hit_times.Fill(ev.GetPMTCal(ipmt).GetTime())
hit_times.Draw()
raw_input("RET to EXIT")

plot_pmt_pos.py

#!/usr/bin/env python
import ROOT
import rat
import sys


plot = ROOT.TGraph2D()
plot_point = 0
# run is the run tree entry and ds is the TTree entry
run = rat.runreader(file_name)
pmt_prop = run.GetPMTProp()
# Loop over the pmts
for ipmt in range(pmt_prop.GetPMTCount()):
    plot.SetPoint(plot_point, pmt_prop.GetPos(ipmt).x(), pmt_prop.GetPos(ipmt).y(), pmt_prop.GetPos(ipmt).z())
    plot_point += 1

plot.Draw("A*")
raw_input("RET to EXIT");

PlotNhit.cc

#include <RAT/DS/Root.hh>
#include <RAT/DS/EV.hh>
#include <RAT/DS/Run.hh>

#include <TH1D.h>
#include <TTree.h>
#include <TFile.h>

#include <time.h>
using namespace std;

void
LoadRootFile( const char* lpFile,
              TTree **tree,
              RAT::DS::Root **rDS,
              RAT::DS::Run **rRun );

void 
PlotNhit( const char* inFile )
{
  // Load the root file first
  RAT::DS::Root* rDS;
  RAT::DS::Run* rRun;
  TTree *tree;
  LoadRootFile( inFile, &tree, &rDS, &rRun );
  
  time_t codeStart = time( NULL );
  TH1D* nhitPlot = new TH1D( "nhit", "nhit", 3000, 0, 3000 );

  for( int iEvent = 0; iEvent < tree->GetEntries(); iEvent++ )
    {
      if( iEvent % 100 == 0 )
        cout << iEvent << " finished at " << time( NULL ) - codeStart << endl;
      tree->GetEntry( iEvent );
      for( int iEV = 0; iEV < rDS->GetEVCount(); iEV++ )
        nhitPlot->Fill( rDS->GetEV( iEV )->GetNhits() );
    }
}

void
LoadRootFile( const char* lpFile,
              TTree **tree,
              RAT::DS::Root **rDS,
              RAT::DS::Run **rRun )
{
  TFile *file = new TFile( lpFile );
  (*tree) = (TTree*)file->Get( "T" );
  TTree *runTree = (TTree*)file->Get( "runT" );

  *rDS = new RAT::DS::Root();

  (*tree)->SetBranchAddress( "ds", &(*rDS) );

  *rRun = new RAT::DS::Run();

  runTree->SetBranchAddress( "run", &(*rRun) );
  runTree->GetEntry();
}

PlotPMTHitTimes.cc

#include <RAT/DS/Root.hh>
#include <RAT/DS/EV.hh>
#include <RAT/DS/Run.hh>

#include <TH1D.h>
#include <TTree.h>
#include <TFile.h>

#include <time.h>
using namespace std;

void
LoadRootFile( const char* lpFile,
              TTree **tree,
              RAT::DS::Root **rDS,
              RAT::DS::Run **rRun );

void 
PlotPMTHitTimes( const char* inFile )
{
  // Load the root file first
  RAT::DS::Root* rDS;
  RAT::DS::Run* rRun;
  TTree *tree;
  LoadRootFile( inFile, &tree, &rDS, &rRun );
  
  time_t codeStart = time( NULL );

  // PMT Cal is always between 0 and 500 ns
  TH1D* hitTimes = new TH1D( "hitTimes", "hitTimes", 500, 0.0, 500.0 );

  for( int iEvent = 0; iEvent < tree->GetEntries(); iEvent++ )
    {
      if( iEvent % 100 == 0 )
        cout << iEvent << " finished at " << time( NULL ) - codeStart << endl;
      tree->GetEntry( iEvent );
      // Loop over the triggered events
      for( int iEV = 0; iEV < rDS->GetEVCount(); iEV++ )
        {
          RAT::DS::EV* rEV = rDS->GetEV( iEV );
          // Loop over the PMTCals
          for( int ipmt = 0; ipmt < rEV->GetPMTCalCount(); ipmt++ )
            {
              RAT::DS::PMTCal* rPMTCal = rEV->GetPMTCal( ipmt );
              hitTimes->Fill( rPMTCal->GetTime() );
            }
        }
    }
  hitTimes->Draw();
}

void
LoadRootFile( const char* lpFile,
              TTree **tree,
              RAT::DS::Root **rDS,
              RAT::DS::Run **rRun )
{
  TFile *file = new TFile( lpFile );
  (*tree) = (TTree*)file->Get( "T" );
  TTree *runTree = (TTree*)file->Get( "runT" );

  *rDS = new RAT::DS::Root();

  (*tree)->SetBranchAddress( "ds", &(*rDS) );

  *rRun = new RAT::DS::Run();

  runTree->SetBranchAddress( "run", &(*rRun) );
  runTree->GetEntry();
}

PlotPMTPos.cc

#include <RAT/DS/Root.hh>
#include <RAT/DS/EV.hh>
#include <RAT/DS/Run.hh>

#include <TH1D.h>
#include <TTree.h>
#include <TFile.h>

#include <time.h>
using namespace std;

void
LoadRootFile( const char* lpFile,
              TTree **tree,
              RAT::DS::Root **rDS,
              RAT::DS::Run **rRun );

void 
PlotPMTPos( const char* inFile )
{
  // Load the root file first
  RAT::DS::Root* rDS;
  RAT::DS::Run* rRun;
  TTree *tree;
  LoadRootFile( inFile, &tree, &rDS, &rRun );
  
  time_t codeStart = time( NULL );

  TGraph2D* pmtPositions = new TGraph2d();
  int graphPoint = 0;
  // Get the PMTProperties class instance
  RAT::DS::PMTProperties* pmtProp = rRun->GetPMTProp();
  // Loop over the PMTs
  for( int ipmt = 0; ipmt < pmtProp->GetPMTCount(); ipmt++ )
    pmtPositions->SetPoint( graphPoint++, pmtProp->GetPos( ipmt ).x(), pmtProp->GetPos( ipmt ).y(), pmtProp->GetPos( ipmt ).z() );
  // Draw
  pmtPositions.Draw("A*");
}

void
LoadRootFile( const char* lpFile,
              TTree **tree,
              RAT::DS::Root **rDS,
              RAT::DS::Run **rRun )
{
  TFile *file = new TFile( lpFile );
  (*tree) = (TTree*)file->Get( "T" );
  TTree *runTree = (TTree*)file->Get( "runT" );

  *rDS = new RAT::DS::Root();

  (*tree)->SetBranchAddress( "ds", &(*rDS) );

  *rRun = new RAT::DS::Run();

  runTree->SetBranchAddress( "run", &(*rRun) );
  runTree->GetEntry();
}

PosGenerators.mac

# To generate events at a set point use the point PosGenerator
/generator/add combo gun:point
/generator/pos/set 0.0 0.0 0.0

# To generate events uniformly within a volume use the fill PosGenerator
/generator/add combo gun:fill
/generator/pos/set 0.0 0.0 0.0 # Fills the volume which contains this point
/generator/pos/set scint # Fills the volume called scint

# To generate events uniformly on the surface of a volume use the paint PosGenerator
/generator/add combo gun:paint
/generator/pos/set 0.0 0.0 0.0 # Paints the inner and outer surfaces of the volume which contains this point
/generator/pos/set scint # Paints the inner and outer surfaces of the volume called scint
/generator/pos/set scint 0- # Paints the inner surface with a thickness of 0 mm (on surface) of the volume called scint
/generator/pos/set scint 1+ # Paints the outer surface with a thickness of 1 mm of the volume called scint

run_electrons_snoplus.mac

#//////////////////////////////////////////////////////////////////
# Last svn revision: $Id$ 
#//////////////////////////////////////////////////////////////////
#/**
#*       Example Macro for simulating electrons in SNO+
#*       
#*       run_electrons_snoplus.mac
#*       Date:           01/12/2008
#*       contact:        J. Wilson, Oxford
#*/

# Don't care about hadrons or muons so quicker not to initialise these processes
/glg4debug/glg4param omit_muon_processes  1.0
/glg4debug/glg4param omit_hadronic_processes  1.0

# Select the snoplus geometry file (should be default)
/rat/db/set DETECTOR geo_file "geo/snoplus.geo"

/run/initialize

# BEGIN EVENT LOOP
/rat/proc frontend
/rat/proc trigger
/rat/proc eventbuilder
/rat/proc count
/rat/procset update 10
/rat/proc calibratePMT
/rat/proclast outroot
/rat/procset file "snoplus_electrons_labppo.root"
# END EVENTLOOP

# Use the simple gun vertex generator, combined with the  fill position generator and poisson timing
/generator/add combo gun:fill:poisson
# Select particle type and momentum direction (0,0,0 = isotropic) and energy in MeV 
/generator/vtx/set e- 0 0 0 1.0
# Define a position in the centre of the scintillator region to select this one to be filled.
/generator/pos/set 0.0 0.0 0.0
# Rate of 1 per second
/generator/rate/set 1

# The above simulates electrons in the scintillator, if we want them in the AV and H2O too add more generators
# (one for each defined geometry region) to list
/generator/add combo gun:fill:poisson
/generator/vtx/set e- 0 0 0 5.0
# Define a position in the AV to select this one to be filled.
/generator/pos/set 5999.0 0.0 0.0
# Rate of 0.028 per second (scale by relative volume to scintillator region = (6000^3-5945^3)/(5945^3))
# Not quite exact since we didn't consider the neck though
/generator/rate/set 0.028
#/generator/rate/set 0.

/generator/add combo gun:fill:poisson
/generator/vtx/set e- 0 0 0 1.0
# Define a position in the H2Osub region to select this one to be filled.
/generator/pos/set 6050.0 0.0 0.0
# Rate of 1.693 per second (scale by relative volume to scintillator region = (8300^3-6000^3)/(5945^3))
# Not quite exact since we didn't consider the neck though
/generator/rate/set 1.693
#/generator/rate/set 0.0

/generator/add combo gun:fill:poisson
/generator/vtx/set e- 0 0 0 1.0
# Define a position in the H2O region (the bit that contains the PMTs) to select this one to be filled.
/generator/pos/set 6050.0 0.0 0.0
# Rate of 0.306 per second (scale by relative volume to scintillator region = (8600^3-8300^3)/(5945^3))
# Not quite exact since we didn't consider the neck though
/generator/rate/set 0.306
#/generator/rate/set 0.0

# generate 10 events
/run/beamOn 10
exit

ScintFitter.mac

# To the scintFitter over events just add this line to the processor loop
/rat/proc scintFitter

Speedy.mac

# To simulate events more quickly add these lines before the /run/initialize
# Turn off muon and hadronic physics (simulating electrons so this is ok)
/glg4debug/glg4param omit_muon_processes 1.0
/glg4debug/glg4param omit_hadronic_processes 1.0

# Change to the simple detector geometry (no complex ropes, belly plates etc...)
/rat/db/set DETECTOR geo_file "geo/snoplus_Simple.geo"

# Thin the number of photons by a factor of 3.5
/rat/db/set MC thin_factor 3.5

# Change to the disc optical model
/rat/db/set GEO[innerPMT] grey_disc [1]
/rat/db/set GEO[outerPMT] grey_disc [1]

# Standard DAQ processor loop follows
/run/initialize
/rat/proc frontend
/rat/proc trigger
/rat/proc eventbuilder
/rat/proc calibratePMT

# Prune some less interesting information i.e. the mc side
/rat/proc prune
/rat/procset prune "mc.particle,mc.hit,mc.pmt.photon"

# Save the output
/rat/proc outroot
/rat/procset file "Speedy.root"

# Generate events
/generator/add combo gun:fill
/generator/vtx/set e- 0 0 0 1.0
/generator/pos/set 0.0 0.0 0.0
/generator/rate/set 1
/run/beamOn 10
exit

write_track_history.py

#!/usr/bin/env python
import ROOT
import rat
import sys

for ds in rat.dsreader(sys.argv[1]):
    mc = ds.GetMC()
    track_count = mc.GetMCTrackCount()
    for track_num in range(0,track_count):
        mc_track = mc.GetMCTrack(track_num)
        print mc_track.GetTrackID()
        for track_step in range(0,mc_track.GetMCTrackStepCount()):
            mc_track_step = mc_track.GetMCTrackStep(track_step)
            print "\t", mc_track_step.GetProcess()

WriteTrackHistory.cc

#include <RAT/DS/Root.hh>
#include <RAT/DS/MC.hh>
#include <RAT/DS/MCTrack.hh>
#include <RAT/DS/MCTrackStep.hh>
#include <RAT/DS/Run.hh>

#include <TH1D.h>
#include <TTree.h>
#include <TFile.h>

#include <time.h>
using namespace std;

void
LoadRootFile( const char* lpFile,
              TTree **tree,
              RAT::DS::Root **rDS,
              RAT::DS::Run **rRun );

void 
WriteTrackHistory( const char* inFile )
{
  // Load the root file first
  RAT::DS::Root* rDS;
  RAT::DS::Run* rRun;
  TTree *tree;
  LoadRootFile( inFile, &tree, &rDS, &rRun );
  
  time_t codeStart = time( NULL );

  for( int iEvent = 0; iEvent < tree->GetEntries(); iEvent++ )
    {
      if( iEvent % 100 == 0 )
        cout << iEvent << " finished at " << time( NULL ) - codeStart << endl;
      tree->GetEntry( iEvent );
      RAT::DS::MC* rMC = rDS->GetMC();
      for( int iMCTrack = 0; iMCTrack < rMC->GetMCTrackCount(); iMCTrack++ )
        {
          RAT::DS::MCTrack* rMCTrack = rMC->GetMCTrack( iMCTrack );
          cout << rMCTrack->GetTrackID() << endl;
          for( int iMCTrackStep = 0; iMCTrackStep < rMCTrack->GetMCTrackStepCount(); iMCTrackStep++ )
            cout << "\t" << rMCTrack->GetMCTrackStep( iMCTrackStep )->GetProcess() << endl;
        }
    }
}

void
LoadRootFile( const char* lpFile,
              TTree **tree,
              RAT::DS::Root **rDS,
              RAT::DS::Run **rRun )
{
  TFile *file = new TFile( lpFile );
  (*tree) = (TTree*)file->Get( "T" );
  TTree *runTree = (TTree*)file->Get( "runT" );

  *rDS = new RAT::DS::Root();

  (*tree)->SetBranchAddress( "ds", &(*rDS) );

  *rRun = new RAT::DS::Run();

  runTree->SetBranchAddress( "run", &(*rRun) );
  runTree->GetEntry();
}