Merge pull request #21 from mscott201/nuPRISM/bugFix/muonDecayBug

Fixing muon decay bug in current WCSim
nuPRISM · Oct 31, 2016 · dfc0014 · dfc0014
2 parents 4a53fb5 + 7f17e2e
commit dfc0014
Show file tree

Hide file tree

Showing 2 changed files with 15 additions and 11 deletions.
diff --git a/src/WCSimPhysicsList.cc b/src/WCSimPhysicsList.cc
@@ -323,11 +323,15 @@ void WCSimPhysicsList::ConstructGeneral()
   while( (*theParticleIterator)() ){
     G4ParticleDefinition* particle = theParticleIterator->value();
     G4ProcessManager* pmanager = particle->GetProcessManager();
+    G4String particleName = particle->GetParticleName();
     if (theDecayProcess->IsApplicable(*particle)) { 
       pmanager ->AddProcess(theDecayProcess);
       // set ordering for PostStepDoIt and AtRestDoIt
       pmanager ->SetProcessOrdering(theDecayProcess, idxPostStep);
-      pmanager ->SetProcessOrdering(theDecayProcess, idxAtRest);
+      // Mu minus decay at rest is handled by MuonMinusCaptureAtRest process
+      if(particleName != "mu-"){
+          pmanager ->SetProcessOrdering(theDecayProcess, idxAtRest);
+      }
     }
   }
 }

diff --git a/src/WCSimTrackingAction.cc b/src/WCSimTrackingAction.cc
@@ -10,7 +10,7 @@
 WCSimTrackingAction::WCSimTrackingAction()
 {
   ProcessList.insert("Decay") ;
-  //ProcessList.insert("MuonMinusCaptureAtRest") ;
+//  ProcessList.insert("muMinusCaptureAtRest");
 //   ProcessList.insert("conv");
   ParticleList.insert(111); // pi0
   ParticleList.insert(211); // pion+
@@ -44,7 +44,6 @@ void WCSimTrackingAction::PostUserTrackingAction(const G4Track* aTrack)
   // added by M Fechner
   const G4VProcess* creatorProcess = aTrack->GetCreatorProcess();
   //  if ( creatorProcess )
-  //    G4cout << "process name " << creatorProcess->GetProcessName() << G4endl;
 
 
   WCSimTrackInformation* anInfo;
@@ -57,20 +56,21 @@ void WCSimTrackingAction::PostUserTrackingAction(const G4Track* aTrack)
   // is the particle in the set ?
   // is it a gamma 
   // due to lazy evaluation of the 'or' in C++ the order is important
-  if( aTrack->GetParentID()==0 || 
-      ((creatorProcess!=0) && ProcessList.count(creatorProcess->GetProcessName()) ) || 
-      (ParticleList.count(aTrack->GetDefinition()->GetPDGEncoding()) )
-      || (aTrack->GetDefinition()->GetPDGEncoding()==22 && aTrack->GetTotalEnergy() > 50.0*MeV)
+  if( aTrack->GetParentID()==0 
+      || ((creatorProcess!=0) && ProcessList.count(creatorProcess->GetProcessName()))
+      || (ParticleList.count(aTrack->GetDefinition()->GetPDGEncoding()))
+      || (aTrack->GetDefinition()->GetPDGEncoding()==22 && aTrack->GetTotalEnergy() > 1.0*MeV)
+      || (creatorProcess->GetProcessName() == "muMinusCaptureAtRest" && aTrack->GetTotalEnergy() > 1.0*MeV)
       )
   {
     // if so the track is worth saving
     anInfo->WillBeSaved(true);
 
-    //      G4cout << "track # " << aTrack->GetTrackID() << " is worth saving\n";
-    //      G4cout << "It is a " <<aTrack->GetDefinition()->GetParticleName() << G4endl;
+          //G4cout << "track # " << aTrack->GetTrackID() << " is worth saving\n";
+          //G4cout << "It is a " <<aTrack->GetDefinition()->GetParticleName() << G4endl;
   }
-  else
-    anInfo->WillBeSaved(false);
+  else 
+      anInfo->WillBeSaved(false);
 
   if (aTrack->GetDefinition()->GetPDGEncoding()==111)
     pi0List.insert(aTrack->GetTrackID()); // list of all pi0-s