Commit d0929469 authored by MARCOS ANTONIO DE OLIVEIRA DEROS's avatar MARCOS ANTONIO DE OLIVEIRA DEROS
Browse files

added more detectors and write and file with the output

parent 976b9581
......@@ -34,6 +34,8 @@
#include "globals.hh"
#include <map>
#include "G4ParticleDefinition.hh"
#include <iostream>
#include <fstream>
class B1RunAction;
......@@ -48,8 +50,8 @@ class B1EventAction : public G4UserEventAction
virtual void BeginOfEventAction(const G4Event* event);
virtual void EndOfEventAction(const G4Event* event);
void PrintParticles(std::map<const G4String, int>& container);
void PrintParticles(std::map<const G4String, int>& container, std::ofstream& gamma, std::ofstream& mu_minus, std::ofstream& mu_plus, G4double position);
void WriteHistogram(const G4String name, G4int Detec);
void AddEdep(G4double edep) { fEdep += edep; };
private:
......
import matplotlib.pyplot as plt
import csv
x = []
y_gamma = []
y_mu_minus = []
y_mu_plus = []
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]))
plt.scatter(x, y_gamma)
plt.title("Geant4 hadronic shower graphs")
plt.ylabel("Nº gamma")
plt.xlabel("distance (km)")
plt.xlim(xmin=0)
plt.ylim(ymin=0)
plt.show()
line1 = plt.scatter(x, y_mu_plus, c='b', label='1')
line2 = plt.scatter(x,y_mu_minus, c='r', label='1')
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()
......@@ -75,9 +75,9 @@ G4VPhysicalVolume* B1DetectorConstruction::Construct()
//
//caracteristicas do cilíndro -> world
G4double raio_i = 0;
G4double raio_e = 80.*m;
G4double h_detector = 1.*m;
G4double h = 100.*m;
G4double raio_e = 9*m;
G4double h_detector = 2.*cm;
G4double h = 9*m;
G4double theta_0 = 0.*deg;
G4double theta_f = 360.*deg;
G4Material* world_mat = nist->FindOrBuildMaterial("G4_AIR");
......@@ -108,34 +108,388 @@ G4VPhysicalVolume* B1DetectorConstruction::Construct()
G4Tubs* detector = new G4Tubs("detector", raio_i, raio_e, h_detector, theta_0, theta_f);
G4LogicalVolume* logicDetector = new G4LogicalVolume(detector,
G4LogicalVolume* logicDetector1 = new G4LogicalVolume(detector,
world_mat,
"detector");
"detector1");
//associando à um sensitive G4VSensitiveDetecto
G4String SDname;
auto sensitive = new B1SD(SDname = "Sensivel");
auto sdman = G4SDManager::GetSDMpointer();
sdman->AddNewDetector(sensitive);
logicDetector->SetSensitiveDetector(sensitive);
G4LogicalVolume* logicDetector2 = new G4LogicalVolume(detector,
world_mat,
"detector2");
G4LogicalVolume* logicDetector3 = new G4LogicalVolume(detector,
world_mat,
"detector3");
G4LogicalVolume* logicDetector4 = new G4LogicalVolume(detector,
world_mat,
"detector4");
G4LogicalVolume* logicDetector5 = new G4LogicalVolume(detector,
world_mat,
"detector5");
G4LogicalVolume* logicDetector6 = new G4LogicalVolume(detector,
world_mat,
"detector6");
G4LogicalVolume* logicDetector7 = new G4LogicalVolume(detector,
world_mat,
"detector7");
G4LogicalVolume* logicDetector8 = new G4LogicalVolume(detector,
world_mat,
"detector8");
G4LogicalVolume* logicDetector9 = new G4LogicalVolume(detector,
world_mat,
"detector9");
G4LogicalVolume* logicDetector10 = new G4LogicalVolume(detector,
world_mat,
"detector10");
G4VPhysicalVolume* physDetector = new G4PVPlacement(0,
G4ThreeVector(0,0,99.*m),
logicDetector,
"detector",
G4LogicalVolume* logicDetector11 = new G4LogicalVolume(detector,
world_mat,
"detector11");
G4LogicalVolume* logicDetector12 = new G4LogicalVolume(detector,
world_mat,
"detector12");
G4LogicalVolume* logicDetector13 = new G4LogicalVolume(detector,
world_mat,
"detector13");
G4LogicalVolume* logicDetector14 = new G4LogicalVolume(detector,
world_mat,
"detector14");
G4LogicalVolume* logicDetector15 = new G4LogicalVolume(detector,
world_mat,
"detector15");
G4LogicalVolume* logicDetector16 = new G4LogicalVolume(detector,
world_mat,
"detector16");
G4LogicalVolume* logicDetector17 = new G4LogicalVolume(detector,
world_mat,
"detector17");
G4LogicalVolume* logicDetector18 = new G4LogicalVolume(detector,
world_mat,
"detector18");
G4LogicalVolume* logicDetector19 = new G4LogicalVolume(detector,
world_mat,
"detector19");
//colocando mais detetores
new G4PVPlacement(0,
G4ThreeVector(0,0,-8.99*m),
logicDetector1,
"detector1",
logicWorld,
false,
1,
checkOverlaps);
fScoringDetector = logicDetector;
new G4PVPlacement(0,
G4ThreeVector(0,0,-7.99*m),
logicDetector2,
"detector2",
logicWorld,
false,
1,
checkOverlaps);
new G4PVPlacement(0,
G4ThreeVector(0,0,-6.99*m),
logicDetector3,
"detector3",
logicWorld,
false,
1,
checkOverlaps);
new G4PVPlacement(0,
G4ThreeVector(0,0,-5.99*m),
logicDetector4,
"detector4",
logicWorld,
false,
1,
checkOverlaps);
new G4PVPlacement(0,
G4ThreeVector(0,0,-4.99*m),
logicDetector5,
"detector5",
logicWorld,
false,
1,
checkOverlaps);
new G4PVPlacement(0,
G4ThreeVector(0,0,-3.99*m),
logicDetector6,
"detector6",
logicWorld,
false,
1,
checkOverlaps);
new G4PVPlacement(0,
G4ThreeVector(0,0,-2.99*m),
logicDetector7,
"detector7",
logicWorld,
false,
1,
checkOverlaps);
new G4PVPlacement(0,
G4ThreeVector(0,0,-1.99*m),
logicDetector8,
"detector8",
logicWorld,
false,
1,
checkOverlaps);
new G4PVPlacement(0,
G4ThreeVector(0,0,-0.99*m),
logicDetector9,
"detector9",
logicWorld,
false,
1,
checkOverlaps);
new G4PVPlacement(0,
G4ThreeVector(0,0,0.01*m),
logicDetector10,
"detector10",
logicWorld,
false,
1,
checkOverlaps);
new G4PVPlacement(0,
G4ThreeVector(0,0,1.01*m),
logicDetector11,
"detector11",
logicWorld,
false,
1,
checkOverlaps);
new G4PVPlacement(0,
G4ThreeVector(0,0,2.01*m),
logicDetector12,
"detector12",
logicWorld,
false,
1,
checkOverlaps);
new G4PVPlacement(0,
G4ThreeVector(0,0,3.01*m),
logicDetector13,
"detector13",
logicWorld,
false,
1,
checkOverlaps);
new G4PVPlacement(0,
G4ThreeVector(0,0,4.01*m),
logicDetector14,
"detector14",
logicWorld,
false,
1,
checkOverlaps);
new G4PVPlacement(0,
G4ThreeVector(0,0,5.01*m),
logicDetector15,
"detector15",
logicWorld,
false,
1,
checkOverlaps);
new G4PVPlacement(0,
G4ThreeVector(0,0,6.01*m),
logicDetector16,
"detector16",
logicWorld,
false,
1,
checkOverlaps);
new G4PVPlacement(0,
G4ThreeVector(0,0,7.01*m),
logicDetector17,
"detector17",
logicWorld,
false,
1,
checkOverlaps);
new G4PVPlacement(0,
G4ThreeVector(0,0,8.01*m),
logicDetector18,
"detector18",
logicWorld,
false,
1,
checkOverlaps);
new G4PVPlacement(0,
G4ThreeVector(0,0,8.99*m),
logicDetector19,
"detector19",
logicWorld,
false,
1,
checkOverlaps);
//fScoringDetector = logicDetector;
//associando à um sensitive G4VSensitiveDetecto
auto sdman = G4SDManager::GetSDMpointer();
G4String SDname1;
auto sensitive1 = new B1SD(SDname1 = "SD1");
sdman->AddNewDetector(sensitive1);
logicDetector1->SetSensitiveDetector(sensitive1);
G4String SDname2;
auto sensitive2 = new B1SD(SDname2 = "SD2");
sdman->AddNewDetector(sensitive2);
logicDetector2->SetSensitiveDetector(sensitive2);
G4String SDname3;
auto sensitive3 = new B1SD(SDname3 = "SD3");
sdman->AddNewDetector(sensitive3);
logicDetector3->SetSensitiveDetector(sensitive3);
G4String SDname4;
auto sensitive4 = new B1SD(SDname4 = "SD4");
sdman->AddNewDetector(sensitive4);
logicDetector4->SetSensitiveDetector(sensitive4);
G4String SDname5;
auto sensitive5 = new B1SD(SDname5 = "SD5");
sdman->AddNewDetector(sensitive5);
logicDetector5->SetSensitiveDetector(sensitive5);
G4String SDname6;
auto sensitive6 = new B1SD(SDname6 = "SD6");
sdman->AddNewDetector(sensitive6);
logicDetector6->SetSensitiveDetector(sensitive6);
G4String SDname7;
auto sensitive7 = new B1SD(SDname7 = "SD7");
sdman->AddNewDetector(sensitive7);
logicDetector7->SetSensitiveDetector(sensitive7);
G4String SDname8;
auto sensitive8 = new B1SD(SDname8 = "SD8");
sdman->AddNewDetector(sensitive8);
logicDetector8->SetSensitiveDetector(sensitive8);
G4String SDname9;
auto sensitive9 = new B1SD(SDname9 = "SD9");
sdman->AddNewDetector(sensitive9);
logicDetector9->SetSensitiveDetector(sensitive9);
G4String SDname10;
auto sensitive10 = new B1SD(SDname10 = "SD10");
sdman->AddNewDetector(sensitive10);
logicDetector10->SetSensitiveDetector(sensitive10);
G4String SDname11;
auto sensitive11 = new B1SD(SDname11 = "SD11");
sdman->AddNewDetector(sensitive11);
logicDetector11->SetSensitiveDetector(sensitive11);
G4String SDname12;
auto sensitive12 = new B1SD(SDname12 = "SD12");
sdman->AddNewDetector(sensitive12);
logicDetector12->SetSensitiveDetector(sensitive12);
G4String SDname13;
auto sensitive13 = new B1SD(SDname13 = "SD13");
sdman->AddNewDetector(sensitive13);
logicDetector13->SetSensitiveDetector(sensitive13);
G4String SDname14;
auto sensitive14 = new B1SD(SDname14 = "SD14");
sdman->AddNewDetector(sensitive14);
logicDetector14->SetSensitiveDetector(sensitive14);
G4String SDname15;
auto sensitive15 = new B1SD(SDname15 = "SD15");
sdman->AddNewDetector(sensitive15);
logicDetector15->SetSensitiveDetector(sensitive15);
G4String SDname16;
auto sensitive16 = new B1SD(SDname16 = "SD16");
sdman->AddNewDetector(sensitive16);
logicDetector16->SetSensitiveDetector(sensitive16);
G4String SDname17;
auto sensitive17 = new B1SD(SDname17 = "SD17");
sdman->AddNewDetector(sensitive17);
logicDetector17->SetSensitiveDetector(sensitive17);
G4String SDname18;
auto sensitive18 = new B1SD(SDname18= "SD18");
sdman->AddNewDetector(sensitive18);
logicDetector18->SetSensitiveDetector(sensitive18);
G4String SDname19;
auto sensitive19 = new B1SD(SDname19= "SD19");
sdman->AddNewDetector(sensitive19);
logicDetector19->SetSensitiveDetector(sensitive19);
//colocando cor vermelha para detector
G4VisAttributes* worldVisAtt1 = new G4VisAttributes(G4Colour(1.0,0.0,0.0));
worldVisAtt1->SetVisibility(true);
logicDetector->SetVisAttributes(worldVisAtt1);
logicDetector1->SetVisAttributes(worldVisAtt1);
logicDetector2->SetVisAttributes(worldVisAtt1);
logicDetector3->SetVisAttributes(worldVisAtt1);
logicDetector4->SetVisAttributes(worldVisAtt1);
logicDetector5->SetVisAttributes(worldVisAtt1);
logicDetector6->SetVisAttributes(worldVisAtt1);
logicDetector7->SetVisAttributes(worldVisAtt1);
logicDetector8->SetVisAttributes(worldVisAtt1);
logicDetector9->SetVisAttributes(worldVisAtt1);
logicDetector10->SetVisAttributes(worldVisAtt1);
logicDetector11->SetVisAttributes(worldVisAtt1);
logicDetector12->SetVisAttributes(worldVisAtt1);
logicDetector13->SetVisAttributes(worldVisAtt1);
logicDetector14->SetVisAttributes(worldVisAtt1);
logicDetector15->SetVisAttributes(worldVisAtt1);
logicDetector16->SetVisAttributes(worldVisAtt1);
logicDetector17->SetVisAttributes(worldVisAtt1);
logicDetector18->SetVisAttributes(worldVisAtt1);
logicDetector19->SetVisAttributes(worldVisAtt1);
return physWorld;
}
......
......@@ -39,7 +39,10 @@
#include "G4SDManager.hh"
#include "G4SystemOfUnits.hh"
#include "G4ios.hh"
#include "B1Analysis.hh"
#include <iostream>
#include <fstream>
using namespace std;
......@@ -71,43 +74,444 @@ void B1EventAction::EndOfEventAction(const G4Event* event)
// accumulate statistics in run action
fRunAction->AddEdep(fEdep);
static int CHCID = -1;
if (CHCID < 0) {
CHCID = G4SDManager::GetSDMpointer()->GetCollectionID("HitCollection");//Descobrir isso
}
G4SDManager * SDman = G4SDManager::GetSDMpointer();
G4String detName;
auto Col1 = SDman->GetCollectionID("SD1");
auto Col2 = SDman->GetCollectionID("SD2");
auto Col3 = SDman->GetCollectionID("SD3");
auto Col4 = SDman->GetCollectionID("SD4");
auto Col5 = SDman->GetCollectionID("SD5");
auto Col6 = SDman->GetCollectionID("SD6");
auto Col7 = SDman->GetCollectionID("SD7");
auto Col8 = SDman->GetCollectionID("SD8");
auto Col9 = SDman->GetCollectionID("SD9");
auto Col10 = SDman->GetCollectionID("SD10");
auto Col11 = SDman->GetCollectionID("SD11");
auto Col12 = SDman->GetCollectionID("SD12");
auto Col13 = SDman->GetCollectionID("SD13");
auto Col14 = SDman->GetCollectionID("SD14");
auto Col15 = SDman->GetCollectionID("SD15");
auto Col16 = SDman->GetCollectionID("SD16");
auto Col17 = SDman->GetCollectionID("SD17");
auto Col18 = SDman->GetCollectionID("SD18");
auto Col19 = SDman->GetCollectionID("SD19");
G4HCofThisEvent* HCE = event->GetHCofThisEvent();
B1HitsCollection* HitsCol = 0;
B1HitsCollection* HitsCol1 = 0;
B1HitsCollection* HitsCol2 = 0;
B1HitsCollection* HitsCol3 = 0;
B1HitsCollection* HitsCol4 = 0;
B1HitsCollection* HitsCol5 = 0;
B1HitsCollection* HitsCol6 = 0;
B1HitsCollection* HitsCol7 = 0;
B1HitsCollection* HitsCol8 = 0;
B1HitsCollection* HitsCol9 = 0;
B1HitsCollection* HitsCol10 = 0;
B1HitsCollection* HitsCol11 = 0;
B1HitsCollection* HitsCol12 = 0;
B1HitsCollection* HitsCol13 = 0;
B1HitsCollection* HitsCol14 = 0;
B1HitsCollection* HitsCol15 = 0;
B1HitsCollection* HitsCol16 = 0;
B1HitsCollection* HitsCol17 = 0;
B1HitsCollection* HitsCol18 = 0;
B1HitsCollection* HitsCol19 = 0;
if(HCE) {
HitsCol = (B1HitsCollection*)(HCE->GetHC(CHCID));
HitsCol1 = (B1HitsCollection*)(HCE->GetHC(Col1));
HitsCol2 = (B1HitsCollection*)(HCE->GetHC(Col2));
HitsCol3 = (B1HitsCollection*)(HCE->GetHC(Col3));
HitsCol4 = (B1HitsCollection*)(HCE->GetHC(Col4));
HitsCol5 = (B1HitsCollection*)(HCE->GetHC(Col5));
HitsCol6 = (B1HitsCollection*)(HCE->GetHC(Col6));
HitsCol7 = (B1HitsCollection*)(HCE->GetHC(Col7));
HitsCol8 = (B1HitsCollection*)(HCE->GetHC(Col8));
HitsCol9 = (B1HitsCollection*)(HCE->GetHC(Col9));
HitsCol10 = (B1HitsCollection*)(HCE->GetHC(Col10));
HitsCol11 = (B1HitsCollection*)(HCE->GetHC(Col11));
HitsCol12 = (B1HitsCollection*)(HCE->GetHC(Col12));
HitsCol13 = (B1HitsCollection*)(HCE->GetHC(Col13));
HitsCol14 = (B1HitsCollection*)(HCE->GetHC(Col14));
HitsCol15 = (B1HitsCollection*)(HCE->GetHC(Col15));
HitsCol16 = (B1HitsCollection*)(HCE->GetHC(Col16));
HitsCol17 = (B1HitsCollection*)(HCE->GetHC(Col17));
HitsCol18 = (B1HitsCollection*)(HCE->GetHC(Col18));
HitsCol19 = (B1HitsCollection*)(HCE->GetHC(Col19));
}
//G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
std::ofstream gammafile("gamma.txt");
std::ofstream mu_plusfile("mu+.txt");
std::ofstream mu_minusfile("mu-.txt");
if(HitsCol1) {
G4double pos_detec = 0.01;
int n_hit = HitsCol1->entries();
G4cout << "My detector 1 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 = (*HitsCol1)[i1];
const G4String name = hit->getParticleInTarget();
//G4cout << name << G4endl;
fparticles[name]++;
WriteHistogram(name, 1);
}
PrintParticles(fparticles, gammafile, mu_plusfile, mu_minusfile, pos_detec);
fparticles.clear();
}
if(HitsCol2) {
G4double pos_detec = 1.01;
int n_hit = HitsCol2->entries();
G4cout << "My detector 2 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 = (*HitsCol2)[i1];
const G4String name = hit->getParticleInTarget();
//G4cout << name << G4endl;
fparticles[name]++;
WriteHistogram(name, 2);
}
PrintParticles(fparticles, gammafile, mu_plusfile, mu_minusfile, pos_detec);
fparticles.clear();
}