ShipSoft / FairShip

SHiP experiment framework based on FairRoot
http://ship.web.cern.ch/ship/
23 stars 107 forks source link

Geant4 only proton simulation crashing (tested for charm detector configuration) #381

Closed antonioiuliano2 closed 1 month ago

antonioiuliano2 commented 3 years ago

Dear all,

unfortunately I have to report another issue: I have been required to increase the statistics of proton simulations to study the MC fluctuations for the proton analysis note. Therefore, I have launched the same simulation used in the past, but to no avail. Simulation crashes with the following error:

Target element Pb Z= 82 A= 207 Unrecoverable error in the method ApplyYourself of protonInelastic TrackID= 1 ParentID= 0 proton Ekin(GeV)= 399.02; direction= (0.000129552,-2.57172e-05,1) Position(mm)= (3.60512,0.398464,917.556); material lead PhysicalVolume ApplyYourself failed

Fatal Exception core dump *** G4Track (0x44abdcb0) - track ID = 1, parent ID = 0 Particle type : proton - creator process : not available Kinetic energy : 399.02 GeV - Momentum direction : (0.000129552,-2.57172e-05,1) Step length : 1.8443 cm - total energy deposit : 24.5542 MeV Pre-step point : (3.60343,0.398834,899.113) - Physical volume : volPasLead (lead)

You should be able to replicate it by launching:

python $FAIRSHIP/muonShieldOptimization/run_MufluxfixedTarget.py --CharmdetSetup 1 -G -n 100 -o sim_charm/prova/

It seems to happen when the -G (G4 only) option is called. It does not always happen (i.e., not at the first event), I do not know the reason.

Best regards, Antonio

ThomasRuf commented 3 years ago

welcome to the club, see https://bugzilla-geant4.kek.jp/show_bug.cgi?id=2299

Since other people cannot reproduce it, somebody of us needs to debug the problem. By the way, why do you choose to run with G4 only?

antonioiuliano2 commented 3 years ago

Dear Thomas,

I run with Geant4 only because I need to track the primary proton too, not only the daughters. This way, I can study how the track/vertex reconstruction treats it.

Otherwise, there would be no emulsion hits from the 400 GeV protons, only from the daughters.

ThomasRuf commented 3 years ago

The same problem occurs with MuDISGenerator, where we want also to have the hits of the incoming muon. The solution which I implemented there, add a copy of the incoming muon with reversed momentum. Now, in your case, applying the same trick, will probably work not so well, since the proton could interact on the way back, which is unphyiscal. You would need to add a second particle with reversed momentum which does not interact hadronically.

antonioiuliano2 commented 3 years ago

Thank you, it is a nice suggestion, we can work on it in the future.

Concerning this issue, however, may I ask you if you got this error with the other fixed target simulation options, or only with G4only?

ThomasRuf commented 3 years ago

Don't remember for sure, but I think it was also with G4only, while doing simulations for miniShip, but I can be wrong.

antonioiuliano2 commented 3 years ago

Thank you very much for your answer.

Since as you say we need to debug the problem ourselves, as a first step I want to collect all the information we know about this issue/bug/error:

Call for FTFP Target element Fe Z= 26 A= 56 Unrecoverable error in the method ApplyYourself of pi+Inelastic TrackID= 7113 ParentID= 42 pi+ Ekin(GeV)= 23.1785; direction= (0.0284275,0.0210293,0.999375) Position(mm)= (281.882,65.5225,8319.82); material iron PhysicalVolume ApplyYourself failed

I will update you as soon as I get other info,

ThomasRuf commented 3 years ago

If I look at the code, for G4only, a proton with Pz=fMom at position xOff, yOff, start[2] is added to the stack: cpg->AddTrack(2212,0.,0.,fMom,xOff/cm,yOff/cm,start[2],-1, kTRUE, -1., 0.,1.);

If I compare this with cpg->AddTrack(id,px,py,pz, x/cm,y/cm,z/cm, im,wanttracking,e,tof,wspill,procID); I was a bit puzzled by the -1., which should be the energy of the proton. But in FairPrimaryGenerator, there is: if (e < 0) { Double_t mass = pdgBase->GetParticle(pdgid)->Mass(); e = TMath::Sqrt(mom.Mag2() + mass * mass); So, should be ok. I would then not understand why it only happens for G4only, and not when running with Pythia first. There are diffractive events in Pythia, which can also produce very high momentum protons. Also, your pion is not very high momentum.

When running muon DIS events through a full Geant4 simulation for SND (many particles with energies up 2-3TeV!), I do not observe this issue.

antonioiuliano2 commented 3 years ago

After two test runs, I can now at least confirm an error with the standard fixed target simulation (i.e. without the -G option, evtgen): python /afs/cern.ch/work/a/aiuliano/public/SHIPBuild/FairShip/muonShieldOptimization/run_MufluxfixedTarget.py --CharmdetSetup 1 -n 100 -o sim_charm/prova2/

The particle which causes the crash is another pion:

Target element Fe Z= 26 A= 56 Unrecoverable error in the method ApplyYourself of pi+Inelastic TrackID= 5 ParentID= 0 pi+ Ekin(GeV)= 43.2667; direction= (0.0323967,-0.00319487,0.99947) Position(mm)= (145.437,-20.0922,8045.83); material iron PhysicalVolume ApplyYourself failed

Fatal Exception core dump *** G4Track (0x471faad0) - track ID = 5, parent ID = 0 Particle type : pi+ - creator process : not available Kinetic energy : 43.2667 GeV - Momentum direction : (0.0323967,-0.00319487,0.99947) Step length : 4.2818 mm - total energy deposit : 4.55282 MeV Pre-step point : (145.298,-20.0783,8041.55) - Physical volume : VPassive (iron)

-------- EEEE -------- G4Exception-END --------- EEEE -------

olantwin commented 1 year ago

welcome to the club, see https://bugzilla-geant4.kek.jp/show_bug.cgi?id=2299

Since other people cannot reproduce it, somebody of us needs to debug the problem. By the way, why do you choose to run with G4 only?

This issue should now be fixed in Geant4 versions starting with 10.7.3. Once we update Geant4, we can close this issue, unless there was a second problem, @antonioiuliano2 ?

olantwin commented 1 month ago

Fixed upstream since at least stack 24.06.

If other issues persist, @antonioiuliano2 , please open a second issue.