varaible density

parent 6e81a114
......@@ -85,10 +85,10 @@ int main(int argc,char** argv)
// Initialize visualization
// (if you don't want the image, just delete the G4VisManager objetc and everything that concerns it)
G4VisManager* visManager = new G4VisExecutive;
// G4VisManager* visManager = new G4VisExecutive;
// G4VisExecutive can take a verbosity argument - see /vis/verbose guidance.
// G4VisManager* visManager = new G4VisExecutive("Quiet");
visManager->Initialize();
// visManager->Initialize();
// Get the pointer to the User Interface manager
G4UImanager* UImanager = G4UImanager::GetUIpointer();
......@@ -113,7 +113,7 @@ int main(int argc,char** argv)
// owned and deleted by the run manager, so they should not be deleted
// in the main() program !
delete visManager; //(delete this line too if you don't want image).
// delete visManager; //(delete this line too if you don't want image).
delete runManager;
}
......
......@@ -36,6 +36,9 @@
#include "G4VisAttributes.hh"
#include <vector>
#include <cmath>
#include <iostream>
#include <fstream>
class B1SD;
class G4VPhysicalVolume;
......
......@@ -79,15 +79,17 @@ G4VPhysicalVolume* B1DetectorConstruction::Construct()
//
//caracteristicas do cilíndro -> world
G4double raio_i = 0;
G4double raio_e = 70*m;
G4double raio_e = 1*km;
G4double h_detector = 1*cm;
G4double h = 200*m;
G4double h = 1*km;
G4double theta_0 = 0.*deg;
G4double theta_f = 360.*deg;
//making world of vacuum
G4Material* world_mat = nist->FindOrBuildMaterial("G4_Galactic");
G4Material* Oxygen = nist->FindOrBuildMaterial("G4_O");
G4Material* Nitrogen = nist->FindOrBuildMaterial("G4_N");
G4Material* Argon = nist->FindOrBuildMaterial("G4_Ar");
G4String name;
G4Tubs* cilindro = new G4Tubs("world", raio_i, raio_e, h, theta_0, theta_f);
logicWorld =
......@@ -112,53 +114,48 @@ G4VPhysicalVolume* B1DetectorConstruction::Construct()
//making air materiais
//Nitrogen
G4double a_N = 14.01*g/mole;
G4String name, symbol;
G4double z;
G4Element* ele_N = new G4Element(name="Nitrogen", symbol="N", z=7., a_N);
//Oxygen
G4double a_o = 16.00*g/mole;
G4Element* ele_O = new G4Element(name="Oxygen", symbol="O", z=8., a_o);
//Argon
G4double a_Ar = 39.948*g/mole;
G4Element* ele_Ar = new G4Element(name="Argon", symbol="Ar", z=18., a_Ar);
G4int num_layers = 200;
G4double layer_delta_h = 2*h/num_layers;
std::vector<G4Material *> layers_material;
layers_material.reserve(num_layers);
std::ofstream density_file("density.csv");
for(G4int j = 0; j < num_layers; j++) {
//fazer if clause para verificar valor de delta h progressivamente.
//dependendo do h + delta_h -> parâmetros diferentes.
G4double density;
G4double atm_position = layer_delta_h*j;
if (atm_position >= 0*km && atm_position < 4*km) {
density = density_function(-86.5562, 1222.6562, 994186.38, atm_position);
G4cout << "a" << G4endl;
density = density_function(-186.5562*g/cm2, 1222.6562*g/cm2, 994186.38*cm, atm_position);
} else if (atm_position >= 4*km && atm_position < 10*km) {
density = density_function(-94.919, 1144.9069, 878153.55, atm_position);
G4cout << "b" << G4endl;
density = density_function(-94.919*g/cm2, 1144.9069*g/cm2, 878153.55*cm, atm_position);
} else if (atm_position >= 10*km && atm_position < 40*km) {
G4cout << "c" << G4endl;
density = density_function(0.61289, 1305.5948, 636143.04, atm_position);
density = density_function(0.61289*g/cm2, 1305.5948*g/cm2, 636143.04*cm, atm_position);
} else if (atm_position >= 40*km && atm_position < 100*km) {
G4cout << "d" << G4endl;
density = density_function(0.0, 540.1778, 772170.16, atm_position);
}
G4cout << density << G4endl;
G4cout << density/(g/cm2) << G4endl;
density_file << density/(g/cm2) << " " << atm_position/(km) << "\n";
G4cout << j << G4endl;
layers_material[j] = new G4Material(name="Air" + std::to_string(j), density, 3);
layers_material[j]->AddElement(ele_N, 78.1*perCent);
layers_material[j]->AddElement(ele_O, 21.0*perCent);
layers_material[j]->AddElement(ele_Ar, 0.9*perCent);
layers_material[j] = new G4Material(name="Air" + std::to_string(j), density/1000000, 3);
layers_material[j]->AddMaterial(Nitrogen, 78.1*perCent);
layers_material[j]->AddMaterial(Oxygen, 21.0*perCent);
layers_material[j]->AddMaterial(Argon, 0.9*perCent);
}
density_file.close();
std::vector<G4LogicalVolume* > air_layers;
air_layers.reserve(num_layers);
G4Tubs* cilindro_ar = new G4Tubs("air_layer", raio_i, raio_e, layer_delta_h/2 , theta_0, theta_f);
for(G4int j = 0; j < num_layers; j++) {
for(G4int j = 0; j < (num_layers - 1); j++) {
air_layers[j] = new G4LogicalVolume(cilindro_ar, //its solid
layers_material[j], //its material
"layer" + std::to_string(j + 1));
......@@ -166,6 +163,7 @@ G4VPhysicalVolume* B1DetectorConstruction::Construct()
//calculing the position in relation to the mother volume
G4double positon = h - (layer_delta_h/2) - (layer_delta_h)*j;
new G4PVPlacement(0, //no rotation
G4ThreeVector(0,0,positon), //at (0,0,0)
air_layers[j], //its logical volume
......@@ -174,21 +172,20 @@ G4VPhysicalVolume* B1DetectorConstruction::Construct()
false, //no boolean operation
j); //copynumber
}
G4double safe_distance = 0.02*m;
//Coloque aqui o número de detetores
G4int num_detector = 32;
number_detectors = num_detector;
number_detectors = num_layers - 1;
G4double total_safe = h_detector + safe_distance;
G4double delta_h = (2*h - 2*safe_distance - 4*h_detector)/(num_detector - 1);
G4double initial_pos = -h + total_safe;
G4double final_pos = h - total_safe;
size = h;
delta = delta_h;
delta = layer_delta_h;
detector = new G4Tubs("detector", raio_i, raio_e, h_detector, theta_0, theta_f);
/*
G4LogicalVolume* logicDetector = new G4LogicalVolume(detector, //its solid
world_mat, //its material
......@@ -217,23 +214,23 @@ logicDetector.reserve(num_detector);
}
*/
//Coloca os deteroes dentro do cilindro e coloca a sua posição em um vetor com push_back()
//fScoringDetector = logicDetector;
//Associando os detetores à classe sensitive detector
//Associando os camdas à classe sensitive detector
auto sdman = G4SDManager::GetSDMpointer();
for(G4int i = 0; i< num_detector; i++) {
for(G4int i = 0; i< (num_layers - 1); i++) {
G4String SDname = "SD" + std::to_string(i + 1);
auto sensitive = new B1SD(SDname);
sdman->AddNewDetector(sensitive);
logicDetector[i]->SetSensitiveDetector(sensitive);
air_layers[i]->SetSensitiveDetector(sensitive);
}
*/
//colocando cor vermelha para detectores
G4VisAttributes* worldVisAtt1 = new G4VisAttributes(G4Colour(1.0,0.0,0.0));
......@@ -247,7 +244,7 @@ logicDetector.reserve(num_detector);
logicDetector[i]->SetVisAttributes(worldVisAtt1 );
}
*/
for(G4int i = 0; i < num_layers; i++) {
for(G4int i = 0; i < (num_layers - 1); i++) {
air_layers[i]->SetVisAttributes(worldVisAtt2);
}
......
//
//
// ********************************************************************
// * License and Disclaimer *
// * *
......@@ -94,8 +94,9 @@ void B1EventAction::EndOfEventAction(const G4Event* event)
std::ofstream nu_mu_file("nu_mu.txt");
std::ofstream anti_nu_e_file("anti_nu_e.txt");
std::ofstream anti_nu_mu_file("anti_nu_mu.txt");
std::ofstream particles_file("particles.csv");
/*
std::vector<G4int> col;
std::vector<B1HitsCollection*> HitsCol;
......@@ -110,9 +111,11 @@ void B1EventAction::EndOfEventAction(const G4Event* event)
}
if(HitsCol[i]) {
G4double pos_detec = SetPosition(i + 1);
G4double calc = 2*half_height - (variation/2) - (variation)*i;
G4double pos_detec = calc/(km);
int n_hit = HitsCol[i]->entries();
G4cout << "My detector " + std::to_string(i+1) +" has " << n_hit << "hits" << G4endl;
particles_file << n_hit << " " << pos_detec << "\n";
B1Hits* hit = new B1Hits;
std::map<const G4String, int> fparticles;
for(int i1 = 0; i1 < n_hit; i1++) {
......@@ -126,7 +129,7 @@ void B1EventAction::EndOfEventAction(const G4Event* event)
fparticles.clear();
}
}
*/
gammafile.close();
mu_plusfile.close();
mu_minusfile.close();
......@@ -134,6 +137,7 @@ void B1EventAction::EndOfEventAction(const G4Event* event)
nu_mu_file.close();
anti_nu_e_file.close();
anti_nu_mu_file.close();
particles_file.close();
}
......
......@@ -79,7 +79,7 @@ void B1PrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent)
fParticleGun->SetParticlePosition(G4ThreeVector(0,0,-100.*m));
fParticleGun->SetParticlePosition(G4ThreeVector(0,0,-1.*km));
fParticleGun->GeneratePrimaryVertex(anEvent);
}
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment