Skip to content

Muon DIS with Pythia6 run through Geant4#280

Open
siilieva wants to merge 8 commits intomasterfrom
muonDIS_with_pythia6_run_through_G4
Open

Muon DIS with Pythia6 run through Geant4#280
siilieva wants to merge 8 commits intomasterfrom
muonDIS_with_pythia6_run_through_G4

Conversation

@siilieva
Copy link

@siilieva siilieva commented Aug 25, 2025

Create a FairTask that injects a custom process into Geant4 (since we run the MC engine through a FairRunSim instance)
In this case the process is muonDIS and is being generated using Pythia6.

The muonDIS procedure is technically a C++ copy of the muon dis generation code from muonDIS.py

Also, the famous MeanMaterianBudget is made public and is widely used in the code. A tiny update was made there allowing the geo navigator to step forward shall a particle be moving along a volume border.

Tuning was done to make sure the proposed simulation chain works correctly and efficiently! See related presentations: 2July2025 and 30July2025 and 10September2025.

On top of the previously available option to set z range for the DIS vertex, the new chain allows to fully simulate events only if DIS vertex is in a specific x-y range AND/OR a volume. This saves a lot of computing resources.

Revision of indices of TPythia6 arrays was done to correct a hiccup with the cross section, but that doesn't affect previous MC chain, just made obtaining this more difficult.

Creating draft to show result on SoftPhys meeting on 27.08.2025 first, then will go public.

ps. candidate for cherry-picking in SHiP
ps. an internal/technical note is to be prepared

@siilieva siilieva force-pushed the muonDIS_with_pythia6_run_through_G4 branch from 031eb67 to fe64122 Compare August 25, 2025 17:13
siilieva added 8 commits September 12, 2025 17:57
It adds to the Geant4 process manager a custom process
that is using Pythia6 to generate muonDIS.
to reuse the function globally in sndsw.
Using standalone CMake and don't rely on the FairRoot CMake macros
as we usually do.
Better control of severity level and make the log print in colour.
In rare cases, a track is moving right along the border of two volumes.
This makes the geo navigation stuck. To unblock it, one now moves along
the particle trajectory with a step of 1e-6cm, smaller than any position
measurement precision in SND. This way, one allows the navigator to move
away from any border at a reasonable computing cost.
When going from original p6 code(Fortran) to ROOT's TPyhia6 some arrays, like the XSEC one,
have their indices swapped, also counting starts from 0 in TPythia6 case. With this fix one can
read the P6 xsec directly from the code, no need to look into the produced txt file.
Another part of the fix is to allow using P6 or P8 when creating the DIS xsec file using existing
CL options.
This is essential when running muonDIS with activated emulsions.
@siilieva siilieva force-pushed the muonDIS_with_pythia6_run_through_G4 branch from 9244683 to 590ce52 Compare September 12, 2025 16:09
@siilieva siilieva marked this pull request as ready for review September 15, 2025 07:16
Copy link

@olantwin olantwin left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

A first batch of comments. Impressive work!

return kSUCCESS;
}

ClassImp(MuonDISProcessInjector)

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This isn't needed.

muMinus = G4MuonMinus::Definition();
muPlus = G4MuonPlus::Definition();

DISProcess = new G4MuonDISProcess();

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Does the process-manager free this? Since the process is added to two different ones, probably not?

Comment on lines +20 to +21
char *fNucleon; //! nucleon type
char *fVolumeName; //! geometry volume name to place the DIS

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'd prefer using std::string here (and if necessary accessing the C-string when needed).

Comment on lines +26 to +27
const G4ParticleDefinition *muPlus;
const G4ParticleDefinition *muMinus;

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why are these const? They are updated in Init, I'm surprised this compiles.

parser.add_argument('--y_range', help='y position range placing muonDIS interaction as two values: start and end', type=float, nargs=2, default=[-1e10,1e10])
parser.add_argument('--z_range', help='z position range placing muonDIS interaction as two values: start and end', type=float, nargs=2, default=[-1e10,1e10])
parser.add_argument('--volume', help='Name of geometry volume where to simulate the DIS; e.g. volTarget, volMuUpstreamDet, volMuDownstreamDet, volVetoPlane, Wall', default="")
parser.add_argument('--dis_xsec_file', help='Name of DIS cross section file', default="/eos/experiment/sndlhc/MonteCarlo/Pythia6/MuonDIS/muDIScrossSec.root")

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Please use the full path here. This will only work on lxplus or on systems where /eos/ is mounted via FUSE.

if ((fZStart - pos.z()) > 1e-3) {
// do nothing
} else {
string geoVolumePath =gGeoManager->GetPath();

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Can we please auto-format all the new files?

fYEnd = y_range->at(1) * cm;
fZStart = z_range->at(0) * cm;
fZEnd = z_range->at(1) * cm;
LOG(INFO) << "G4MuonDISProcess: Setting xyz ranges[cm] for muon DIS\nx: "

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You mix iostream and FairLogger output. I'd suggest using one or the other.

The include of FairLogger is hidden via the include of MuonDISGenerator?

//myPythia->Pylist(1); // list table of properties for all particles
myPythia->SetMSTU(11, 6); // file number to which all program output is directed

delete myPythia;

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Reusing Pythia6 instances is impossible?

}

std::tuple<double, double, double>
G4MuonDISProcess::rotate(double ctheta, double stheta, double cphi, double sphi, double px, double py, double pz)

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is this really something we need to reimplement?

if options.debug: checking4overlaps = True
if options.debug:
checking4overlaps = True
ROOT.fair.Logger.SetConsoleSeverity("debug")

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What's the difference between this and SetLogScreenLevel? I'd use the same method for both cases. Ideally right next to each other, so that it's easy to understand the logging behaviour. Maybe where you set the log to be colourful.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants