Major changes in DetectorConstruction, minor changes in SteppingAction and added some parameters

include/MedtechAnalysis.hh

@@ -14,6 +14,7 @@
 #include "G4AutoLock.hh"
 #include "TTree.h"
 #include "TFile.h"
+#include "Parameters.hh"
 
 class MedtechAnalysis
 {

include/MedtechDetectorConstruction.hh

@@ -21,6 +21,9 @@
 #include "G4SubtractionSolid.hh"
 #include "Parameters.hh"
 #include "G4Tet.hh"
+#include "G4UnionSolid.hh"
+#include "G4Cons.hh"
+#include "G4Tubs.hh"
 
 class MedtechDetectorConstruction : public G4VUserDetectorConstruction
 {

include/Parameters.hh

@@ -25,6 +25,10 @@ private:
     
     double box_size;
     
+    double tube_size;
+    double cone_size;
+    double radius;
+    
 public:
     /* Static access method. */
     static Parameters* getInstance();
@@ -42,6 +46,15 @@ public:
     
     void SetBoxSize(double);
     double GetBoxSize();
+    
+    void SetRadius(double);
+    double GetRadius();
+    
+    void SetConeSize(double);
+    double GetConeSize();
+    
+    void SetTubeSize(double);
+    double GetTubeSize();
 };
 
 #endif /* Parameters_hh */

medtech.cc

@@ -14,7 +14,7 @@
 #include "G4UImanager.hh"
 #include "G4VisExecutive.hh"
 #include "G4UIExecutive.hh"
-#include "QGSP_BERT.hh"
+#include "QGSP_BIC.hh"
 #include "MedtechDetectorConstruction.hh"
 #include "MedtechActionInitialization.hh"
 #include "Parameters.hh"
@@ -23,35 +23,50 @@
 int main(int argc,char** argv)
 {
     Parameters *param = Parameters::getInstance();
-    
+    bool visualization = false;
-    int NoE=0;
+    int NoE=1;
     
     for(int i=1; i < argc-1; i++)
     {
-        if(!strcmp("-n", argv[i]))
+        if(!strcmp("-n", argv[i])) //number of events
         {
             NoE = atoi(argv[i+1]);
         }
         
-        else if(!strcmp("-p", argv[i]))
+        else if(!strcmp("-p", argv[i])) //particle energy
         {
             param -> SetParticleEnergy(atoi(argv[i+1]));
         }
         
-        else if(!strcmp("-h", argv[i]))
+        else if(!strcmp("-h", argv[i])) //H degree
         {
             param -> SetHDegree(atoi(argv[i+1]));
         }
         
-        else if(!strcmp("-v", argv[i]))
+        else if(!strcmp("-v", argv[i])) //V degree
         {
             param -> SetVDegree(atoi(argv[i+1]));
         }
         
-        else if(!strcmp("-s", argv[i]))
+        else if(!strcmp("-s", argv[i])) //windows size
         {
             param -> SetBoxSize(atof(argv[i+1]));
         }
+        
+        else if(!strcmp("-x", argv[i])) //visualization
+        {
+            visualization = false;
+        }
+        
+        else if(!strcmp("-t", argv[i])) //tube size
+        {
+            param -> SetTubeSize(atof(argv[i+1]));
+        }
+        
+        else if(!strcmp("-c", argv[i])) //Cone size
+        {
+            param -> SetConeSize(atof(argv[i+1]));
+        }
     }
     //parameters from command line
     
@@ -63,17 +78,16 @@ int main(int argc,char** argv)
 #else
     G4RunManager* runManager = new G4RunManager;
 #endif
-
+    runManager -> SetUserInitialization(new QGSP_BIC);
-    runManager -> SetUserInitialization(new QGSP_BERT());
+    runManager -> SetUserInitialization(new MedtechDetectorConstruction);
-    runManager -> SetUserInitialization(new MedtechDetectorConstruction());
+    runManager -> SetUserInitialization(new MedtechActionInitialization);
-    runManager -> SetUserInitialization(new MedtechActionInitialization());
     
     G4VisManager* visManager = new G4VisExecutive;
     visManager->Initialize();
     G4UImanager* UImanager = G4UImanager::GetUIpointer();
     
     
-    if (NoE!=0)
+    if (visualization)
     {
         //batch mode
         runManager->Initialize();

src/MedtechAnalysis.cc

@@ -9,8 +9,11 @@
 
 MedtechAnalysis::MedtechAnalysis()
 {
+    char filename[20];
+    Parameters *parameter = Parameters::getInstance();
+    snprintf(filename, 30, "data_c%dt%d.root", (int)parameter -> GetConeSize(), (int)parameter -> GetTubeSize());
     tree = new TTree("tree", "tree");
-    file = new TFile("data.root","RECREATE");
+    file = new TFile(filename,"RECREATE");
     instance = this;
     MedtechAnalysisMutex = G4MUTEX_INITIALIZER;
     

src/MedtechDetectorConstruction.cc

@@ -25,14 +25,13 @@ G4VPhysicalVolume* MedtechDetectorConstruction::Construct()
     G4NistManager *manager = G4NistManager::Instance();
     G4Material *air = manager -> FindOrBuildMaterial("G4_AIR");
     G4Material *targetMaterial = manager -> FindOrBuildMaterial("G4_Pb");
-    //G4Material *collMaterial = manager -> FindOrBuildMaterial("G4_W");
     G4Material *boxMaterial = manager -> FindOrBuildMaterial("G4_W");
     G4Material *tetMaterial = manager -> FindOrBuildMaterial("G4_Fe");
     //End of Material Section
     
     //Physical world
     G4int world_xyz = 1*m; // world is 1m x 1m x 2m for each
-    G4Box *solidWorld = new G4Box("World",0.5*world_xyz,0.5*world_xyz,4*m); //define a box for world
+    G4Box *solidWorld = new G4Box("World",world_xyz,world_xyz,4*m); //define a box for world
     
     G4LogicalVolume *logicWorld = new G4LogicalVolume(solidWorld, air, "World", 0, 0, 0); //define a logical volume for world
     
@@ -45,7 +44,7 @@ G4VPhysicalVolume* MedtechDetectorConstruction::Construct()
                                                      0); //no magnetic field
     //End of Physical world
     
-    //Target box
+    //Target box == screen
     G4Box *solidTarget = new G4Box("Target", 1*m, 1*m, 0.1*cm);
     
     G4LogicalVolume *logicTarget = new G4LogicalVolume(solidTarget, targetMaterial, "Collimator", 0, 0, 0);
@@ -59,9 +58,24 @@ G4VPhysicalVolume* MedtechDetectorConstruction::Construct()
                                                       0); //no particular field
     //End of Target box
     
+    //Water after screen
+    /*G4Box *solidWater = new G4Box("Water", 1*m, 1*m, 0.1*cm);
+    
+    G4LogicalVolume *logicTarget = new G4LogicalVolume(solidTarget, targetMaterial, "Collimator", 0, 0, 0);
+    
+    G4PVPlacement *physicalTarget = new G4PVPlacement(0, //no rotation
+                                                      G4ThreeVector(0, 0, 1.1*m), //at (0,0,0)
+                                                      logicTarget, //it's logical volume
+                                                      "Target", //it's name
+                                                      logicWorld, //it's mother volume
+                                                      false, //no boolean operations
+                                                      0); //no particular field
+     */
+    //End of water after screen
+    
     //Collimator
     G4double box_size = parameter -> GetBoxSize();
-    G4Box *solidB1 = new G4Box("B1", 25*cm, 25*cm, 1*cm);
+    G4Box *solidB1 = new G4Box("B1", 1*m, 1*m, 0.1*cm);
     G4Box *solidB2 = new G4Box("B2", box_size*cm, box_size*cm, 1.5*cm);
     
     G4SubtractionSolid *Box = new G4SubtractionSolid("Box", solidB1, solidB2);
@@ -76,12 +90,29 @@ G4VPhysicalVolume* MedtechDetectorConstruction::Construct()
     //End of Collimator
     
     //Cone
-    G4Tet *solidTet = new G4Tet("solidTet", G4ThreeVector(0, 0, 1*cm), G4ThreeVector(-5*cm, -5*cm, 0), G4ThreeVector(5*cm, -5*cm, 0), G4ThreeVector(0, 5*cm, 0));
+    //G4Tet *solidTet = new G4Tet("solidTet", G4ThreeVector(0, 0, parameter -> GetTetSize()*cm), G4ThreeVector(-2.5*cm, -2.5*cm, 0), G4ThreeVector(2.5*cm, -2.5*cm, 0), G4ThreeVector(0, 2.5*cm, 0));
-    G4LogicalVolume *logicTet = new G4LogicalVolume(solidTet, tetMaterial, "Tet", 0, 0, 0);
+    //G4Box *solidTetBox = new G4Box("solidTetBox", 1*cm, 1*cm, 1*cm);
-    G4PVPlacement *physicalTet = new G4PVPlacement(0,
+    //G4UnionSolid *unionTet = new G4UnionSolid("Cone+box", solidTet, solidTetBox);
-                                                   G4ThreeVector(0,0,1*cm),
+    
-                                                   logicTet,
+    //G4Paraboloid *solidCone = new G4Paraboloid("Cone", 2.5*cm, 0*cm, 2.5*cm);
-                                                   "Tet",
+    G4double radius = parameter -> GetRadius();
+    G4double tube_size = parameter -> GetTubeSize();
+    G4double cone_size = parameter -> GetConeSize();
+    
+    G4Cons *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
+    G4Tubs *solidTube = 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
+    
+    G4double move_z = tube_size + cone_size;
+    
+    G4UnionSolid *unionCone = new G4UnionSolid("cone", solidCone, solidTube, 0, G4ThreeVector(0, 0, -move_z*cm));
+    
+    G4double cone_placement = 2*tube_size + cone_size+ 1;
+    
+    G4LogicalVolume *logicCone = new G4LogicalVolume(unionCone, tetMaterial, "Cone", 0, 0, 0);
+    G4PVPlacement *physicalCone = new G4PVPlacement(0,
+                                                   G4ThreeVector(0, 0, cone_placement*cm),
+                                                   logicCone,
+                                                   "Cone",
                                                    logicWorld,
                                                    false,
                                                    0);

src/MedtechSteppingAction.cc

@@ -60,14 +60,28 @@ void MedtechSteppingAction::UserSteppingAction(const G4Step* step)
     if(fTrack -> GetCreatorProcess() != 0)
     {
         procN = fTrack -> GetCreatorProcess() -> GetProcessName();
-        if(postName == "Target")
+        
+        if(particle == G4Gamma::GammaDefinition())
+        {
+            if(postkinE < 1*MeV)
             {
                 fTrack -> SetTrackStatus(fStopAndKill);
+                std::cout << "Low energy killed" <<std::endl;
             }
-        //G4cout << "X: " << postX << " Y: " << postY << " KinE: " << postkinE << " Particle: " << fTrack -> GetDefinition() << G4endl;
+            else
+            {
+                if(postName == "Target")
+                {
+                    //fTrack -> SetTrackStatus(fStopAndKill);
                     MedtechAnalysis *man = MedtechAnalysis::getInstance();
                     man -> Fill(postX, postY, postkinE);
                 }
+            }
+        }
+        
+        //G4cout << "X: " << postX << " Y: " << postY << " KinE: " << postkinE << " Particle: " << fTrack -> GetDefinition() << G4endl;
+        
+    }
     
     // check if we are in scoring volume
     if (volume != fScoringVolume) return;

src/Parameters.cc

@@ -17,7 +17,7 @@ Parameters* Parameters::getInstance()
     return instance;
 }
 
-Parameters::Parameters() : ParticleEnergy(6), hdegree(0), vdegree(0), box_size(10)
+Parameters::Parameters() : ParticleEnergy(6), hdegree(0), vdegree(0), box_size(10), tube_size(2), cone_size(2), radius(2)
 {
     
 }
@@ -74,3 +74,33 @@ double Parameters::GetBoxSize()
 {
     return box_size;
 }
+
+void Parameters::SetRadius(double r)
+{
+    radius = r;
+}
+
+double Parameters::GetRadius()
+{
+    return radius;
+}
+
+void Parameters::SetConeSize(double s)
+{
+    cone_size = s;
+}
+
+double Parameters::GetConeSize()
+{
+    return cone_size;
+}
+
+void Parameters::SetTubeSize(double s)
+{
+    tube_size = s;
+}
+
+double Parameters::GetTubeSize()
+{
+    return tube_size;
+}