...
 
Commits (2)
......@@ -63,13 +63,12 @@ class B1DetectorConstruction : public G4VUserDetectorConstruction
G4LogicalVolume* GetScoringVolume() const { return fScoringVolume; }
G4LogicalVolume* GetScoringDetector() const { return fScoringDetector; }
G4int return_num_detec() {return number_detectors; }
G4double get_mass_overburden(G4int i, G4double delta);
G4double get_delta(){return delta;}
G4double get_size() {return size;}
G4double GetDetecPos(G4int num_detec);
G4double density_function(G4double a_i, G4double b_i, G4double c_i, G4double h) {
return (a_i + b_i*exp(-h/c_i));
}
G4double density_function(G4double h);
std::vector<G4double>& GetPositions() {
return positions;}
......@@ -95,6 +94,10 @@ class B1DetectorConstruction : public G4VUserDetectorConstruction
G4LogicalVolume* logicWorld;
std::vector<G4double> positions;
std::vector<G4LogicalVolume* > logicDetector;
std::vector<G4double> a_atm = {-186.5562, -94.919, 0.61289, 0.0, 0.01128292};
std::vector<G4double> b_atm = {1222.6562, 1144.9069, 1305.5948, 540.1778, 1.0};
std::vector<G4double> c_atm = {994186.38, 878153.55, 636143.04, 772170., 1.0e9};
std::vector<G4double> h_atm = {0, 4.0e5, 1.0e6, 4.0e6, 1.0e7};
G4int number_detectors;
G4double delta;
G4double size;
......
......@@ -81,7 +81,7 @@ G4VPhysicalVolume* B1DetectorConstruction::Construct()
G4double raio_i = 0;
G4double raio_e = 1*km;
G4double h_detector = 1*cm;
G4double h = 1*km;
G4double h = 40*km;
G4double theta_0 = 0.*deg;
G4double theta_f = 360.*deg;
//making world of vacuum
......@@ -115,41 +115,37 @@ G4VPhysicalVolume* B1DetectorConstruction::Construct()
//making air materiais
//Nitrogen
G4int num_layers = 200;
G4int num_layers = 250;
G4double layer_delta_h = 2*h/num_layers;
std::vector<G4Material *> layers_material;
layers_material.reserve(num_layers);
std::ofstream density_file("density.csv");
std::ofstream mass_over("mass.csv");
//Corsika density parameters
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(-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*g/cm2, 1144.9069*g/cm2, 878153.55*cm, atm_position);
} else if (atm_position >= 10*km && atm_position < 40*km) {
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) {
G4double gm2;
density = density_function(0.0, 540.1778, 772170.16, atm_position);
}
G4cout << density/(g/cm2) << G4endl;
density_file << density/(g/cm2) << " " << atm_position/(km) << "\n";
density = density_function(atm_position);
gm2 = get_mass_overburden(j, layer_delta_h);
density_file << density/(kg/m3) << " " << atm_position/(km) << "\n";
mass_over << gm2 << " " << atm_position/(km) << "\n";
G4cout << j << G4endl;
layers_material[j] = new G4Material(name="Air" + std::to_string(j), density/1000000, 3);
layers_material[j] = new G4Material(name="Air" + std::to_string(j), density, 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();
mass_over.close();
std::vector<G4LogicalVolume* > air_layers;
air_layers.reserve(num_layers);
......@@ -260,5 +256,59 @@ void B1DetectorConstruction::ConstructSDandField() {
*/
}
G4double B1DetectorConstruction::density_function(G4double h) {
G4int l; //layer
G4double res;
G4cout << h << G4endl;
if (h >= h_atm[0]*cm && h < h_atm[1]*cm) {
l = 0;
} else if (h >= h_atm[1]*cm && h < h_atm[2]*cm) {
l = 1;
} else if (h >= h_atm[2]*cm && h < h_atm[3]*cm) {
l = 2;
} else if (h >= h_atm[3]*cm && h < h_atm[4]*cm) {
l = 3;
} else if (h >= h_atm[5]*cm) {
l = 4;
}
if (l == 4) {
res = b_atm[4]*cm2 / c_atm[4]*cm;
} else {
G4double div = b_atm[l]/c_atm[l];
res = div*exp(-h/(c_atm[l]*cm));
}
return res*g/cm3;
}
G4double B1DetectorConstruction::get_mass_overburden(G4int i, G4double delta) {
G4double h = (delta*i);
G4int l;
G4double res;
if (h >= h_atm[0]*cm && h < h_atm[1]*cm) {
l = 0;
} else if (h >= h_atm[1]*cm && h < h_atm[2]*cm) {
l = 1;
} else if (h >= h_atm[2]*cm && h < h_atm[3]*cm) {
l = 2;
} else if (h >= h_atm[3]*cm && h < h_atm[4]*cm) {
l = 3;
} else if (h >= h_atm[5]*cm) {
l = 4;
}
if (l == 4) {
res = (a_atm[4]*cm2) - ((b_atm[4]*cm2)/(c_atm[4]*cm))*h;
} else {
res = (a_atm[l]) + (b_atm[l])*exp(-h/(c_atm[l]*cm));
G4cout << "TT" << G4endl;
G4cout << res << G4endl;
}
return res;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
......@@ -112,10 +112,12 @@ void B1EventAction::EndOfEventAction(const G4Event* event)
if(HitsCol[i]) {
G4double calc = 2*half_height - (variation/2) - (variation)*i;
G4double pos_detec = calc/(km);
//G4double pos_detec = calc/(km);
G4double pos_detec = detectorConstruction->get_mass_overburden(i, variation);
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";
particles_file << n_hit << " " << pos_detec/(g/cm2) << "\n";
G4cout << pos_detec/(g/cm2) << G4endl;
B1Hits* hit = new B1Hits;
std::map<const G4String, int> fparticles;
for(int i1 = 0; i1 < n_hit; i1++) {
......
......@@ -79,7 +79,7 @@ void B1PrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent)
fParticleGun->SetParticlePosition(G4ThreeVector(0,0,-1.*km));
fParticleGun->SetParticlePosition(G4ThreeVector(0,0,-40.*km));
fParticleGun->GeneratePrimaryVertex(anEvent);
}
......