...
 
Commits (8)
......@@ -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;
}
......
......@@ -35,6 +35,10 @@
#include "globals.hh"
#include "G4VisAttributes.hh"
#include <vector>
#include <cmath>
#include <iostream>
#include <fstream>
class B1SD;
class G4VPhysicalVolume;
......@@ -59,10 +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 h);
std::vector<G4double>& GetPositions() {
return positions;}
......@@ -88,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;
......
......@@ -52,14 +52,16 @@ class B1EventAction : public G4UserEventAction
virtual void BeginOfEventAction(const G4Event* event);
virtual void EndOfEventAction(const G4Event* event);
void PrintParticles(std::map<const G4String, int>& container, std::ofstream& gamma, std::ofstream& mu_minus, std::ofstream& mu_plus,
std::ofstream& nu_e, std::ofstream& nu_mu, std::ofstream& anti_nu_e, std::ofstream& anti_nu_mu, G4double position);
void PrintParticles(std::map<const G4String, int>& container, std::ofstream& gamma, std::ofstream& mu,
std::ofstream& pi, std::ofstream& e, G4double position);
void WriteHistogram(const G4String name, G4int Detec);
void AddEdep(G4double edep) { fEdep += edep; };
G4double SetPosition(G4int detec);
private:
B1RunAction* fRunAction;
G4int mu_event;
G4int y;
G4double fEdep;
std::vector<G4double> dec_pos;
B1DetectorConstruction* detectorConstruction;
......
This diff is collapsed.
//
//
// ********************************************************************
// * License and Disclaimer *
// * *
......@@ -87,19 +87,18 @@ void B1EventAction::EndOfEventAction(const G4Event* event)
half_height = detectorConstruction->get_size();
variation = detectorConstruction->get_delta();
std::ofstream gammafile("gamma.txt");
std::ofstream mu_plusfile("mu+.txt");
std::ofstream mu_minusfile("mu-.txt");
std::ofstream nu_e_file("nu_e.txt");
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 gammafile("gamma.csv");
std::ofstream mu_file("mu.csv");
std::ofstream e_file("e.txt");
std::ofstream particles_file("particles.csv");
std::ofstream pi_file("pi.csv");
std::vector<G4int> col;
std::vector<B1HitsCollection*> HitsCol;
col.reserve(num);
HitsCol.reserve(num);
col.reserve(num + 1);
HitsCol.reserve(num + 1);
for(G4int i=0; i < num; i++) {
col[i] = SDman->GetCollectionID("SD" + std::to_string(i + 1));
......@@ -109,9 +108,13 @@ 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);
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/(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++) {
......@@ -121,29 +124,59 @@ void B1EventAction::EndOfEventAction(const G4Event* event)
fparticles[name]++;
//WriteHistogram(name, 1);
}
PrintParticles(fparticles, gammafile, mu_plusfile, mu_minusfile, nu_e_file, nu_mu_file, anti_nu_e_file, anti_nu_mu_file, pos_detec);
PrintParticles(fparticles, gammafile, mu_file, pi_file, e_file, pos_detec);
fparticles.clear();
}
}
/*
G4int mu_num = 0;
col[num + 1] = SDman->GetCollectionID("SDs");
HitsCol[num + 1] = (B1HitsCollection*)(HCE->GetHC(col[num + 1]));
int n_hit = HitsCol[num + 1]->entries();
G4cout << "surface detector has " << n_hit << "hits" << G4endl;
B1Hits* hit = new B1Hits;
for(int i1 = 0; i1 < n_hit; i1++) {
B1Hits* hit = (*HitsCol[num + 1])[i1];
const G4String name = hit->getParticleInTarget();
if (name == "mu+" || name == "mu-") {
G4cout << "muuuuuuuuu" << G4endl;
mu_num++;
}
}
G4cout << mu_num << G4endl;
mu_event = mu_event + mu_num;
y = y + 1;
mu_number << y << " " << mu_event;
mu_num = 0;
*/
gammafile.close();
mu_plusfile.close();
mu_minusfile.close();
nu_e_file.close();
nu_mu_file.close();
anti_nu_e_file.close();
anti_nu_mu_file.close();
mu_file.close();
particles_file.close();
pi_file.close();
e_file.close();
}
void B1EventAction::PrintParticles(std::map<const G4String, int>& container, std::ofstream& gamma, std::ofstream& mu_minus, std::ofstream& mu_plus,
std::ofstream& nu_e, std::ofstream& nu_mu, std::ofstream& anti_nu_e, std::ofstream& anti_nu_mu, G4double position) {
void B1EventAction::PrintParticles(std::map<const G4String, int>& container, std::ofstream& gamma, std::ofstream& mu,
std::ofstream& pi, std::ofstream& e, G4double position) {
std::map<const G4String, int>::iterator it;
// G4cout << "Número de párticulas identificadas no detetor: " << G4endl;
G4bool isGamma, ismu_minus, ismu_plus, isnu_e, isnu_mu, isanti_nu_e, isanti_nu_mu = false;
G4bool isGamma, ismu, is_e, is_pi = false;
G4int n_e_minus, n_e_pl, n_mu_minus, n_mu_plus, n_pi_plus, n_pi_minus;
n_e_minus = 0;
n_e_pl = 0;
n_mu_minus = 0;
n_mu_plus = 0;
n_pi_plus = 0;
n_pi_minus = 0;
for(it = container.begin() ;it != container.end(); it ++)
{
G4cout << n_e_pl << G4endl;
G4cout << "N " << it->first << " : " << it->second << G4endl;
//FUTURO -> substituir por switch;
if(it->first == "gamma") {
......@@ -151,47 +184,56 @@ void B1EventAction::PrintParticles(std::map<const G4String, int>& container, std
gamma << it->second << " " << position << "\n";
isGamma = true;
} else if (it->first == "mu+") {
mu_plus << it->second << " " << position << "\n";
ismu_plus= true;
n_mu_plus = it->second;
ismu = true;
} else if (it->first == "mu-") {
mu_minus << it->second << " " << position << "\n";
ismu_minus = true;
} else if (it->first == "nu_e") {
nu_e << it->second << " " << position << "\n";
isnu_e = true;
} else if (it->first == "nu_mu") {
nu_mu << it->second << " " << position << "\n";
isnu_mu = true;
} else if (it->first == "anti_nu_e") {
anti_nu_e << it->second << " " << position << "\n";
isanti_nu_e = true;
} else if (it->first == "anti_nu_mu") {
anti_nu_mu << it->second << " " << position << "\n";
isanti_nu_mu = true;
n_mu_minus = it->second;
ismu = true;
} else if (it->first == "e+") {
n_e_pl = it->second;
G4cout << "aa" << n_e_pl << G4endl;
is_e = true;
} else if (it->first == "e-") {
n_e_minus = it->second;
G4cout << "aa" << n_e_minus << G4endl;
is_e = true;
} else if (it->first == "pi+") {
n_pi_plus = it->second;
is_pi = true;
} else if (it->first == "pi-") {
n_pi_minus = it->second;
is_pi = true;
}
}
}
if (!isGamma) {
gamma << "0 " << position << "\n";
}
if (!ismu_plus) {
mu_plus << "0 " << position << "\n";
}
if (!ismu_minus) {
mu_minus << "0 " << position << "\n";
gamma << 0 << " " << position << "\n";
}
if (!isnu_e) {
nu_e << "0 " << position << "\n";
}
if (!isnu_mu) {
nu_mu << "0 " << position << "\n";
if (!ismu) {
mu << 0 << " " << position << "\n";
} else {
G4int n = n_mu_minus + n_mu_plus;
mu << n << " " << position << "\n";
}
if (!isanti_nu_e) {
anti_nu_e << "0 " << position << "\n";
if(!is_e) {
e << 0 << " " << position << "\n";
} else {
G4int n = n_e_pl + n_e_minus;
G4cout << n_e_pl << G4endl;
G4cout << n_e_minus << G4endl;
G4cout << n << G4endl;
e << n << " " << position << "\n";
}
if (!isanti_nu_mu) {
anti_nu_mu << "0 " << position << "\n";
if(!is_pi) {
pi << 0 << " " << position << "\n";
} else {
G4int n = n_pi_plus + n_pi_minus;
pi << n << " " << position << "\n";
}
}
void B1EventAction::WriteHistogram(const G4String name, G4int Detec) {
......
......@@ -79,7 +79,7 @@ void B1PrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent)
fParticleGun->SetParticlePosition(G4ThreeVector(0,0,-100.*m));
fParticleGun->SetParticlePosition(G4ThreeVector(0,0,-50.*km));
fParticleGun->GeneratePrimaryVertex(anEvent);
}
......
......@@ -5,7 +5,7 @@
# Use these open statements to open selected visualization
#
# Use this open statement to create an OpenGL view:
/vis/open OGL 600x600-0+0
/vis/open OGL
#
# Use this open statement to create an OpenInventor view:
#/vis/open OI
......