Hands on 4: A first physics measurments

• Introduction:
• User Actions II:
  • Exercise 1 Step 1
  • Exercise 1 Step 2
  • Exercise 1 Step 3
  • Exercise 1 Step 4
  • Exercise 2
• Using g4analysis to store results
  • Exercise 3 Step 1 (Optional)
  • Exercise 3 Step 2 (Optional)
  • Exercise 3 Step 3 (Optional)
  • Exercise 3 Step 4 (Optional)
  • Exercise 3 Step 5 (Optional)
• Physics Lists
  • Exercise 4.a
  • Exercise 4.b

Back to Agenda

Introduction

Important Note: The code of this tutorial is an adaptation of Geant4 example B5. Thus you can review almost all concepts from this tutorial in the example that can be found under: <g4-source-tree>/examples/basic/B5.

In this third hands-on you will learn:

• How to define a user run to accumulate quantities over events.
• How to process this information to make a simple physics measurement using G4UderRunAction.
• How to store physics results on output files with g4analysis for further analysis.
• How to modify physics list using constructors.

Follow the instructions of Hands On 1 to configure with cmake the example and build it.
Try out the application:

$ source <where-G4-was-installed>/bin/geant4.[c]sh   $ mkdir build-HandsOn4   $ cmake -DGeant4_DIR=<tutorial>/lib/Geant4-10.0.1 <tutorial>/HandsOn4   $ [g]make [-jN]   $ ./SlacTut

The geometry is the same obtained at the end of the previous hands on.

We will not modify the geometry anymnore, also several sensitive detectors and hits classes have been added to the setup. Take a moment to look at the classes wich name ends in SD and Hit. In particular it is important that you understand how the calorimeter hits work. In the exercise number 2 you will calculate a very simple physics quantity (a partial shower shape) from the energy released in calorimeters.
The goal of these exercises is to show how to interact with Geant4 kernel to extract physics quantities. In complex applications you will probably rely on experimental framework to implement the analysis and recording of data. Additionaly you may need digitization (i.e. simulation of detector read-out) and interface to persistency libraries for data storage.
Geant4 does not provide/recommend specific utilities or for these opeartions because these are strongly user-dependentent. However we provide light-weight histogramming and ntuple utilities. These are compatible with AIDA and ROOT output format (they do not require any library installed on the system). They can also dump ntuples in tabular form in text files (CSV) to import data in (virtually) any analysis system (pylab, R, Octave, Excel, Matlab, Mathematica, ...). If you do not have neither an AIDA-compliant tool or ROOT installed on your system you will not be able to display histograms, but you will still be able to read ntuples written in CSV format.

For your interest here are some links (my personal preferences):

• GNU Octave (open-source MATLAB-inspired)
• R (statistical analysis, de facto standard?)
• ROOT (HEP oriented)
• Java Analysis Studio 3 (AIDA compliant tool with additional support for ROOT files, HEP orineted)
• Better Analysis Workflow with [I]Python (HEP oriented, python based, includes how to read from ROOT files TTrees).

User Actions II

In Exercise 3 of Hands On 3 you have printed on screen, for each simulated event, the hits collected in the hodosope. In this exercise we will show how to accumualte some information (the energy deposited in the calorimeters) over the entire run. We will also show how to merge (e.g. reduce, combine) the results in a multi-threaded application.

Goal of these two exercises is to calculate the average energy released in the electromagnetic and hadronic calorimeter and the average partial shower shape (a shower shape is a quantity that somehow describes the charactersitics of spatial dimensions of particle showers in calorimeters. In this example we will calcualte the fraction of energy released in the electromagnetic calorimeter. These quantities are useful to determine properties of the impinging particle. For example an electron or gamma has the em fraction very close to 1, while an hadron will have a smaller em fraction a muon will have even a smaller value. It is possible to develop algorithms to identify the impinging particle from these quantities. Clearly this is an over-simplified example...).

