search

alcap-org / g4sim

Simulation toolkit based on Geant4 and ROOT
http://wuchen1106.github.io/g4sim/
2 stars 2 forks source link

Incorrect time of capture products (Possible Geant4 Bug) #48

Open AndrewEdmonds11 opened 9 years ago

AndrewEdmonds11 commented 9 years ago

Ok, so I think I've found a major problem with Geant4, which I will try to explain to the best of my ability.

What The Problem Is

I was looking at the time of protons originating in the target in the thick silicon with the following draw command (NB x-axis is in ns):

tree->Draw("M_Ot>>hist", "M_particleName == \"proton\" && M_volName ==\"ESi1\" && M_ovolName==\"Target\"")

timefit

and then fitting an exponential to the results we get a lifetime of 1/0.0487 = 20.104!!

Delving into Geant4

(This is where my explanation might get confusing)

We are using all the physics that are in QGSP_BERT.icc. Within this the G4StoppingPhysics constructor is created.

In G4StoppingPhysics.cc (located in physics_lists/constructors/stopping/src), various physics processes are created, including G4MuonMinusCapture.

In the constructor for this class we have the following code:

G4MuonMinusCapture::G4MuonMinusCapture(G4HadronicInteraction* hiptr)
  : G4HadronStoppingProcess ("muMinusCaptureAtRest")
{
  std::cout << "AE: Created G4MuonMinusCapture " << std::endl;
  SetBoundDecay(new G4MuonMinusBoundDecay()); // Owned by InteractionRegistry                                                                         
  if (!hiptr) {
    hiptr = new G4CascadeInterface(); // Owned by InteractionRegistry                                                                                 
  }
  RegisterMe(hiptr);
}

Note that this class inherits from G4HadronStoppingProcess, which implements a general way with dealing with these processes.

In its AtRestDoIt, G4HadronStoppingProcesses does the following:

It is in G4MuonMinusBoundDecay that the capture rates and decay rates are calculated and capture/decay is randomly chosen. If capture is chosen, then nothing happens in this class and step 3 above is performed. Otherwise, the decay is performed and step 3 above is skipped.

It's also worth noting that there are a few other muon processes in the source code but I don't think an of it gets called (including useful stuff like DoMuCapture())

Printing Out Information

Anyway, I put some useful comments in various places and printed out what was going on.

At the start of each process, I've printed out 1/lambda (the lifetime of the muonic atom), 1/lambdac (the effective lifetime of the capture) and t (the time that should be added to the global time).

Then I print the incoming projectile (always a mu- phew) and whether it decays or captures (and the capture process it is trying to do).

Finally, I print the secondary particle and it's initial time.

Here are two stopped muons:

AE: Trying G4MuonMinusBoundDecay: 1/lambda= 872.8409729841 1/lambdac= 1431.6392269148 t= 881.76357556529
AE: mu-
AE: Decaying
AE: nuclearCapture? 0
AE: time of secondary (e-): 0
AE: time of secondary (e-): 0
AE: time of secondary (e-): 0
AE: time of secondary (gamma): 0
AE: time of secondary (gamma): 0
AE: time of secondary (gamma): 0
AE: time of secondary (gamma): 0
AE: time of secondary (gamma): 0
AE: time of secondary (gamma): 0
AE: time of secondary (e-): 881.764
AE: time of secondary (anti_nu_e): 881.764
AE: time of secondary (nu_mu): 881.764
AE: Trying G4MuonMinusBoundDecay: 1/lambda= 872.8409729841 1/lambdac= 1431.6392269148 t= 243.18241174079
AE: mu-
AE: Capturing (maybe)
AE: nuclearCapture? 1
AE: G4HadronicInteraction model = BertiniCascade
AE: time of secondary (e-): 0
AE: time of secondary (gamma): 0
AE: time of secondary (gamma): 0
AE: time of secondary (e-): 0
AE: time of secondary (gamma): 0
AE: time of secondary (gamma): 0
AE: time of secondary (gamma): 0
AE: time of secondary (nu_mu): 0
AE: time of secondary (neutron): 0
AE: time of secondary (gamma): 0
AE: time of secondary (neutron): 0
AE: time of secondary (gamma): 0
AE: time of secondary (Mg25): 0

