Take a look at: https://www.cs.rice.edu/~vo9/deeplearning/2019/slides/Rivanna.pdf


Updated Instructions for Generating Monte Carlo Simulations in RUS format for SpinQuest (Last updated on Feb 23, 2025, by Forhad)

If you would like to generate simulated events for SpinQuest experiments and save the output file in RUS format, you should follow these steps.

1. Access the Rivanna computer (https://ood.hpc.virginia.edu/pun/sys/dashboard)

2. Clone the repository: git clone https://github.com/uva-spin/UVA_RUS_Basic

3. Go to the repository:
cd UVA_RUS_Basic

4. Set up the environment:
source /project/ptgroup/spinquest/this-e1039.sh

5. Run the simulation macro locally for testing. In the Fun4Sim.C macro, we have used: 

const bool count_only_good_events = true;
se->run(nevent, count_only_good_events);

This means that the Fun4All macro will keep running until we get accepted events, as required by the SQGeomAcc condition.The legacy generator has been used as the default for J/Psi productions because it is faster compared to Pythia8. You can use the package based on your interest. For example, let's use Drell-Yan events, where the beam interaction point is at the target location:

cd DY_Target
root -b 'Fun4Sim.C(10)'

6. Once the job runs locally and looks alright, you can submit a few jobs on the grid before submitting large jobs:
./jobscript.sh test 2 10

7. Once everything looks alright, submit large jobs. For example:
./jobscript.sh DY_Target 100 1000

8.   Once all the jobs are completed, you can use the hadd or TChain command to merge the ROOT files. Here is an example using the hadd command:

  1. Basic Usage: To merge multiple ROOT files into a single output file, use: hadd output_file.root input_file1.root input_file2.root input_file3.root
  2. Merging Large File Lists with a Minimum File Size Constraint: To merge only root.root files larger than 20 KB from a directory, use:   hadd output_file.root $(find output_file_path -path "*/RUS.root" -size +20k)

Configurations of Event Generation Properly  (Last updated on Jan 30, 2025, by Forhad)

Wall time when dimuons are accepted in a top-bottom plane-wise (default).

Vertex LocationChannelGenerator

Acceptance Conditions
SQGeomAcc::PAIR_TBBT

Total Events

Rivanna HPC Cluster (Wall Time) in sec
TargetDYPythia8True100
1038
TargetDYLegacy E906True1001638
TargetJPsiPythia8True102804
TargetJPsiLegacy E906True101394


Wall time when dimuons are accepted in SpinQuest spectrometer

Vertex LocationChannelGenerator

Acceptance Conditions
SQGeomAcc::PAIR

Total Events

Rivanna HPC Cluster (Wall Time) in sec
TargetDYPythia8True100
593
TargetDYLegacy E906True1001386
TargetJPsiPythia8True101182
TargetJPsiLegacy E906True10945

Here are the updated instruction of generating Monte Carlo for SpinQuest simulations and analysis

            The jobscript.sh should look like this:

#!/bin/bash

dir_macros=$(dirname $(readlink -f $BASH_SOURCE))
jobname=$1
njobs=$2
nevents=$3
echo "njobs=$njobs"
echo "nevents=$nevents"
work=/scratch/$USER/MC/$jobname
mkdir -p $work
chmod -R 01755 $work
cd $dir_macros

for (( id=1; id<=$njobs; id++ ))
do 
  mkdir -p $work/$id/
  chmod -R 01755 $work/$id
  cd $work/$id/
  cp $dir_macros/*.C .
  cp $dir_macros/*.cfg .
  cp $dir_macros/*.txt .
  cp $dir_macros/*.slurm .
  sed -i "s/1234/$nevents/" Fun4Sim.C
  echo "submitting job number = $id"
  sbatch grid.slurm

   Which will submit the job according to the grid.slurm script:

#!/bin/sh
#SBATCH --ntasks=1
#SBATCH --ntasks-per-node=1
#SBATCH --time=6:00:00
#SBATCH --output=slurm.out
#SBATCH --error=slurm.err
#SBATCH --partition=standard
#SBATCH -A spinquest
export DISPLAY=:0.0
root -l Fun4Sim.C
root -l filter.C

Additional Information about changing the Setup

The script that you need to change for any particular configuration is the Fun4Sim.C

Most of this configuration file you will not need to change, but depending on what you want to do here are some pointers

Choice of the Generators: 

This is how you would select what physics generator you want to use (for normal DY or J/Psi use pythia):

const bool gen_pythia8  = true;
const bool gen_cosmic   = false;
const bool gen_gun      = false;
const bool gen_particle = false;
const bool read_hepmc   = false;
const bool gen_e906legacy = false;

gen_pythia8:

   1. If you want to use Pythia8, you can load the config file using this function. pythia8->set_config_file("phpythia8_DY.cfg") or pythia8->set_config_file("phpythia8_Jpsi.cfg"); 
depending on your work.  
   2. In the configuration file, you can change the configuration based on your need.  For example, you can set the mass cut at Pythia8 by putting this line: PhaseSpace:mHatMin = 4.0. 
   3. In addition to the pythia configuration file, you can add more variables using PHPy8ParticleTrigger. A basic setup has shwn below, but you can find more function here.

// pythia8
  if(gen_pythia8) {    
    PHPythia8 *pythia8 = new PHPythia8();
    //pythia8->Verbosity(99);
    pythia8->set_config_file("phpythia8_DY.cfg");
    if(SQ_vtx_gen) pythia8->enableLegacyVtxGen();
    else{
      pythia8->set_vertex_distribution_mean(0, 0, target_coil_pos_z, 0);
    } 
    pythia8->set_embedding_id(1);
    se->registerSubsystem(pythia8);

    pythia8->set_trigger_AND();

    PHPy8ParticleTrigger* trigger_mup = new PHPy8ParticleTrigger();
    trigger_mup->AddParticles("-13");
    //trigger_mup->SetPxHighLow(7, 0.5);
    //trigger_mup->SetPyHighLow(6, -6);
    trigger_mup->SetPzHighLow(120, 10);
    pythia8->register_trigger(trigger_mup);

    PHPy8ParticleTrigger* trigger_mum = new PHPy8ParticleTrigger();
    trigger_mum->AddParticles("13");
    //trigger_mum->SetPxHighLow(-0.5, 7);
    //trigger_mum->SetPyHighLow(6, -6);
    trigger_mum->SetPzHighLow(120, 10);
    pythia8->register_trigger(trigger_mum);
  }
  }

gen_e906legacy:


This generator is adopted from the SeaQuest experiment. If you interested in generating dimuon mass, especially in the low mass regions, you can use the script, because generating 
more statistics using gen_pythia8 is computationally much more expensive.  And If you use gen_e906legacy generator, you will need to save the event weight information, since 
each event has a generator weight.  For example, the DY has been used below:

if(gen_e906dim){
    SQPrimaryParticleGen *e906legacy = new  SQPrimaryParticleGen();
    const bool pythia_gen = false;
    const bool drellyan_gen = true;
    const bool JPsi_gen = false;
    const bool Psip_gen = false;  

    if(drellyan_gen){
      e906legacy->set_xfRange(0.1, 0.5); //[-1.,1.]
      e906legacy->set_massRange(3.0, 10.0);
      e906legacy->enableDrellYanGen();
    }
    if(Psip_gen){ 
      e906legacy->set_xfRange(0.1, 0.5); //[-1.,1.]
      e906legacy->enablePsipGen();
    }
    if(JPsi_gen){
      e906legacy->set_xfRange(0.1, 0.5); //[-1.,1.]
      e906legacy->enableJPsiGen();
    }

}



Generating Single or Multi muons from using the multi particle gun

For making a single  or multiple muons you can use multi particle gun. The portion of the script below can be found on the SimChainDev module. Set the gen_particle true if you would like to generate single muons from multi particle guns. Some user function are available on the doxygen page.

if(gen_particle) {

    PHG4SimpleEventGenerator *genp = new PHG4SimpleEventGenerator("MUP");

    //genp->set_seed(123);
    genp->add_particles("mu+", nmu);  // mu+,e+,proton,pi+,Upsilon
    if (legacyVtxGen) genp->enableLegacyVtxGen();
    else{
    genp->set_vertex_distribution_function(PHG4SimpleEventGenerator::Uniform,
        PHG4SimpleEventGenerator::Uniform,
        PHG4SimpleEventGenerator::Uniform);
    genp->set_vertex_distribution_mean(0.0, 0.0, -300);
    genp->set_vertex_distribution_width(100.0, 100.0, 400.0);
    genp->set_vertex_size_function(PHG4SimpleEventGenerator::Uniform);  [For Gaussian distributions use "Gaus"]
    genp->set_vertex_size_parameters(0.0, 0.0);

    }

Basic Configurations for the simulations:

const double target_coil_pos_z = -300;
  const int nmu = 1;
  int embedding_opt = 0;
  const bool legacy_rec_container = true;

  const bool do_collimator = true;
  const bool do_target = true;
  const bool do_e1039_shielding = true;
  const bool do_fmag = true;
  const bool do_kmag = true;
  const bool do_absorber = true;
  const bool do_dphodo = true;
  const bool do_station1DC = false;       //station-1 drift chamber should be turned off by default

  const double target_l = 7.9; //cm
  const double target_z = (7.9-target_l)/2.; //cm
  const int use_g4steps = 1;

  const double FMAGSTR = -1.044;
  const double KMAGSTR = -1.025;


Here you can change the muon to mu+ or mu-, or change the or change the vertex interaction volume, for example this one creates a cylinder 1m radius around the beamline from 700 cm upstream of the dump to 100 cm inside the dump.


 

To save the accepted tracks/Dimuons: 

If you would like to save only accepted single muon or dimuon, you can use SQGeomAcc.

  /// Save only events that are in the geometric acceptance.
  SQGeomAcc* geom_acc = new SQGeomAcc();
  geom_acc->SetMuonMode(SQGeomAcc::PAIR); // PAIR, PAIR_TBBT, SINGLE, SINGLE_T, etc.
  geom_acc->SetPlaneMode(SQGeomAcc::HODO_CHAM); // HODO, CHAM or HODO_CHAM
  geom_acc->SetNumOfH1EdgeElementsExcluded(4); // Exclude 4 elements at H1 edges
  se->registerSubsystem(geom_acc);