...
 
Commits (3)
......@@ -58,7 +58,8 @@ class B1DetectorConstruction : public G4VUserDetectorConstruction
G4LogicalVolume* GetScoringVolume() const { return fScoringVolume; }
G4LogicalVolume* GetScoringDetector() const { return fScoringDetector; }
G4int return_num_detec() {return number_detectors; }
G4int return_num_detec_longitudinal() {return number_detectors_longitudinal; }
G4int return_num_detec_radial() {return number_detectors_radial;}
G4double get_delta(){return delta;}
G4double get_size() {return size;}
......@@ -87,8 +88,10 @@ class B1DetectorConstruction : public G4VUserDetectorConstruction
G4LogicalVolume* fScoringDetector;
G4LogicalVolume* logicWorld;
std::vector<G4double> positions;
std::vector<G4LogicalVolume* > logicDetector;
G4int number_detectors;
//std::vector<std::vector<G4LogicalVolume* > > logicDetector_array;
std::vector<G4Tubs* >radial_tubs;
G4int number_detectors_longitudinal;
G4int number_detectors_radial;
G4double delta;
G4double size;
G4Tubs* detector;
......
......@@ -10,90 +10,92 @@ y_nu_mu = []
y_anti_nu_e = []
y_anti_nu_mu = []
#number of events in geant4
n_of_events = 2
#choose the graphs you want putting true:
gamma_graph= True
mu_graph = True
neutrino_graph = True
with open("gamma.txt", "r") as csvfile:
plots = csv.reader(csvfile, delimiter=" ")
for row in plots:
y_gamma.append(float(row[0]))
x.append(float(row[1]))
with open("mu+.txt", "r") as csvfile:
plots = csv.reader(csvfile, delimiter=" ")
for row in plots:
y_mu_plus.append(float(row[0]))
with open("mu-.txt", "r") as csvfile:
plots = csv.reader(csvfile, delimiter=" ")
for row in plots:
y_mu_minus.append(float(row[0]))
with open("nu_e.txt", "r") as csvfile:
plots = csv.reader(csvfile, delimiter=" ")
for row in plots:
y_nu_e.append(float(row[0]))
with open("nu_mu.txt", "r") as csvfile:
plots = csv.reader(csvfile, delimiter=" ")
for row in plots:
y_nu_mu.append(float(row[0]))
with open("anti_nu_e.txt", "r") as csvfile:
plots = csv.reader(csvfile, delimiter=" ")
for row in plots:
y_anti_nu_e.append(float(row[0]))
with open("anti_nu_mu.txt", "r") as csvfile:
plots = csv.reader(csvfile, delimiter=" ")
for row in plots:
y_anti_nu_mu.append(float(row[0]))
if (gamma_graph):
#gamma
plt.scatter(x, y_gamma, s=20)
plt.title("Geant4 hadronic shower graphs")
plt.ylabel("Nº gamma")
plt.xlabel("distance (km)")
plt.xlim(xmin=0)
plt.ylim(ymin=0)
plt.show()
for i in range(0, n_of_events):
with open("gamma" + str(i) + ".txt", "r") as csvfile:
plots = csv.reader(csvfile, delimiter=" ")
for row in plots:
y_gamma.append(float(row[0]))
x.append(float(row[1]))
plt.scatter(x, y_gamma, s=20)
plt.title("Geant4 hadronic shower graphs")
plt.ylabel("Nº gamma")
plt.xlabel("distance (km)")
plt.xlim(xmin=0)
plt.ylim(ymin=0)
plt.show()
y_gamma.clear()
x.clear()
if (mu_graph):
#muons
line1 = plt.scatter(x, y_mu_plus, c='b', label='1', s=20, marker = '^')
line2 = plt.scatter(x,y_mu_minus, c='r', label='1', s=20)
plt.title("Geant4 hadronic shower graphs")
plt.ylabel("Nº muons")
plt.xlabel("distance (km)")
plt.xlim(xmin=0)
plt.ylim(ymin=0)
plt.legend((line1, line2), ('mu-', 'mu+'))
plt.show()
#muons
for i in range(0, n_of_events):
with open("mu+" + str(i) + ".txt", "r") as csvfile:
plots = csv.reader(csvfile, delimiter=" ")
for row in plots:
y_mu_plus.append(float(row[0]))
x.append(float(row[1]))
with open("mu-" + str(i) + ".txt", "r") as csvfile:
plots = csv.reader(csvfile, delimiter=" ")
for row in plots:
y_mu_minus.append(float(row[0]))
line1 = plt.scatter(x, y_mu_plus, c='b', label='1', s=20, marker = '^')
line2 = plt.scatter(x,y_mu_minus, c='r', label='1', s=20)
plt.title("Geant4 hadronic shower graphs")
plt.ylabel("Nº muons")
plt.xlabel("distance (km)")
plt.xlim(xmin=0)
plt.ylim(ymin=0)
plt.legend((line1, line2), ('mu-', 'mu+'))
plt.show()
y_mu_plus.clear()
y_mu_minus.clear()
x.clear()
if (neutrino_graph):
#neutrinos
line1 = plt.scatter(x, y_nu_e, c='b', label='1', s=20)
line2 = plt.scatter(x,y_anti_nu_e, c='r', label='1', s=20)
plt.title("Geant4 hadronic shower graphs")
plt.ylabel("Nº electron neutrinos")
plt.xlabel("distance (km)")
plt.xlim(xmin=0)
plt.ylim(ymin=0)
plt.legend((line1, line2), ('nu_e', 'anti_nu_e'))
plt.show()
line1 = plt.scatter(x, y_nu_mu, c='b', label='1', s=20)
line2 = plt.scatter(x,y_anti_nu_mu, c='r', label='1', s=20, marker = '^')
plt.title("Geant4 hadronic shower graphs")
plt.ylabel("Nº muon neutrinos")
plt.xlabel("distance (km)")
plt.xlim(xmin=0)
plt.ylim(ymin=0)
plt.legend((line1, line2), ('nu_mu', 'anti_nu_mu'))
plt.show()
for i in range(0, n_of_events):
with open("anti_nu_e" + str(i) + ".txt", "r") as csvfile:
plots = csv.reader(csvfile, delimiter=" ")
for row in plots:
y_anti_nu_e.append(float(row[0]))
x.append(float(row[1]))
with open("anti_nu_mu" + str(i) + ".txt", "r") as csvfile:
plots = csv.reader(csvfile, delimiter=" ")
for row in plots:
y_anti_nu_mu.append(float(row[0]))
line1 = plt.scatter(x, y_nu_e, c='b', label='1', s=20)
line2 = plt.scatter(x,y_anti_nu_e, c='r', label='1', s=20)
plt.title("Geant4 hadronic shower graphs")
plt.ylabel("Nº electron neutrinos")
plt.xlabel("distance (km)")
plt.xlim(xmin=0)
plt.ylim(ymin=0)
plt.legend((line1, line2), ('nu_e', 'anti_nu_e'))
plt.show()
line1 = plt.scatter(x, y_nu_mu, c='b', label='1', s=20)
line2 = plt.scatter(x,y_anti_nu_mu, c='r', label='1', s=20, marker = '^')
plt.title("Geant4 hadronic shower graphs")
plt.ylabel("Nº muon neutrinos")
plt.xlabel("distance (km)")
plt.xlim(xmin=0)
plt.ylim(ymin=0)
plt.legend((line1, line2), ('nu_mu', 'anti_nu_mu'))
plt.show()
......@@ -78,7 +78,9 @@ G4VPhysicalVolume* B1DetectorConstruction::Construct()
//caracteristicas do cilíndro -> world
G4double raio_i = 0;
G4double raio_e = 50*m;
G4double h_detector = 1*cm;
G4double h_detector = .25*m;
G4double h = 100*m;
G4double theta_0 = 0.*deg;
G4double theta_f = 360.*deg;
......@@ -108,35 +110,75 @@ G4VPhysicalVolume* B1DetectorConstruction::Construct()
fScoringVolume = logicWorld;
G4double safe_distance = 0.02*m;
//Coloque aqui o número de detetores
G4int num_detector = 32;
number_detectors = num_detector;
G4double safe_distance = 2*m;
//Coloque aqui o número de detetores ao longo do cilindro
G4int num_detector_longitudinal = 10;
number_detectors_longitudinal = num_detector_longitudinal;
G4double total_safe = h_detector + safe_distance;
G4double delta_h = (2*h - 2*safe_distance - 4*h_detector)/(num_detector - 1);
//para deixar detetores igualmente espaçados
G4double total_safe = h_detector + safe_distance;
G4double delta_h = (2*h - 2*safe_distance - 4*h_detector)/(num_detector_longitudinal - 1);
G4double initial_pos = -h + total_safe;
G4double final_pos = h - total_safe;
size = h;
delta = delta_h;
for(G4int i = 0; i < num_detector_longitudinal; i++) {
positions.push_back(initial_pos + i*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
"logicDetector");
*/
G4int num_radial_detectors = 5 ;
number_detectors_radial = num_radial_detectors;
G4int delta_r = raio_e/num_radial_detectors;
radial_tubs.reserve(num_radial_detectors);
//inicializando a matriz
std::vector<std::vector<G4LogicalVolume* >> logicDetector_array(num_radial_detectors, std::vector<G4LogicalVolume* >(num_detector_longitudinal));
//logicDetector_array = logicDetectors;
for(G4int i = 0; i < num_radial_detectors; i++) {
radial_tubs[i] = new G4Tubs("radial_detector" + std::to_string(i), (delta_r*i), (delta_r*i + delta_r), h_detector, theta_0, theta_f);
//logicDetector_array[i] = new G4LogicalVolume[num_detector_longitudinal];
for(G4int j = 0; j < num_detector_longitudinal; j++) {
logicDetector_array[i][j] = new G4LogicalVolume(radial_tubs[i], //its solid
world_mat, //its material
"logicRadialDetector" + std::to_string(i) + std::to_string(j));
new G4PVPlacement(0, //no rotation
G4ThreeVector(0,0,positions[j]), //at (0,0,0)
logicDetector_array[i][j], //its logical volume
("Radialdetector" + std::to_string(i) + std::to_string(j)), //its name
logicWorld, //its mother volume
false, //no boolean operation
j, //copy number
checkOverlaps); //overlaps checking
}
}
//Criando os Logical volume
/*
logicDetector.reserve(num_detector);
for(G4int i=0; i < num_detector; i++) {
for(G4int i=0; i < num_detector_longitudinal; i++) {
logicDetector[i] = new G4LogicalVolume(detector, //its solid
world_mat, //its material
"logicDetector" + std::to_string(i + 1));
positions.push_back(initial_pos + i*delta_h);
std::cout << positions[i] << '\n';
new G4PVPlacement(0, //no rotation
......@@ -150,41 +192,57 @@ 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
auto sdman = G4SDManager::GetSDMpointer();
//Associando os detetores à classe sensitive detector
//DEPOIS: Juntar tudo que puder em apenas um "for"
auto sdman = G4SDManager::GetSDMpointer();
for(G4int i = 0; i< num_detector; i++) {
G4String SDname = "SD" + std::to_string(i + 1);
for(G4int i = 0; i< num_radial_detectors; i++) {
for (G4int j = 0; j < num_detector_longitudinal; j++) {
G4String SDname = "SD" + std::to_string(i) + std::to_string(j);
auto sensitive = new B1SD(SDname);
sdman->AddNewDetector(sensitive);
logicDetector[i]->SetSensitiveDetector(sensitive);
logicDetector_array[i][j]->SetSensitiveDetector(sensitive);
}
}
//colocando cor vermelha para detectores
//colocando cor para detectores
std::vector<G4VisAttributes*> visAtributes;
visAtributes.reserve(num_radial_detectors);
G4VisAttributes* worldVisAtt1 = new G4VisAttributes(G4Colour(1.0,0.0,0.0));
worldVisAtt1->SetVisibility(true);
for(G4int i = 0; i < num_detector; i++) {
logicDetector[i]->SetVisAttributes(worldVisAtt1 );
for(G4int i = 0; i < num_radial_detectors; i++) {
visAtributes[i] = new G4VisAttributes(G4Colour(1.0,5.0*i,1.0*i));
visAtributes[i]->SetVisibility(true);
for (G4int j = 0; j < num_detector_longitudinal; j++) {
logicDetector_array[i][j]->SetVisAttributes(visAtributes[i]);
}
}
return physWorld;
}
void B1DetectorConstruction::ConstructSDandField() {
//creanting uniform magnetic field
/*
G4MagneticField* magField = new G4UniformMagField(G4ThreeVector(.0000005*LT,0.,0.));
G4FieldManager* FieldMgr =new G4FieldManager(magField);
logicWorld->SetFieldManager(FieldMgr, true);
*/
}
......
......@@ -81,41 +81,44 @@ void B1EventAction::EndOfEventAction(const G4Event* event)
G4SDManager * SDman = G4SDManager::GetSDMpointer();
G4HCofThisEvent* HCE = event->GetHCofThisEvent();
G4int num_event = event->GetEventID();
detectorConstruction = (B1DetectorConstruction*) G4RunManager::GetRunManager()->GetUserDetectorConstruction();
G4int num = detectorConstruction->return_num_detec();
dec_pos = detectorConstruction->GetPositions();
G4int num_longitudinal = detectorConstruction->return_num_detec_longitudinal();
G4int num_radial = detectorConstruction->return_num_detec_radial();
dec_pos = detectorConstruction->GetPositions();
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" + std::to_string(num_event) + ".txt");
std::ofstream mu_plusfile("mu+" + std::to_string(num_event) + ".txt");
std::ofstream mu_minusfile("mu-" + std::to_string(num_event) + ".txt");
std::ofstream nu_e_file("nu_e" + std::to_string(num_event) + ".txt");
std::ofstream nu_mu_file("nu_mu" + std::to_string(num_event) + ".txt");
std::ofstream anti_nu_e_file("anti_nu_e" + std::to_string(num_event) + ".txt");
std::ofstream anti_nu_mu_file("anti_nu_mu" + std::to_string(num_event) + ".txt");
std::vector<G4int> col;
std::vector<B1HitsCollection*> HitsCol;
std::vector<std::vector <G4int> > col(num_radial, std::vector <G4int>(num_longitudinal));
std::vector<std::vector <B1HitsCollection*> > HitsCol(num_radial, std::vector<B1HitsCollection*>(num_longitudinal));
col.reserve(num);
HitsCol.reserve(num);
for(G4int i=0; i < num_radial; i++) {
for(G4int j = 0; j < num_longitudinal; j++ ) {
for(G4int i=0; i < num; i++) {
col[i] = SDman->GetCollectionID("SD" + std::to_string(i + 1));
col[i][j] = SDman->GetCollectionID("SD" + std::to_string(i) + std::to_string(j));
if(HCE) {
HitsCol[i] = (B1HitsCollection*)(HCE->GetHC(col[i]));
HitsCol[i][j] = (B1HitsCollection*)(HCE->GetHC(col[i][j]));
}
if(HitsCol[i]) {
if(HitsCol[i][j]) {
G4double pos_detec = SetPosition(i + 1);
int n_hit = HitsCol[i]->entries();
G4cout << "My detector " + std::to_string(i+1) +" has " << n_hit << "hits" << G4endl;
int n_hit = HitsCol[i][j]->entries();
G4cout << "My detector " + std::to_string(i) + std::to_string(j) + " has " << n_hit << "hits" << G4endl;
B1Hits* hit = new B1Hits;
std::map<const G4String, int> fparticles;
for(int i1 = 0; i1 < n_hit; i1++) {
B1Hits* hit = (*HitsCol[i])[i1];
B1Hits* hit = (*HitsCol[i][j])[i1];
const G4String name = hit->getParticleInTarget();
//G4cout << name << G4endl;
fparticles[name]++;
......@@ -125,6 +128,7 @@ void B1EventAction::EndOfEventAction(const G4Event* event)
fparticles.clear();
}
}
}
gammafile.close();
mu_plusfile.close();
......@@ -133,6 +137,7 @@ void B1EventAction::EndOfEventAction(const G4Event* event)
nu_mu_file.close();
anti_nu_e_file.close();
anti_nu_mu_file.close();
}
......@@ -140,10 +145,11 @@ void B1EventAction::PrintParticles(std::map<const G4String, int>& container, std
std::ofstream& nu_e, std::ofstream& nu_mu, std::ofstream& anti_nu_e, std::ofstream& anti_nu_mu, 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;
for(it = container.begin() ;it != container.end(); it ++)
{
G4bool isGamma, ismu_minus, ismu_plus, isnu_e, isnu_mu, isanti_nu_e, isanti_nu_mu = false;
G4cout << "N " << it->first << " : " << it->second << G4endl;
//FUTURO -> substituir por switch;
if(it->first == "gamma") {
......@@ -169,29 +175,30 @@ void B1EventAction::PrintParticles(std::map<const G4String, int>& container, std
anti_nu_mu << it->second << " " << position << "\n";
isanti_nu_mu = true;
}
}
if (!isGamma) {
gamma << "0 " << position << "\n";
}
if (!ismu_plus) {
mu_plus << "0 " << position << "\n";
}
if (!ismu_minus) {
mu_minus << "0 " << position << "\n";
}
if (!isnu_e) {
nu_e << "0 " << position << "\n";
}
if (!isnu_mu) {
nu_mu << "0 " << position << "\n";
}
if (!isanti_nu_e) {
anti_nu_e << "0 " << position << "\n";
}
if (!isanti_nu_mu) {
anti_nu_mu << "0 " << position << "\n";
}
if (!isGamma) {
gamma << "0 " << position << "\n";
}
if (!ismu_plus) {
mu_plus << "0 " << position << "\n";
}
if (!ismu_minus) {
mu_minus << "0 " << position << "\n";
}
if (!isnu_e) {
nu_e << "0 " << position << "\n";
}
if (!isnu_mu) {
nu_mu << "0 " << position << "\n";
}
if (!isanti_nu_e) {
anti_nu_e << "0 " << position << "\n";
}
if (!isanti_nu_mu) {
anti_nu_mu << "0 " << position << "\n";
}
}
}
void B1EventAction::WriteHistogram(const G4String name, G4int Detec) {
......@@ -226,4 +233,6 @@ G4double B1EventAction::SetPosition(G4int detec) {
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......