Example 1: Writing Selected MC events

Suppose you want to do the following:-

• Run a standard Monte Carlo.

• Histogram the number of PMT hits.

• Output any event with more than 40 hits.

how do you go about doing it?

As the memory requirements are very modest you don't need to extend any ZEBRA memory so can leave USER_DEFINE_MEMORY_SIZE alone. The first routine to change in USER_DEFINE_TASK, and the first question is what SUs are required? For a complete list of SUs see section 11.7, but for now, 3 are needed:- MCO Monte Carlo.
OUT Event Output.
HBK HBOOK Support

MCO and OUT are understandable enough, after all we want to run the Monte Carlo and output events, but why HBK?. All HBK does is to initialise and terminate HBOOK, and this could easily be done in the user code. However the SU ANL (Event Analyser) also uses HBOOK which could lead to a conflict; which is to initialise HBOOK? In this application ANL is not being used, but suppose the code continues to be developed and that eventually ANL is needed. So this is at least a potential problem, and the solution is to let neither do the job but for both to request HBK. So HBOOK gets initialised once if, and only if, it is required. The USER_DEFINE_TASK looks like this:-

        subroutine user_define_task
 
        include  'shl_com.inc'
        include  'su_mnemonics.inc'
  
        call su_add( KSU_MCO )
        call su_add( KSU_OUT )
        call su_add( KSU_HBK )
 
        cuser_pro_seq = 'MCO OUT END'
 
        end

There are several important points to appreciate here:-

• Setting cuser_pro_seq (passed via the /SHL_COMx/ common) forces execution of the processors in the order given and prevents the initialisation code from processing the JOBP titles bank or any event loop program. However it will still look for a JOB bank.

• SNOMAN will need more than the 3 SUs shown above. For example MCO requires GEO (Detector Geometry). However the SU control knows of this dependency and will ensure that all required SUs are activated. Your job as the user is to list the processors and support systems you will explicitly need. Strictly speaking this means that you should request ZBR (Zebra) and MTT (Management of Titles), but these are core SUs; they always get loaded.

• MCO and OUT are processors. Activating them via calls to su_add means that they will get initialised and terminated. The point at which they are called in the execution phase is determined by the processing sequence. HBK is a support system, so it is not included in the processing sequence. Execution depends on the SUs or user code that use it.

So, to summarise, the USER_DEFINE_TASK has set up an event loop consisting of Monte carlo and Event Output. HBOOK will be initialised and the JOBP bank or any defined event loop program will be ignored. It is not actually necessary to include calls to su_add for the processors specified in cuser_pro_seq as SNOMAN knows that they must be activated.

Now you are ready to start the main job, writing the USER routine. Your USER has to do 3 things:-

1. Book the required histogram.

2. For each event, fill the histogram and (somehow) filter events with more than 40 pmt hits.

3. Output the histogram.

All this could be accomplished within USER but, no doubt, following good programming practice you would rather do this in 3 routines: USER_INI, USER_EXE and USER_TRM and make USER just a dispatch routine to these. Where should the USER interrupt the program flow? It is best to call USER_INI just before the SHL (Shell) SU processing is called. The shell is the overall program framework so, once processing is about to begin, i.e. its initialisation is complete and processing is starting all required SUs will have been initialised and the user code now has everything ready for it. Note: before 3.00, the suggested position for USER initialisation was after SHL initialisation but, due to structural changes, all SUs now receive several initialisation calls, the very first of which just asks it to confirm that it is non-dummy, see section 11.6.3. So adding user initialisation code after this first initialisation means it would occur far too soon.

The best place for USER_TRM is just before SHL termination starts, so all SUs are still active.

