divaldo
/
Medtech
1
0
Fork 0
Code Issues 1 Pull Requests Releases Wiki Activity

Latest changes. 1B electrons

This commit is contained in:
Gitea 2018-04-12 14:43:05 +02:00
parent 50cad925d6
commit eccdc31423
10 changed files with 122 additions and 20 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

5
include/MedtechAnalysis.hh
View File

@ -24,20 +24,21 @@ private:
/* Private constructor to prevent instancing. */
TTree *tree;
TTree *electrontree;
TFile *file;
G4Mutex MedtechAnalysisMutex;
/* Static access method. */
MedtechAnalysis();
double x, y, e;
double x, y, z, e;
public:
~MedtechAnalysis();
MedtechAnalysis(const MedtechAnalysis&) = delete;
MedtechAnalysis& operator=(const MedtechAnalysis&) = delete;
void Fill(double, double, double);
void Fill(int, double, double, double, double);
void Close();
static MedtechAnalysis* getInstance();

2
include/MedtechDetectorConstruction.hh
View File

@ -27,6 +27,8 @@
#include "G4VPVParameterisation.hh"
#include "G4RotationMatrix.hh"
#include "G4PVReplica.hh"
#include "G4Orb.hh"
#include "G4Sphere.hh"
class MedtechDetectorConstruction : public G4VUserDetectorConstruction
{

1
include/MedtechSteppingAction.hh
View File

@ -13,6 +13,7 @@
#include "MedtechEventAction.hh"
#include "G4Gamma.hh"
#include "MedtechAnalysis.hh"
#include "G4Electron.hh"
class MedtechSteppingAction : public G4UserSteppingAction
{

9
include/Parameters.hh
View File

@ -29,6 +29,9 @@ private:
double cone_size;
double radius;
double globe_size;
double globe_placement;
public:
/* Static access method. */
static Parameters* getInstance();
@ -55,6 +58,12 @@ public:
void SetTubeSize(double);
double GetTubeSize();
void SetGlobeSize(double);
double GetGlobeSize();
void SetGlonePlacement(double);
double GetGlobePlacement();
};
#endif /* Parameters_hh */

14
medtech.cc
View File

@ -23,7 +23,7 @@
int main(int argc,char** argv)
{
Parameters *param = Parameters::getInstance();
bool visualization = false;
bool visualization = true;
int NoE=1;
for(int i=1; i < argc-1; i++)
@ -67,12 +67,22 @@ int main(int argc,char** argv)
{
param -> SetConeSize(atof(argv[i+1]));
}
else if(!strcmp("-g", argv[i])) //Globe size
{
param -> SetGlobeSize(atof(argv[i+1]));
}
else if(!strcmp("-gp", argv[i])) //Globe placement
{
param -> SetGlonePlacement(atof(argv[i+1]));
}
}
//parameters from command line
#ifdef G4MULTITHREADED
G4MTRunManager* runManager = new G4MTRunManager;
unsigned nthreads = sysconf(_SC_NPROCESSORS_ONLN); //Use all available cores
unsigned int nthreads = (int)sysconf(_SC_NPROCESSORS_ONLN); //Use all available cores
std::cout << "Using all cores (" << nthreads << ")\n";
runManager->SetNumberOfThreads(nthreads);
#else

18
src/MedtechAnalysis.cc
View File

@ -13,13 +13,20 @@ MedtechAnalysis::MedtechAnalysis()
Parameters *parameter = Parameters::getInstance();
snprintf(filename, 30, "data_c%dt%d.root", (int)parameter -> GetConeSize(), (int)parameter -> GetTubeSize());
tree = new TTree("tree", "tree");
electrontree = new TTree("electrontree","electrontree");
file = new TFile(filename,"RECREATE");
instance = this;
MedtechAnalysisMutex = G4MUTEX_INITIALIZER;
tree -> Branch("x", &x, "x/D");
tree -> Branch("y", &y, "y/D");
tree -> Branch("z", &z, "z/D");
tree -> Branch("e", &e, "e/D");
electrontree -> Branch("x", &x, "x/D");
electrontree -> Branch("y", &y, "y/D");
electrontree -> Branch("e", &e, "e/D");
electrontree -> Branch("z", &z, "z/D");
}
MedtechAnalysis::~MedtechAnalysis()
@ -38,15 +45,23 @@ MedtechAnalysis* MedtechAnalysis::getInstance()
return instance;
}
void MedtechAnalysis::Fill(double x1, double y1, double e1)
void MedtechAnalysis::Fill(int num, double x1, double y1, double z1, double e1)
{
G4AutoLock lock(&MedtechAnalysisMutex);
x = x1;
y = y1;
z = z1;
e = e1;
if(num == 1)
{
tree -> Fill();
}
if(num == 2)
{
electrontree -> Fill();
}
lock.unlock();
}
@ -54,6 +69,7 @@ void MedtechAnalysis::Fill(double x1, double y1, double e1)
void MedtechAnalysis::Close()
{
tree -> Write();
electrontree -> Write();
file -> Close();
}