During Exercise 1 you will modify the application to accumulate the energy released in calorimeters in each event. You will modify 2 files: Run.hh and Run.cc, implementing a user-defined G4Run/

During Exercise 2 you will modify the file RunAction.cc that implements the user defined G4UserRunAction. You willretrieve the information collected in the first exercise and dump on screen the results of your data analysis: energy in calorimeters and shower shape.

During the simulation an instance of a G4Run exists and is managed by Geant4 kernel. User can extend this class to accumulate user data.

Exercise 1 Step 1

Create a user-defined run class

Modify the file Run.hh that defines a class inheriting from G4Run. Extend the class to contain the information to be stored: the total energy deposited in calorimeters and the accumulated shower shape (all of double type). Since you will need to access hits collections from calorimeters, add two integer data members to keep track of the hits collection ids.

Extra question: what are the data members of the base class G4Run?

Solution

Exercise 1 Step 2

Accumualte physics quantities

Modify file Run.cc in the method RecordEvent. This method is called by Geant4 kernel at the end of each event passing current event pointer. In this method retrieve the hits collections of calorimeters, loop on all hits and calculate physics quantities. In the constructor of Run class initialize the class data members to an initial value (0 for energy and shape and -1 for ids).

Hint 1: Note that the initial value of -1 for hits id allows you to be efficient in searching the hits by string: if id==-1 you need to search the collections, if not you already did this opeartion and you can skip it.

Solution

Exercise 1 Step 3

Implement reduction for multi-threading.

This step is optional if you do not have a multi-threading enabled application. However in this case the code is so simple that it's worth to do it (and in a sequential application this code will simply be not executed).

Why you need this? Remember in a multi-threaded application each worker thread has its own instance of class G4Run. Events are distributed and you finish up with many run objects (one per worker thread). Geant4 provides a way to merge these sub-runs into a single one. This is done implementing a Merge method. Geant4 kernel works in way that the worker threads will call the Merge method of the master run object passing a pointer to the worker run object. This simple animation explains what is happening under the hood (note that Geant4 kernel will take care of synchronizing the threads):

Solution

Exercise 1 Step 4

Create an instance of user-defined run class at each new run.

Now that you have extended G4Run you need to tell Geant4 kernel to use it instead of the default one. To do so you need to modify method RunAction::GenerateRun and return an instance of Run instead of the default (this method is called by Geant4 at the beginning of each run). The method is implemented in RunAction.cc file.

Solution

Exercise 2

Calculate physics quantities and pring them on screen.

Now that Run class has been modified to include user data we can print out our simple data analysis at the end of the run. To do that we modify the method EndOfRunAction of the RunAction class (RunAction.cc file). Retrieve from the run object the information you need and calculate the average energy release in calorimeters and the shower shape.

Hint 1: Note that Geant4 will pass you an object of type G4Run (the base class). You need to make an appropriate cast to access your data.

Hint 2: The total number of events is a data member of base class G4Run. Check in online documentation how to get it.

Hint 3: The quantity have been stored in Geant4 natural units. A useful function G4BestUnit can be used to print on screen a variable with a dimension. For example:

G4double someValue = 0.001*GeV;
G4cout<< G4BestUnit( someValue , "Energy" )<<G4endl; //Will print "1 MeV"

Solution

Using g4analysis to store results

In these exercises we will use G4AnalysisManager to store in ntuples and histograms the content of hits collections. The goal of the g4analysis module is to provide light-weight support for simple storage of data. You may skip this example if you already know that you will not use g4analysis. You need an AIDA-compliant tool or ROOT to visualize histograms.
With ROOT this is shows how histograms look like:

While this is how they look like in JAS3:

Hint 1: Histograms are automatically merged from all worker threads. With a concept similar to what shown in Exercise 1 Step 3 histograms are summed at the end of the run. A single histograms file "SlacTutorial.[root|xml]" exists. For ntuples it has not so much sense the automatic merging, because what you will do for analysis is to process one file after the other. In ROOT terminology you will create a TChain ( ITupleFactory::createChained in AIDA); with CSV format the merging is as simple as: cat *_t*.csv > merged.csv.

Exercise 3 Step 1 (Optional)

Define content of the output file(s).

Create output file(s) and define their content: four histograms and one nutple.

• First histogram: number of hits in first arm drift-chambers.
• Second histogram: number of hits in second arm drift-chambers.
• Third histogram: 2D histo containing X vs Y of all hits in first arm drift chamber.
• Foruth histogram: 2D histo containing X vs Y of all hits in second arm drift chamber.
• One ntuple containing the following data:
  • Integer column: number of hits in first drift chamber
  • Integer column: number of hits in second drift chamber
  • Double column: energy in electromagnetic calorimeter
  • Double column: energy in hadronic calorimeter
  • Double column: time on first hodscope (whatever time if multiple hits)
  • Double column: time on second hodosope (whatever time if multiple hits)

The content of the output file can be defined in the constructor of RunAction since it does not change between runs.

Solution

Exercise 3 Step 2 (Optional)

Open the output file at each new run.

Defining an output file and its content, it is not enough, you need to explicitly open it when needed. The best is to open it at the beginning of a new run. In more complex setups you can change the file name at each new run (e.g. via UI commands), so you can produce one file output for each run.

Solution

Exercise 3 Step 3 (Optional)

Write out the file.

Output files must be explicitly written to disk and closed. It is a good idea to do that at the end of the run.

Solution

Exercise 3 Step 4 (Optional)

Fill histograms and ntuple with corresponding data.

At the end of each event you should retrieve informaiton from hits collection and fill the histograms and ntuple objects.

You can access filled hits at the end of each event in EventAction class.

Solution

Exercise 3 Step 5 (Optional)

Select output file format.

You can select the output file format including one of g4root.hh, g4xml.hh or g4csv.hh headers file. To keep code clean a very simple file Analysis.hh is used to define the persistency format.

Solution

Physics Lists

Assembling physics models into processes and then into physics lists it is a error-prone process. It may be necessary, for your specific use-case to assemble your own physics list. However most often you can start from an already pre-packaged physics list and modify it according to your needs. The easier way to proceed is to add components or replace them from physics lists using constructors.

In these exercises you will learn how to instantiate one of the pre-packaged physics lists via the physics list factory mechanism and how to modify an existing physics list via constructors.

For example in this exercise an additional construct to the physics list is used: G4StepLimiterPhysics that limits the step size to a user defined value in logical volumes that have a limiter attached (in our case the magnetic field volume, see DetectorConstruction.cc file).

Exercise 4.a

Use physics list factory mechanism.

Modify the tutorial.cc file and replace the explicit instance of FTFP_BERT with the use of G4PhysListFactory class.

Hint 1: Do not forget to include the correct .hh file.
Hint 2: G4PhysListFactory::ReferencePhysList() method can be used to instantiate the physics list corresponding to the value of the environment value PHYSLIST, if the variable is not defined FTFP_BERT is used. Run few times the application changing the physics list.
Hint 3: For any given physics list, you can activate a different electro-magnetic physics adding a prefix to physics list name:

• _LIV: Use Livermore low-E EM physics
• _PEN: Use Penelope low-E EM physics
• _EMV: Use "fast options", good for crystals, not for sampling calorimeters (emOption1)
• _EMX: "Experimental" photon and bremsstralhung options, also use fast options as in EMV (emOption2).
• _EMY: The most accurate set of parameters from standard physics (emOption3).
• _EMZ: Combination of low-E Penelope and Livermore with standard (emOption4)

Solution

Exercise 4.b

Replace a constructor for a specific physics via G4VModularPhysicsList interface.

Try replace the default hadronics ion physics modelling with the alternative one based on QMD.

Hint 1: Do not forget to include the G4IonQMDPhsyics header file.
Hint 2: Pay attention to the application output just after the Geant4 banner, it will print out modifications to the physics list.
Hint 3: The method G4VModularPhysicsList::ReplacePhysics is what you need.

Solution