I think the first set of secondaries (e- and gamma) are from the muon falling to the 1s state and we can see that only the products of decaying muons have time added to it.

Conclusion

I think there's a bug in geant4. I've had a quick search of the Geant4 hypernews but haven't found anything yet.

I will search more tomorrow.

benkrikler commented 9 years ago

I can't comment beyond this, but I agree that looks weird, and from your checks like a G4 bug. Do you know how the timing in the decay case is worked out? Is it a random number with an exponential PDF?

AndrewEdmonds11 commented 9 years ago

It gets the capture rate in G4MuonMinusBoundDecay::GetMuonCaptureRate() which has an array of values taken from Suzuki, Measday and Roalsvig and if it's not in that array it calculates it based on Goulard and Primakoff.

The decay rate is calculated in G4MuonMinusBoundDecay::GetMuonDecayRate() based on Mukhopadhyay.

The total rate is then capture rate + decay rate and a random time is created in G4MuonMinusBoundDecay::ApplyYourself():

G4double time = -std::log(G4UniformRand()) / lambda;

The choice between decaying and capturing is based on a uniform random value between 0 and the total rate. If this number is less than the capture rate then we do that, otherwise we do the decay.

AndrewEdmonds11 commented 9 years ago

After further looking through the code, it seems that this "BertiniCascade" is fine since it uses G4CascadeMuMinusPChannel which implements the mu + p --> n + nu_mu.

So I think the only problem is with the time of these products.

I've requested an account to the Geant4 hypernews to post the question there but am waiting for approval. In the mean time I'll try and find a kludge around it.

AndrewEdmonds11 commented 9 years ago

Found a way to kludge this :smile:

In processes/hadronic/util/include/G4HadFinalState.hh, add the following to the public area:

// So that we can get the correct times for muon capture products                                                                                   
  void ModifyTimes(G4double time) {
    for (std::vector<G4HadSecondary>::iterator i_sec = theSecs.begin(); i_sec != theSecs.end(); ++i_sec) {
      i_sec->SetTime(time);
    }
  }

and then in processes/hadronic/stopping/src/G4HadronStoppingProcess.cc add the line highlighted below at L.229:

while(!resultNuc);

    // Add times for NMC                                                                                                                              
    resultNuc->ModifyTimes(thePro.GetGlobalTime()); // this line here!!!
    edep = resultNuc->GetLocalEnergyDeposit();
    nSecondaries += resultNuc->GetNumberOfSecondaries();
    result->AddSecondaries(resultNuc);
    resultNuc->Clear();

Here's a before and after of a 50k, gen_mum_tgt run, with tree->Draw("M_t>>hist1", "M_volName==\"ESi1\" && M_ovolName==\"Target\" && M_t<5000 && M_particleName==\"proton\"", "") (NB x-axis in ns).

Before woutmodification

After wmodification

Unfortunately, I can't commit this to our repo because we always retrieve the tarball from geant during installation.

Anyway, once I get an account for the hypernews I will flag it up there.

benkrikler commented 9 years ago

This is great! This will be crucial for the main COMET and Mu2e simulations as well I imagine... We're on the most recent version of Geant4, right?

Unfortunately, I can't commit this to our repo because we always retrieve the tarball from geant during installation.

You could make a patch from the diff of your changes, add that patch to the repo and alter the install script to apply it automatically. Or just share the patch around so people can apply it themselves since there's only a few of us.

wuchen1106 commented 9 years ago

A late response: In my previous simulations I was using Geant4.9.6.p02, where the time distribution of particles after muon capture is correct. I believe this is a new bug in Geant4.10.X which should be corrected by a new patch hopefully.

AndrewEdmonds11 commented 9 years ago

Yes, I think that's right. I reported it to Geant and there is a patch here.