51
src/MedtechDetectorConstruction.cc
View File

@ -27,7 +27,7 @@ G4VPhysicalVolume* MedtechDetectorConstruction::Construct()
G4Material *targetMaterial = manager -> FindOrBuildMaterial("G4_Pb");
G4Material *boxMaterial = manager -> FindOrBuildMaterial("G4_W");
G4Material *coneMaterial = manager -> FindOrBuildMaterial("G4_Fe");
G4Material *waterMaterial = manager -> FindOrBuildMaterial("G4_WATER");
G4Material *globeMaterial = manager -> FindOrBuildMaterial("G4_W");
//End of Material Section
//Physical world
@ -114,12 +114,12 @@ G4VPhysicalVolume* MedtechDetectorConstruction::Construct()
{
cone_placement = abs(cone_size);
solidCone = new G4Cons("cone", 0, radius*cm, 0, 0, cone_size*cm, 0, 2*CLHEP::pi); //name, inside R -Z, outside, R -Z, inside R Z, outside R Z, half length in z, vmi, vmi
solidCone = new G4Cons("cone", 0, radius*cm, 0, 0, cone_size*mm, 0, 2*CLHEP::pi); //name, inside R -Z, outside, R -Z, inside R Z, outside R Z, half length in z, vmi, vmi
//move_in_z_to_zero = cone_size; //put the cone to zero
cone_placement = cone_size + 2*tube_size; //use tube size + threshold
logicCone = new G4LogicalVolume(solidCone, coneMaterial, "Cone", 0, 0, 0);
physicalCone = new G4PVPlacement(0,
G4ThreeVector(0, 0, cone_placement*cm),
G4ThreeVector(0, 0, cone_placement*mm),
logicCone,
"Cone",
logicWorld,
@ -129,17 +129,18 @@ G4VPhysicalVolume* MedtechDetectorConstruction::Construct()
else if(cone_size < 0) //cone_size < 0
{
cone_placement = abs(cone_size) + 2*tube_size;
G4double abs_cone_size = abs(cone_size);
cone_placement = abs_cone_size + 2*tube_size;
G4RotationMatrix *rot = new G4RotationMatrix();
rot -> rotateX(180 * deg);
solidCone = new G4Cons("cone", 0, radius*cm, 0, 0, -cone_size*cm, 0, 2*CLHEP::pi); //name, inside R -Z, outside, R -Z, inside R Z, outside R Z, half length in z, vmi, vmi
solidTube = new G4Tubs("tube", 0, radius*cm, -cone_size*cm, 0, 2*CLHEP::pi); //name, inner R, outter R, Half length in Z, starting angle, angle of the segment in rad
solidCone = new G4Cons("cone", 0, radius*cm, 0, 0, -cone_size*mm, 0, 2*CLHEP::pi); //name, inside R -Z, outside, R -Z, inside R Z, outside R Z, half length in z, vmi, vmi
solidTube = new G4Tubs("tube", 0, radius*cm, -cone_size*mm, 0, 2*CLHEP::pi); //name, inner R, outter R, Half length in Z, starting angle, angle of the segment in rad
G4SubtractionSolid *Cone = new G4SubtractionSolid("Cone", solidTube, solidCone, rot, G4ThreeVector(0, 0, move_in_z_to_zero*cm));
logicCone = new G4LogicalVolume(Cone, coneMaterial, "Cone", 0, 0, 0);
physicalCone = new G4PVPlacement(0,
G4ThreeVector(0, 0, cone_placement*cm),
G4ThreeVector(0, 0, cone_placement*mm),
logicCone,
"Cone",
logicWorld,
@ -155,10 +156,10 @@ G4VPhysicalVolume* MedtechDetectorConstruction::Construct()
if(tube_size > 0)
{
G4Tubs *fixsolidTube = new G4Tubs("tube", 0, radius*cm, tube_size*cm, 0, 2*CLHEP::pi); //name, inner R, outter R, Half length in Z, starting angle, angle of the segment in rad
G4Tubs *fixsolidTube = new G4Tubs("tube", 0, radius*cm, tube_size*mm, 0, 2*CLHEP::pi); //name, inner R, outter R, Half length in Z, starting angle, angle of the segment in rad
G4LogicalVolume *fixlogicalTube = new G4LogicalVolume(fixsolidTube, coneMaterial, "fixTube", 0, 0, 0);
fixTube = new G4PVPlacement(0,
G4ThreeVector(0, 0, cone_placement*cm),
G4ThreeVector(0, 0, cone_placement*mm),
fixlogicalTube,
"Cone",
logicWorld,
@ -166,6 +167,38 @@ G4VPhysicalVolume* MedtechDetectorConstruction::Construct()
0);
}
//Globe
G4PVPlacement *Globe = 0;
G4double globe_size = parameter -> GetGlobeSize() / 2;
G4double globe_placement = globe_size + cone_placement + 100;
if(globe_size > 0)
{
G4Orb *solidGlobe = new G4Orb("Globe", globe_size*mm);
G4LogicalVolume *logicalGlobe = new G4LogicalVolume(solidGlobe, globeMaterial, "Globe", 0, 0, 0);
Globe = new G4PVPlacement(0,
G4ThreeVector(0, 0, 0.5*m),
logicalGlobe,
"Globe",
logicWorld,
false,
0);
}
//End of globe
//köpeny a másik köré
G4PVPlacement *sphere = 0;
G4double sphere_outter_size = 200 / 2;
G4double sphere_inner_size = 190 / 2;
G4Sphere *solidSphere = new G4Sphere("sphere", sphere_inner_size, sphere_outter_size, 0, 360*deg, 0, 360*deg); //name, inner R, outter R, starting phi angle, Delta phi angle, starting theta angle, delta theta angle
G4LogicalVolume * logicalSphere = new G4LogicalVolume(solidSphere, air, "sphere", 0, 0, 0);
sphere = new G4PVPlacement(0,
G4ThreeVector(0, 0, 0.5*m),
logicalSphere,
"sphere",
logicWorld,
false,
0);
return physicalWorld;
}

4
src/MedtechPrimaryGeneratorAction.cc
View File

@ -11,7 +11,7 @@ MedtechPrimaryGeneratorAction::MedtechPrimaryGeneratorAction()
{
Parameters *param = Parameters::getInstance();
G4int numberOfParticles = 1000;
G4int numberOfParticles = 1000000000;
particleGun = new G4ParticleGun(numberOfParticles);
G4ParticleTable *particleTable = G4ParticleTable::GetParticleTable();
G4ParticleDefinition *particle = particleTable -> FindParticle("e-");
@ -28,6 +28,6 @@ MedtechPrimaryGeneratorAction::~MedtechPrimaryGeneratorAction()
void MedtechPrimaryGeneratorAction::GeneratePrimaries(G4Event *anEvent)
{
particleGun -> SetParticlePosition(G4ThreeVector(0, 0, 0));
particleGun -> SetParticlePosition(G4ThreeVector(0, 0, -20*cm));
particleGun -> GeneratePrimaryVertex(anEvent);
}

14
src/MedtechSteppingAction.cc
View File

@ -66,7 +66,7 @@ void MedtechSteppingAction::UserSteppingAction(const G4Step* step)
if(postkinE < 1*MeV)
{
fTrack -> SetTrackStatus(fStopAndKill);
std::cout << "Low energy killed" <<std::endl;
//std::cout << "Low energy killed" <<std::endl;
}
else
{
@ -74,10 +74,20 @@ void MedtechSteppingAction::UserSteppingAction(const G4Step* step)
{
//fTrack -> SetTrackStatus(fStopAndKill);
MedtechAnalysis *man = MedtechAnalysis::getInstance();
man -> Fill(postX, postY, postkinE);
//std::cout << "Target" <<std::endl;
man -> Fill(1, postX, postY, postZ, postkinE);
}
}
}
if(particle == G4Electron::ElectronDefinition())
{
if(preName == "sphere")
{
MedtechAnalysis *man = MedtechAnalysis::getInstance();
std::cout << "Globe to sphere" <<std::endl;
man -> Fill(2, postX, postY, postZ, postkinE);
}
}
//G4cout << "X: " << postX << " Y: " << postY << " KinE: " << postkinE << " Particle: " << fTrack -> GetDefinition() << G4endl;

22
src/Parameters.cc
View File

@ -17,7 +17,7 @@ Parameters* Parameters::getInstance()
return instance;
}
Parameters::Parameters() : ParticleEnergy(6), hdegree(0), vdegree(0), box_size(10), tube_size(2), cone_size(2), radius(5)
Parameters::Parameters() : ParticleEnergy(6), hdegree(0), vdegree(0), box_size(10), tube_size(2), cone_size(2), radius(5), globe_size(1), globe_placement(0)
{
}
@ -104,3 +104,23 @@ double Parameters::GetTubeSize()
{
return tube_size;
}
void Parameters::SetGlobeSize(double s)
{
globe_size = s;
}
double Parameters::GetGlobeSize()
{
return globe_size;
}
void Parameters::SetGlonePlacement(double p)
{
globe_placement = p;
}
double Parameters::GetGlobePlacement()
{
return globe_placement;
}