The USER_EXE is a bit more problematic. It has to be called in the main event loop, which just consists of MCO and OUT. Clearly it needs to come between them so could either be called after the MCO process call or before the OUT process call. Both are equally as good so assume the choice to be after MCO. The user routine looks like this:-

      subroutine user( i_su_id ,
     +                 i_process ,
     +                 i_call_status ,
     +                 skip_call )
 
 
      include  'su_com.inc'
 
      integer  i_su_id , i_process , i_call_status
      logical  skip_call
 
      skip_call =.false.
 
      if (             i_su_id .eq. KSU_SHL
     +    .and.      i_process .eq. KSU_PROCESS
     +    .and.  i_call_status .eq. KSU_BEFORE_CALL  ) call user_ini
 
      if (             i_su_id .eq. KSU_MCO
     +    .and.      i_process .eq. KSU_PROCESS
     +    .and.  i_call_status .eq. KSU_AFTER_CALL  ) call user_exe
 
      if (             i_su_id .eq. KSU_SHL
     +    .and.      i_process .eq. KSU_TERMINATE
     +    .and.  i_call_status .eq. KSU_BEFORE_CALL ) call user_trm
 
      end

Setting skip_call false ensures that the program flow is not interrupted.

The USER_INI is very straight forward, it simply opens an HBOOK RZ file and books the histogram:-

      subroutine user_ini
 
      implicit none
 
      integer istat
 
      call hropen(91,'myfile','hbook.hst','N',1024,istat)
      call hbook1(10,'Number of PMTs',100,0.,100.,0.)
 
      end

The USER_EXE is where most of the action takes place. You should only assume that the LMAST is set up, and it does no harm even to check that. There are complete set of protected links on z.inc you can use to navigate round the data structure. Here is the routine:-

        subroutine user_exe
 
        include 'bank_mast.inc'
        include 'bank_ev.inc'
        include 'z.inc'
        include 'zunit.inc'
 
        integer   num_hits
 
*       Check data structure exists.
 
        if ( lmast .eq. 0 ) then
          write(iqlog,*) 'USER_EXE: No event in memory'
          return
        endif
 
*       Loop over all EVs counting all hits.
 
        num_hits = 0
        lev = lq(lmast - KMAST_EV)
 
        do while (lev .ne. 0)
          num_hits = num_hits + iq(lev + KEV_NPM)
          lev = lq(lev)
        enddo
 
*       Fill histogram.
 
        call hfill( 10 , float(num_hits) , 0. , 1. )
 
*       Quit event if not enough hits.
 
        if (num_hits .le. 40 ) call ztell( 800 , 1 )
 
        end

If you are not familiar with code to access the data structure you should look at the examples in:-

...snoman/n_nn/doc/banks/data_structure.tex

The code is simple enough except the call to ZTELL. Calling ZTELL with a second argument of 1 forces it to quit the current event and go onto the next. So if there are not enough hits the event is abandoned and not output. The first number is a ZTELL ID. Codes in the range 600 to 899 are reserved to users. 600 to 799 are errors and will generate an error message, 800 to 899 are signals to quit. When SNOMAN exits you will get a list of all the ZTELL IDs and the number of each that occurred. If you wish you can record an ID without losing control:-

call ztell( 801 , 0 )

increments the ID for 801 but returns to the caller. In case you were wondering, you can abort the entire job by passing a value of 2 for the second argument to ZTELL. For more information about ZTELL see Error Handling in the Programmer's Manual.

All that USER_TRM has to do is to output the histograms and close the file:-

      subroutine user_trm
 
      implicit none
 
      integer  icycle
 
c     Output histograms to file.
      call hrout(0,icycle,' ')
      call hrend('myfile')
 
      end

Having coded and compiled all the user routines, the complete program can be assembled by typing:-

$SNO_TOOLS/link_snolib.scr "*.o" ex1_snoman.exe
(VMS: COPY *.OBJ OBJS.TMP)
(VMS: @SNO_TOOLS:LINK_SNOLIB OBJS.TMP EX1_SNOMAN.EXE)

To complete the job you have to create a command file, say my_run.cmd, containing something like the following commands:-

$event_limit 10
FILE OUT 1 'my_event_output.dat'
@run

This sets a limit of 10 events to process, defines the output file name, and then performs a standard run loading a complete set of titles banks. Now run the program and type:-

@my_run.cmd