Muon multiple scattering in geant4.10 p02

[thnam]

I did not see any muon scattered into the silicon detectors with Al100 target, and suspected that has something to do with the physics list.

On geant4.10 p02, I noticed that the QGSP_BERT physics builder use a new muon multiple scattering model instead of the Urban model:

export PHYSICSLIST=QGSP_BERT           # QGSP_BERT;QGSP_BERT_HP;QGSP_BERT_noRadi;QGSP_BERT_HP_noRadi

msc:   for mu-    SubType= 10
      RangeFactor= 0.2, step limit type: 0, lateralDisplacement: 1, polarAngleLimit(deg)= 180
      ===== EM models for the G4Region  DefaultRegionForTheWorld ======
        WentzelVIUni :  Emin=        0 eV    Emax=       10 TeV  Table with 77 bins Emin=    100 eV    Emax=     10 TeV

Then I test these model using a custom code with manually built physics list. The geometry is similar to the g4sim:

The direction of muons exited from the target is recorded, and its angles with X, Y, Z axes are plotted:

The plot shows that no muons scattered close to X-axis (where the silicon detectors are), hence explains the problem.

Then I tested the Coulomb scattering (discrete process) and the old Urban multiple scattering model:


So it is clear that Urban model gives a better description angular distribution and we should use this instead of the WentzelVIUni.

On a side note, when checking the printed physics list of g4sim I see that Coulomb scattering is used for a lot of particles. That could be one reason why Andy's simulation took so long. In my comparison above, a run with Coulomb scattering is about 6 times longer than a run using either Wentzel or Urban model (8 min vs 1.3 min).

[AndrewEdmonds11]

Thanks, Nam. Good spot.

So it looks like we should change our physics list to use the Urban model. Have you already made the change in the code?

[thnam]

No, I do not know how to do it with the built-in physics list. One possibility is hack the geant4 source code ...

[thnam]

The QGSP_BERT would still give a correct angular distribution though, because the Wentzel model is used in combination with the Coulomb scattering, according to this discussion: http://hypernews.slac.stanford.edu:5090/HyperNews/geant4/get/emprocess/957.html.

msc:   for mu-    SubType= 10
      RangeFactor= 0.2, step limit type: 0, lateralDisplacement: 1, polarAngleLimit(deg)= 180
      ===== EM models for the G4Region  DefaultRegionForTheWorld ======
        WentzelVIUni :  Emin=        0 eV    Emax=       10 TeV  Table with 77 bins Emin=    100 eV    Emax=     10 TeV

muIoni:   for  mu-    SubType= 2
      dE/dx and range tables from 100 eV  to 10 TeV in 77 bins
      Lambda tables from threshold to 10 TeV in 77 bins, spline: 1
      finalRange(mm)= 1, dRoverRange= 0.2, integral: 1, fluct: 1, linLossLimit= 0.01
      ===== EM models for the G4Region  DefaultRegionForTheWorld ======
            ICRU73QO :  Emin=        0 eV    Emax=      200 keV
          BetheBloch :  Emin=      200 keV   Emax=        1 GeV
        MuBetheBloch :  Emin=        1 GeV   Emax=       10 TeV

muBrems:   for  mu-    SubType= 3
      dE/dx and range tables from 100 eV  to 10 TeV in 77 bins
      Lambda tables from threshold to 10 TeV in 77 bins, spline: 1
      ===== EM models for the G4Region  DefaultRegionForTheWorld ======
              MuBrem :  Emin=        0 eV    Emax=       10 TeV

muPairProd:   for  mu-    SubType= 4
      dE/dx and range tables from 100 eV  to 10 TeV in 77 bins
      Lambda tables from threshold to 10 TeV in 77 bins, spline: 1
      Sampling table 17x1001; from 1 GeV to 10 TeV
      ===== EM models for the G4Region  DefaultRegionForTheWorld ======
          muPairProd :  Emin=        0 eV    Emax=       10 TeV

CoulombScat:   for  mu-    SubType= 1
      Lambda table from 100 eV  to 10 TeV in 43 bins, spline: 1
      180 < Theta(degree) < 180; pLimit(GeV^1)= 0.139531
      ===== EM models for the G4Region  DefaultRegionForTheWorld ======
  eCoulombScattering :  Emin=        0 eV    Emax=       10 TeV

[AndrewEdmonds11]

So are your plots just using custom physics lists with WentzelVI, Coulomb and Urban individually?

[thnam]

Yes.
There is a threshold to switch from WentzelVI to Coulomb somewhere in the QGSP_BERT.