Introduction
The MC_LYSO project is a Geant4 simulation of an electromagnetic calorimeter prototype, consisting of a large cylindrical LYSO scintillator crystal coupled to MPPCs (SiPMs) on both faces for signal readout. The ultimate goal is to develop an auxiliary detector for the calibration of MEG-II calorimeter and investigate its potential as a prototype for future calorimeters designed for gamma particles with energies on the order of 50 MeV, offering ultra-precise time and energy resolution.
It is important to note that the output data from this application do not directly represent the final physical observables, due to the absence of the "digitization" process in the Monte Carlo. Thus, the generated data serves as input for the SiPM-response simulation: an early version of that can be found at https://github.com/lorebianco/Bartender_LYSO .
A preliminary version of a possible calibration system for the detector itself has also been implemented, aiming to assess its feasibility and performance. Please refer to the related page Calibration mode.
The MC_LYSO application will undergo continuous updates to assess various aspects.
Geometry Definition
The default setup consists of a cylindrical scintillating crystal (left image) made of LYSO - \( \text{Lu}_{1.9}\text{Y}_{0.1}\text{SiO}_{5}(\text{Ce content: 0.5mol%}) \) - with 115x2 MPPCs (Hamamatsu S13360-6025PE) attached to its Front and Back faces, in accordance with the scheme shown in the right figure.
Each MPPC is constructed as depicted in the central image, comprising a silicon layer (representing the sensitive area of the detector) placed inside an epoxy resin window, which is further enclosed in the FR4 package.
The correct placement of each SiPM is implemented in the MyDetectorConstruction::Construct() function using nested for loops and if conditions based on the panel row to determine the number of detectors to insert. Additionally, by default, two PCBs (made of FR4) and two carbon fiber endcaps are constructed, but they can be removed using interactive commands:
/MC_LYSO/my_construction/isPCB false
/MC_LYSO/my_construction/isEndcap false
Light guides can be inserted between the crystal and SiPMs using the following command:
/MC_LYSO/my_construction/isLightGuide true
The material, in the current version, can be chosen between plexiglass (default) and sapphire:
- Plexiglass:
/MC_LYSO/my_construction/MaterialOfLightGuide 1
- Sapphire:
/MC_LYSO/my_construction/MaterialOfLightGuide 2
Finally, a cylindrical coating (not in contact with the crystal) has been added solely to absorb optical photons for better graphical representation.
Physics List
The essential processes for this Monte Carlo simulation are electromagnetic and optical. The physics list has been implemented in MyPhysicsList through the RegisterPhysics() method, allowing the inclusion of pre-packaged modules from Geant4. The following modules have been loaded:
- G4EmStandardPhysics()
- G4OpticalPhysics()
The decay line for lutetium will likely be included in a future upgrade, requiring the addition of -possibly- G4DecayPhysics() and G4RadioactiveDecayPhysics().
Action Inizialization
To instantiate and register various user action classes on the G4 kernel, MyActionInitialization() has been implemented.
The user actions configured, in addition to the mandatory MyPrimaryGenerator, include MyRunAction, MyEventAction and MySteppingAction.
In sequential mode, the action classes are instantiated once by invoking the MyActionInitialization::Build() method. However, in multi-threading mode, the same method is called for each thread worker, resulting in the definition of all user action classes as thread-local instances.
The run action class is instantiated both thread-local and globally, this is why its instance is created in the MyActionInitialization::BuildForMaster() method, which is specifically invoked only in multi-threading mode.
Primary Generator
The primary generator action class employs the G4ParticleGun. In this case as well, a customizable class has been implemented with new user-defined commands. The default configuration involves generating a gamma at the origin of the world volume, directed along the positive z-axis, thus entering the cylinder perpendicularly through the Front face.
The photon's energy is sampled from a Gaussian distribution, with the mean and sigma customizable using the commands:
/MC_LYSO/my_gun/meanEnergy [value] [unit]
/MC_LYSO/my_gun/sigmaEnergy [value] [unit]
By default, the values are set to 1 MeV and 0.01 MeV. Please note that at significantly higher energies, the graphical viewer may not be able to represent all scintillation photon traces, potentially causing a crash.
Using the command
/MC_LYSO/my_gun/enableSpread true
it's possible to uniformly emit the gamma, still at the origin, with angles such that it enters the front face of the scintillator within a circle. The radius of this circle can be set with:
/MC_LYSO/my_gun/radiusSpread [value] [unit]
Runs and Events
A run consists of a set of events. Through user action classes, operations have been configured at various stages:
- at the beginning and end of each run (class MyRunAction)
- at the beginning and end of each event (class MyEventAction)
- at the end of each step (class MySteppingAction)
The last one is discussed further in Energy Deposition.
At the start of a run, a TTree with its TBranches is created. At the conclusion, the output root file is created and opened so the TTree, filled with Monte Carlo data, is then written to it.
In the method MyEventAction::BeginOfEventAction(), all event-data containers are reset. In MyEventAction::EndOfEventAction(), the data packaged in the containers are copied into a new row of the TTree along with other directly extracted information.
A detailed description of the collected data and the structure of the TTree is provided below in Data Flow and Output File.
Detector Response
The detector response is implemented through the commonly used sensitive detector + hit scheme. In this application, the silicon layers of the MPPCs are considered the detector. Therefore, the logic volume associated with them is declared as a "sensitive detector" (SD) in MyDetectorConstruction::ConstructSDandField(), and they are associated with an instance of the MySensitiveDetector() class.
This class contains a hit collection tasked with containing all hits from scintillation optical photons. A hit is represented by an instance of the MyHit class created in MySensitiveDetector::ProcessHits() when an optical photon undergoes a step in the SD and is detected. Indeed, in this version the photon detection efficiency of the MPPC (24%) is also implemented, which could potentially be moved to the SiPM-response simulation in the future.
A hit is defined as a set of the following info related to the detection of the scintillation photon:
- Detection time
- Detector position (note that the package's position is stored, not the layer's)
- Detector channel
Energy Deposition
We are also interested in extracting the energy deposited inside the crystal. However, since it is not associated with a sensitive detector, the method described above is not applicable. To extract the information, in MyDetectorConstruction, a G4LogicalVolume is defined and associated with the logical volume of the scintillator, representing the scoring volume. This volume is utilized in MySteppingAction::UserSteppingAction() to extract, after each step localized within the crystal, the deposited energy and the position of the maximum \( \frac{dE}{dx} \), identified as the midpoint between the pre and post step points where the maximum \( \frac{dE_{dep}}{step length}\) occurred.
Additionally, the entry position and time of the primary gamma into the scintillator are extracted.
Data Flow and Output File
In G4, there are various methods for extracting and transmitting data between different classes and stages of the simulation. Here is a summary of the implemented data flow and the description of the output ROOT file.
As described, in the constructor of MyRunAction, a TTree lyso is created using G4AnalysisManager, and its branches are defined. The structure is shown in the table.
| Branch | Type | Unit | Description |
|---|---|---|---|
| Event | int | - | Event identifier |
| E_gun | double | MeV | Energy of the primary gamma |
| X_gun | double | mm | Initial x position of the primary gamma |
| Y_gun | double | mm | Initial y position of the primary gamma |
| Z_gun | double | mm | Initial z position of the primary gamma |
| MomX_gun | double | - | X component of the momentum direction of the primary gamma |
| MomY_gun | double | - | Y component of the momentum direction of the primary gamma |
| MomZ_gun | double | - | Z component of the momentum direction of the primary gamma |
| ToA | double | ns | Time of arrival of the primary gamma in the scintillator |
| XoA | double | mm | X position of arrival of the primary gamma in the scintillator |
| YoA | double | mm | Y position of arrival of the primary gamma in the scintillator |
| ZoA | double | mm | Z position of arrival of the primary gamma in the scintillator |
| Edep | double | MeV | Energy deposited within the scintillator |
| MaxEdep | double | MeV/mm | Maximum dE/dx within the scintillator |
| MaxEdepPosX | double | mm | X position of the maximum dE/dx within the scintillator |
| MaxEdepPosY | double | mm | Y position of the maximum dE/dx within the scintillator |
| MaxEdepPosZ | double | mm | Z position of the maximum dE/dx within the scintillator |
| NHits_F | int | - | Number of SiPM hits on the Front face |
| NHits_B | int | - | Number of SiPM hits on the Back face |
| NHits_Tot | int | - | Total number of SiPM hits |
| T_F | vector<double> | ns | Time of detection of SiPMs on the Front face |
| X_F | vector<double> | mm | X position of SiPMs on the Front face |
| Y_F | vector<double> | mm | Y position of SiPMs on the Front face |
| Ch_F | vector<int> | - | Channel of SiPMs on the Front face |
| T_B | vector<double> | ns | Time of detection of SiPMs on the Back face |
| X_B | vector<double> | mm | X position of SiPMs on the Back face |
| Y_B | vector<double> | mm | Y position of SiPMs on the Back face |
| Ch_B | vector<int> | - | Channel of SiPMs on the Back face |
Many of these quantities are simulated at different times within an event: the strategy used then was to include, as class members of MyEventAction, all the variables and structures (referred to as data containers) necessary to store the data. The instance of MyEventAction is then provided as a parameter to the constructors of MyRunAction and MySteppingAction, enabling the sharing of information between user action classes.
At the beginning of the event, the primary gamma is generated, and its characteristics are saved by default by Geant4, thus accessible without data containers. When the particle enters the scintillator, the data related to the physics occurring in the crystal is processed through MySteppingAction::UserSteppingAction(), which extracts various quantities and writes them into the containers.
The scintillation photons are then detected by the sensitive detector, which, through MySensitiveDetector::ProcessHits(), saves the hits in the hit collection.
At the end of the event, the MyEventAction::EndOfEventAction() method is called. Here, access is made to the hit collection, and all information encapsulated in the hits is written into the vector data containers of the class.
With all event data available, the information is saved in a new row of the TTree, and at the beginning of the next event the containers are reset.
To conclude, in MyRunAction::EndOfRunAction(), the TTree lyso is written to the ROOT file output_MCID_run [runID] _t [threadNumber].root.
Summary of the Monte Carlo run
In batch mode only, the function MC_summary() is invoked at the very end of the application. It creates (or updates) a file named MC_summaries.txt where a summary of the simulation is written. The recap is extracted from the macro files used (see the next paragraph) and includes various settings. Additionally, it records the date, user name, duration, and the randomly generated Monte Carlo seed at the executable launch.
An example output is displayed:
How to run
This simulation, as often used for Geant4 applications, supports both interactive mode (with visualization) and batch mode.
Note that the provided macro files, that you can find in the macros directory of the repository, should be considered as valid templates to modify command values according to specific needs, without the necessity of adding new ones.
To run in interactive mode, simply execute the following shell command:
$ ./mc_lyso
Then, the init_vis.mac macro initializes the kernel with the essential command:
/run/initialize
Subsequently, the vis.mac macro is executed, opening the graphical viewer and configuring its settings.
For batch mode, simply provide a macro file:
$ ./mc_lyso run.mac
In the run.mac file, when in multithreading mode, execute the crucial command:
/run/numberOfThreads [value]
in order to set the number of threads to instantiate.
Before initializing the kernel, set the command to execute the construction.mac macro, defining the geometry using UI-commands (see Geometry Definition):
/control/execute construction.mac
At the end, initiate the run with the command:
/run/beamOn [nOfEvents]
An example of the output can be found in Standard mode: output example, in the related page.
Please note that G4, in order to ensure thread safety, generates an output file for each utilized thread. To merge these files, you can employ the hadd function of ROOT:
$ hadd MCID_[MCID].root output_MCID_run[runID]_t*
where [MCID] represents the initial seed of the run, serving as a serial number to identify the specific Monte Carlo. However, before closing, the application prints the name that should be given to the merged file.
