Message: PhysicsList for RadioactiveSource with GPS Not Logged In (login)
 Next-in-Thread Next-in-Thread
 Next-in-Forum Next-in-Forum

Question PhysicsList for RadioactiveSource with GPS 

Forum: Physics List
Date: 06 Aug, 2012
From: <esther>

Dear all,

I have to simulate the decay of a cylindrical 32P-ATP radioactive source whose betas generate Cherenkov radiation, so I switched from G4ParticleGun to G4GeneralParticleSource (my reference example is N06). As a first trial, I tried to simulate a positron in matter, and verified that the results coincide with those using G4ParticleGun. The problem comes when I substitute the e+ with 32P (to simplify, I'm still assuming a point source), which I did the way you can see from my main.cc:

//CherenkovMain.cc--------------------------------------------------------------

  runManager->Initialize();

  G4UImanager* uim = G4UImanager::GetUIpointer();
  uim->ApplyCommand("/run/initialize");
  uim->ApplyCommand("/tracking/verbose 2"); 
  uim->ApplyCommand("/control/verbose 1");

  /*uim->ApplyCommand("/gps/particle e+");*/
  uim->ApplyCommand("/gps/particle ion");
  uim->ApplyCommand("/gps/ion 15 32 0 0");
  uim->ApplyCommand("/gps/position 0. 0. 0. cm");
  uim->ApplyCommand("/gps/ang/type iso");
  uim->ApplyCommand("/gps/energy 400. keV"); //<--- this is just to try, what  I'll use is the spectrum in next lines

  /*uim->ApplyCommand("/gps/ene/type Arb");
  uim->ApplyCommand("/gps/hist/type arb");

  uim->ApplyCommand("/gps/hist/point 0.0000 3.03E-02");
  uim->ApplyCommand("/gps/hist/point 0.0855 4.38E-02");
  uim->ApplyCommand("/gps/hist/point 0.1710 5.50E-02");	
  [...etc etc]
  uim->ApplyCommand("/gps/hist/inter Lin");*/

I get no errors/warnings during compilation, my simulation is running fine, but this is my SteppingVerbose output:

*************************************************************
 Geant4 version Name: geant4-09-04-patch-02    (24-June-2011)
                      Copyright : Geant4 Collaboration
                      Reference : NIM A 506 (2003), 250-303
                            WWW : http://cern.ch/geant4
*************************************************************

G4Cerenkov::G4Cerenkov constructor
NOTE: this is now a G4VProcess!
Required change in UserPhysicsList: 
change: pmanager->AddContinuousProcess(theCerenkovProcess);
to:     pmanager->AddProcess(theCerenkovProcess);
        pmanager->SetProcessOrdering(theCerenkovProcess,idxPostStep);
 AddDiscreteProcess to OpticalPhoton 
/gps/particle ion
/gps/ion 15 32 0 0
/gps/position 0. 0. 0. cm
/gps/ang/type iso
/gps/energy 400. keV

conv:   for  gamma    SubType= 14
      Lambda tables from 1.022 MeV to 10 TeV in 77 bins, spline: 1
      ===== EM models for the G4Region  DefaultRegionForTheWorld ======
        BetheHeitler :   Emin=        0 eV       Emax=         10 TeV

[...etc. etc...]

### Run 0 start.

*********************************************************************************************************
* G4Track Information:   Particle = P32[0.0],   Track ID = 1,   Parent ID = 0
*********************************************************************************************************

Step#      X         Y         Z        KineE    dEStep   StepLeng  TrakLeng    Volume     Process
    0      0 fm      0 fm      0 fm    400 keV     0 eV      0 fm      0 fm      Phantom    initStep
    1     -6 mm   2.97 mm  -1.17 cm    400 keV     0 eV   1.35 cm   1.35 cm      Phantom
    2  -5.13 m    2.54 m     -10 m     400 keV     0 eV   11.5 m    11.5 m    OutOfWorld
Number of optical photons produced in this event : 0

RunAction -> Timer: number of event = 1 User=0.01s Real=0.01s Sys=0s

CHERENKOV INFO: world_dim: 10, 10, 10 m phantom_dim: 6, 100, 100 mm phantom_pos: 0, 0, 0 cm detector_pos: 5, 0, 0 cm CCDpixel_dim: 0, 13.5, 13.5 um num_pixel: 100, 100 num_events: 1

PrimaryGeneratorAction: time: 0 ns number of particles: 1 position: 0,0,0 cm direction: -0.0445045,0.0220211,-0.086801 cm energy: 400 keV

WARNING - Attempt to delete the physical volume store while geometry closed ! WARNING - Attempt to delete the logical volume store while geometry closed ! WARNING - Attempt to delete the solid store while geometry closed ! WARNING - Attempt to delete the region store while geometry closed !

So nothing is happening: no decay, no electrons, no Cherenkov photons.

I'm pretty sure I'm missing something in my PhysicsList, but I couldn't figure out what that is from older posts and geant4 examples. I've looked at /extended/eventgenerator/exgps, /extended/radioactivedecay and /extended/electromagnetic, but they all do something different, so I don't understand what is mandatory and what isn't, and what does which thing as far as it concerns radioactive decay. Geant4 and GPS manuals haven't been useful either, so I include my PhysicsList (for completeness, I am also adding my PrimaryGeneratorAction):

//CherenkovPhysicsList.cc------------------------------------------------------- // $Id: ExN06PhysicsList.cc,v 1.19 2010-10-23 19:14:03 gum Exp $

#include "globals.hh"
#include "CherenkovPhysicsList.hh"
#include "G4ParticleDefinition.hh"
#include "G4ParticleTypes.hh"
#include "G4ParticleTable.hh"
#include "G4ProcessManager.hh"
#include "G4RadioactiveDecay.hh"
#include "G4Cerenkov.hh"
#include "G4OpAbsorption.hh"
#include "G4OpRayleigh.hh"
#include "G4OpBoundaryProcess.hh"
//#include "G4IonConstructor.hh"

CherenkovPhysicsList::CherenkovPhysicsList() :  G4VUserPhysicsList(){
  theRadioactiveDecay          = NULL;
  theCerenkovProcess           = NULL;
  theAbsorptionProcess         = NULL;
  theRayleighScatteringProcess = NULL;
  theBoundaryProcess           = NULL;
}

CherenkovPhysicsList::~CherenkovPhysicsList() {}

void CherenkovPhysicsList::ConstructParticle(){
  ConstructBosons();
  ConstructLeptons();
  ConstructBaryons();

  G4GenericIon::GenericIonDefinition();
  /*G4IonConstructor iConstructor;			//<-----Do I need it???
  iConstructor.ConstructParticle();  */
}

void CherenkovPhysicsList::ConstructBosons(){
  G4Geantino::GeantinoDefinition();
  G4ChargedGeantino::ChargedGeantinoDefinition();
  G4Gamma::GammaDefinition();
  G4OpticalPhoton::OpticalPhotonDefinition();
}

void CherenkovPhysicsList::ConstructLeptons(){
  G4Electron::ElectronDefinition();
  G4Positron::PositronDefinition();
  G4MuonPlus::MuonPlusDefinition();
  G4MuonMinus::MuonMinusDefinition();
  G4NeutrinoE::NeutrinoEDefinition();
  G4AntiNeutrinoE::AntiNeutrinoEDefinition();
  G4NeutrinoMu::NeutrinoMuDefinition();
  G4AntiNeutrinoMu::AntiNeutrinoMuDefinition();
}

void CherenkovPhysicsList::ConstructBaryons(){
  G4Proton::ProtonDefinition();
  G4AntiProton::AntiProtonDefinition();
  G4Neutron::NeutronDefinition();
  G4AntiNeutron::AntiNeutronDefinition();
} 

void CherenkovPhysicsList::ConstructProcess(){
  AddTransportation();
  ConstructGeneral();
  ConstructEM();
  ConstructOp();
}

#include "G4Decay.hh"

void CherenkovPhysicsList::ConstructGeneral(){
  // Add Decay Process
  G4Decay* theDecayProcess = new G4Decay();
  theParticleIterator->reset();
  while( (*theParticleIterator)() ){
    G4ParticleDefinition* particle = theParticleIterator->value();
    G4ProcessManager* pmanager = particle->GetProcessManager();
    if (theDecayProcess->IsApplicable(*particle)) {
      pmanager ->AddProcess(theDecayProcess);
      // set ordering for PostStepDoIt and AtRestDoIt
      pmanager ->SetProcessOrdering(theDecayProcess, idxPostStep);
      pmanager ->SetProcessOrdering(theDecayProcess, idxAtRest);
    }
  }
}

#include "G4ComptonScattering.hh"
#include "G4GammaConversion.hh"
#include "G4PhotoElectricEffect.hh"

#include "G4eMultipleScattering.hh"
#include "G4MuMultipleScattering.hh"
#include "G4hMultipleScattering.hh"

#include "G4eIonisation.hh"
#include "G4eBremsstrahlung.hh"
#include "G4eplusAnnihilation.hh"

#include "G4MuIonisation.hh"
#include "G4MuBremsstrahlung.hh"
#include "G4MuPairProduction.hh"

#include "G4hIonisation.hh"

void CherenkovPhysicsList::ConstructEM(){

// Declare radioactive decay to the GenericIon in the IonTable.

  /*const G4IonTable *theIonTable = G4ParticleTable::GetParticleTable()->GetIonTable(); //<-----Do I need it???

  G4RadioactiveDecay*  theRadioactiveDecay = new G4RadioactiveDecay();
  for (G4int i=0; i<theIonTable->Entries(); i++)
    {
      G4String particleName = theIonTable->GetParticle(i)->GetParticleName();
      if (particleName == "GenericIon")
	{
	  G4ProcessManager* pmanager = theIonTable->GetParticle(i)->GetProcessManager();
	  pmanager->SetVerboseLevel(0);
	  pmanager ->AddProcess(theRadioactiveDecay);
	  pmanager ->SetProcessOrdering(theRadioactiveDecay, idxPostStep);
	  pmanager ->SetProcessOrdering(theRadioactiveDecay, idxAtRest);
	}
    }*/

  theRadioactiveDecay = new G4RadioactiveDecay();
  //theradioactiveDecay->SetHLThreshold(-1.*s);	      //<-----Do I need it???
  //theradioactiveDecay->SetICM(true);
  //theradioactiveDecay->SetARM(false);

  theParticleIterator->reset();
  while( (*theParticleIterator)() ){
    G4ParticleDefinition* particle = theParticleIterator->value();
    G4ProcessManager* pmanager = particle->GetProcessManager();
    G4String particleName = particle->GetParticleName();

    if (particleName == "gamma") {
      pmanager->AddDiscreteProcess(new G4GammaConversion());
      pmanager->AddDiscreteProcess(new G4ComptonScattering());
      pmanager->AddDiscreteProcess(new G4PhotoElectricEffect());

    } else if (particleName == "e-") {      
      //G4LowEnergyIonisation* lowEIon = new G4LowEnergyIonisation();
      //G4LowEnergyBremsstrahlung* lowEBrem = new G4LowEnergyBremsstrahlung();  
      pmanager->AddProcess(new G4eMultipleScattering(),-1, 1, 1);
      pmanager->AddProcess(new G4eIonisation(),       -1, 2, 2);
      pmanager->AddProcess(new G4eBremsstrahlung(),   -1, 3, 3);
      //pmanager->AddProcess(new G4eBremsstrahlung(),   -1, -1, 3);

    } else if (particleName == "e+") {

      pmanager->AddProcess(new G4eMultipleScattering(),-1, 1, 1);
      pmanager->AddProcess(new G4eIonisation(),       -1, 2, 2);
      pmanager->AddProcess(new G4eBremsstrahlung(),   -1, 3, 3);
      //pmanager->AddProcess(new G4eBremsstrahlung(),   -1, -1, 3);
      pmanager->AddProcess(new G4eplusAnnihilation(),  0,-1, 4);

    } else if( particleName == "mu+" || particleName == "mu-"    ) {
     pmanager->AddProcess(new G4MuMultipleScattering(),-1, 1, 1);
     pmanager->AddProcess(new G4MuIonisation(),      -1, 2, 2);
     pmanager->AddProcess(new G4MuBremsstrahlung(),  -1, 3, 3);
     pmanager->AddProcess(new G4MuPairProduction(),  -1, 4, 4);

    } else if (particleName == "GenericIon") {
      //pmanager ->AddProcess(theRadioactiveDecay, 0, -1, 1);
      pmanager ->AddProcess(theRadioactiveDecay);
      pmanager ->SetProcessOrdering(theRadioactiveDecay, idxPostStep);
      pmanager ->SetProcessOrdering(theRadioactiveDecay, idxAtRest);
    } else {
      if ((particle->GetPDGCharge() != 0.0) &&
          (particle->GetParticleName() != "chargedgeantino")) {
       // all others charged particles except geantino
       pmanager->AddProcess(new G4hMultipleScattering(),-1,1,1);
       pmanager->AddProcess(new G4hIonisation(),       -1,2,2);
     }
    }
  }
}

void CherenkovPhysicsList::ConstructOp(){
  theCerenkovProcess           = new G4Cerenkov("Cerenkov");
  theAbsorptionProcess         = new G4OpAbsorption();
  theRayleighScatteringProcess = new G4OpRayleigh();
  theBoundaryProcess           = new G4OpBoundaryProcess();

  theCerenkovProcess->SetMaxNumPhotonsPerStep(20);
  theCerenkovProcess->SetMaxBetaChangePerStep(10.0);
  theCerenkovProcess->SetTrackSecondariesFirst(true);

  /*G4OpticalSurfaceModel model = unified; //by not defining any model, 
					    GLISUR and polished surface finish are automatically assumed*/

  theParticleIterator->reset();
  while( (*theParticleIterator)() ){
    G4ParticleDefinition* particle = theParticleIterator->value();
    G4ProcessManager* pmanager = particle->GetProcessManager();
    G4String particleName = particle->GetParticleName();
    if (theCerenkovProcess->IsApplicable(*particle)) {
      pmanager->AddProcess(theCerenkovProcess);
      pmanager->SetProcessOrdering(theCerenkovProcess,idxPostStep);  
    }
    if (particleName == "opticalphoton") {
      G4cout << " AddDiscreteProcess to OpticalPhoton " << G4endl;
      pmanager->AddDiscreteProcess(theAbsorptionProcess);
      pmanager->AddDiscreteProcess(theRayleighScatteringProcess);
      //pmanager->AddDiscreteProcess(theMieHGScatteringProcess);
      pmanager->AddDiscreteProcess(theBoundaryProcess);
    }
  }
}

void CherenkovPhysicsList::SetCuts(){

  SetCutsWithDefault();
}

//CherenkovPrimaryGeneratorAction.cc-------------------------------------------- // $Id: ExN06PrimaryGeneratorAction.cc,v 1.6 2006-06-29 17:54:27 gunter Exp $

#include "CherenkovPrimaryGeneratorAction.hh"
#include "Randomize.hh"
#include "G4Event.hh"
#include "G4GeneralParticleSource.hh"
#include "G4ParticleTable.hh"
#include "G4ParticleDefinition.hh"

CherenkovPrimaryGeneratorAction::CherenkovPrimaryGeneratorAction(){
  particleGun = new G4GeneralParticleSource();
}

CherenkovPrimaryGeneratorAction::~CherenkovPrimaryGeneratorAction(){
  delete particleGun;
}

void CherenkovPrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent){
  particleGun->GeneratePrimaryVertex(anEvent);
}

Any help/comment/suggestion will be greatly appreciated. Thanks in advance,

Esther

Inline Depth:
 1 1
 All All
Outline Depth:
 1 1
 2 2
 All All
Add message: (add)

1 Idea: Re: PhysicsList for RadioactiveSource with GPS   (Vladimir Ivanchenko - 27 Aug, 2012)
(_ None: Re: PhysicsList for RadioactiveSource with GPS   (esther - 02 Sep, 2012)
 Add Message Add Message
to: "PhysicsList for RadioactiveSource with GPS"

 Subscribe Subscribe

This site runs SLAC HyperNews version 1.11-slac-98, derived from the original HyperNews


[ Geant 4 Home | Geant 4 HyperNews | Search | Request New Forum | Feedback ]