Thank you for your advice.
I don’t use Decay.This is my Primarygenerator code.
//code in Primarygenerator
PrimaryGeneratorAction::PrimaryGeneratorAction()
{
G4int n_particle = 1;
fParticleGun = new G4ParticleGun(n_particle);
// default particle kinematic
G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();
G4String particleName;
G4ParticleDefinition* particle = particleTable->FindParticle(particleName = "gamma");
fParticleGun->SetParticleDefinition(particle);
fParticleGun->SetParticlePosition(G4ThreeVector((27.5-30)*cm, 0, 0));
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
PrimaryGeneratorAction::~PrimaryGeneratorAction()
{
delete fParticleGun;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void PrimaryGeneratorAction::GeneratePrimaries(G4Event* event)
{
// 随机选择Co-60的两个伽马能量之一
G4double energy = (G4UniformRand() < 0.5) ? 1.17 * MeV : 1.33 * MeV;
fParticleGun->SetParticleEnergy(energy);
// 随机生成单位球面上的方向向量（各向同性）
G4double costheta = 2.0 * G4UniformRand() - 1.0;
G4double sintheta = std::sqrt(1.0 - costheta * costheta);
G4double phi = 2.0 * CLHEP::pi * G4UniformRand();
G4double dx = sintheta * std::cos(phi);
G4double dy = sintheta * std::sin(phi);
G4double dz = costheta;
fParticleGun->SetParticleMomentumDirection(G4ThreeVector(dx, dy, dz));
fParticleGun->GeneratePrimaryVertex(event);
}
And then in Runaction.cc
namespace B1
{
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
RunAction::RunAction()
{
// add new units for dose
//
const G4double milligray = 1.e-3 * gray;
const G4double microgray = 1.e-6 * gray;
const G4double nanogray = 1.e-9 * gray;
const G4double picogray = 1.e-12 * gray;
new G4UnitDefinition("milligray", "milliGy", "Dose", milligray);
new G4UnitDefinition("microgray", "microGy", "Dose", microgray);
new G4UnitDefinition("nanogray", "nanoGy", "Dose", nanogray);
new G4UnitDefinition("picogray", "picoGy", "Dose", picogray);
// Register accumulable to the accumulable manager
G4AccumulableManager* accumulableManager = G4AccumulableManager::Instance();
accumulableManager->Register(fEdep);
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void RunAction::BeginOfRunAction(const G4Run*)
{
// inform the runManager to save random number seed
G4RunManager::GetRunManager()->SetRandomNumberStore(false);
// reset accumulables to their initial values
G4AccumulableManager* accumulableManager = G4AccumulableManager::Instance();
accumulableManager->Reset();
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void RunAction::EndOfRunAction(const G4Run* run)
{
G4int nofEvents = run->GetNumberOfEvent();
if (nofEvents == 0) return;
// Merge accumulables
G4AccumulableManager* accumulableManager = G4AccumulableManager::Instance();
accumulableManager->Merge();
// Compute dose = total energy deposit in a run and its variance
//
G4double edep = fEdep.GetValue();
const auto detConstruction = static_cast<const DetectorConstruction*>(
G4RunManager::GetRunManager()->GetUserDetectorConstruction());
G4double mass = detConstruction->GetScoringVolume()->GetMass();
G4double dose = edep / mass;
// Print
//
if (IsMaster()) {
// 计算剂量
G4double c=1132*3.7*std::pow(10,10)*3600/nofEvents*dose*2; //Co60 product two gamma in one radiodecay
G4cout << "--------------------End of Global Run-----------------------"<< G4endl;
G4cout << " The run consists of " << nofEvents << " primaries." << G4endl;
G4cout << " Energy deposit in scoring volume: " << G4BestUnit(edep, "Energy") << G4endl;
G4cout << " Dose in scoring volume: " << G4BestUnit(c, "Dose") << G4endl;
G4cout << "------------------------------------------------------------" << G4endl;
}
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void RunAction::AddEdep(G4double edep)
{
fEdep += edep;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
} // namespace B1
of course I should provide my detector.
G4VPhysicalVolume* DetectorConstruction::Construct()
{
// Get nist material manager
G4NistManager* nist = G4NistManager::Instance();
// Option to switch on/off checking of volumes overlaps
//
G4bool checkOverlaps = true;
//
// World
//
G4double world_sizeXY = 2.*m, world_sizeZ =20* cm;
G4Material* world_mat = nist->FindOrBuildMaterial("G4_AIR");
auto solidWorld =
new G4Box("World", // its name
0.5 * world_sizeXY, 0.5 * world_sizeZ, 0.5 * world_sizeZ); // its size
auto logicWorld = new G4LogicalVolume(solidWorld, // its solid
world_mat, // its material
"World"); // its name
auto physWorld = new G4PVPlacement(nullptr, // no rotation
G4ThreeVector(), // at (0,0,0)
logicWorld, // its logical volume
"World", // its name
nullptr, // its mother volume
false, // no boolean operation
0, // copy number
checkOverlaps); // overlaps checking
//
// 探测小块
//
G4double size1_r = 10 * cm;
G4double size1_h = 10 * cm;
G4ThreeVector pos_tance = G4ThreeVector(27.5 * cm, 0 * cm, 0 * cm);
auto solidtan = new G4Tubs("tance", // its name
0 * cm, 0.5 * size1_r, 0.5 * size1_h, 0, 360 * deg); // its size
auto logictan = new G4LogicalVolume(solidtan, // its solid
world_mat, // its material
"tance"); // its name
G4RotationMatrix* rotation = new G4RotationMatrix();
rotation->rotateY(90 * deg); // 绕 y 轴旋转 90 度
new G4PVPlacement(rotation, // 应用旋转
pos_tance, // 放置位置
logictan, // 逻辑体
"tance", // 名称
logicWorld, // 母体逻辑体
false, // 不进行布尔操作
1, // 副本编号
checkOverlaps); // 检查重叠
// Set Shape2 as scoring volume
//
fScoringVolume = logictan;
//
// always return the physical World
//
return physWorld;
}
Other code is same to exampleB1.
I don’t kown why the result differ so much from thoretical value.