diff options
Diffstat (limited to 'src')
| -rw-r--r-- | src/#G4BeamTestSiHit.cxx# | 115 | ||||
| -rw-r--r-- | src/G4BeamTestDetectorConstruction.cxx | 90 | ||||
| -rw-r--r-- | src/G4BeamTestEventAction.cxx | 24 | ||||
| -rw-r--r-- | src/G4BeamTestSiHit.cxx | 39 | ||||
| -rw-r--r-- | src/G4BeamTestSiSD.cxx | 17 | ||||
| -rw-r--r-- | src/G4BeamTestTank.cxx | 9 | ||||
| -rw-r--r-- | src/G4BeamTestUserTrackingAction.cxx | 14 |
7 files changed, 265 insertions, 43 deletions
diff --git a/src/#G4BeamTestSiHit.cxx# b/src/#G4BeamTestSiHit.cxx# new file mode 100644 index 0000000..663a03f --- /dev/null +++ b/src/#G4BeamTestSiHit.cxx# @@ -0,0 +1,115 @@ +#include <iostream> +#include <iomanip> + +#include "G4BeamTestSiHit.h" +#include "G4UnitsTable.hh" +#include "G4VVisManager.hh" +#include "G4Circle.hh" +#include "G4Colour.hh" +#include "G4VisAttributes.hh" + +G4ThreadLocal G4Allocator<G4BeamTestSiHit>* G4BeamTestSiHitAllocator=0; +static std::fstream testnew("./testnew.txt", std::ofstream::out); + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +G4BeamTestSiHit::G4BeamTestSiHit() + : G4VHit(), + fTrackID(-1), + + fEdep(0.), + fTime(0.), + fPos(G4ThreeVector()) +{ + /* G4cout << "opening file" << G4endl; */ + /* testnew.open(testnew_out); */ +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +G4BeamTestSiHit ::~G4BeamTestSiHit() { + /* G4cout << "closing file" << G4endl; */ + /* testnew.close(); */ +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +G4BeamTestSiHit::G4BeamTestSiHit(const G4BeamTestSiHit& right) + : G4VHit() +{ + fTrackID = right.fTrackID; + // fChamberNb = right.fChamberNb; + fEdep = right.fEdep; + fPos = right.fPos; + fTime = right.fTime; +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +const G4BeamTestSiHit& G4BeamTestSiHit::operator=(const G4BeamTestSiHit& right) +{ + fTrackID = right.fTrackID; + // fChamberNb = right.fChamberNb; + fEdep = right.fEdep; + fPos = right.fPos; + fTime = right.fTime; + + return *this; +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +G4int G4BeamTestSiHit::operator==(const G4BeamTestSiHit& right) const +{ + return ( this == &right ) ? 1 : 0; +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +void G4BeamTestSiHit::Draw() +{ + /* G4VVisManager* pVVisManager = G4VVisManager::GetConcreteInstance(); */ + /* if(pVVisManager) */ + /* { */ + /* G4Circle circle(fPos); */ + /* circle.SetScreenSize(4.); */ + /* circle.SetFillStyle(G4Circle::filled); */ + /* G4Colour colour(1.,0.,0.); */ + /* G4VisAttributes attribs(colour); */ + /* circle.SetVisAttributes(attribs); */ + /* pVVisManager->Draw(circle); */ + /* } */ +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + + +void G4BeamTestSiHit::Print() +{ + G4cout + << " trackID: " << fTrackID + << "Edep: " + << std::setw(7) << G4BestUnit(fEdep,"Energy") + << " Position: " + << std::setw(7) << G4BestUnit(fPos ,"Length") + << "Time: " + << std::setw(7) << G4BestUnit(fTime,"Time") + << G4endl; +} + +void G4BeamTestSiHit::Dataout() +{ +testnew + << " trackID: " << fTrackID + << "Edep: " + << std::setw(7) << G4BestUnit(fEdep,"Energy") + << " Position: " + << std::setw(7) << G4BestUnit( fPos ,"Length") + << "Time: " + << std::setw(7) << G4BestUnit( fTime,"Time") + << G4endl; + + } + + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... diff --git a/src/G4BeamTestDetectorConstruction.cxx b/src/G4BeamTestDetectorConstruction.cxx index ebbd9b8..495330a 100644 --- a/src/G4BeamTestDetectorConstruction.cxx +++ b/src/G4BeamTestDetectorConstruction.cxx @@ -3,6 +3,7 @@ #include <G4LogicalVolume.hh> #include <G4PVPlacement.hh> +#include <G4SDManager.hh> #include <G4Box.hh> @@ -34,12 +35,17 @@ G4VPhysicalVolume* G4BeamTestDetectorConstruction::Construct() /* origin_.set(delaunay.GetOrigin().x(), delaunay.GetOrigin().y(), zSnowBottom + zHalfLength); */ // Determine World dimensions - G4double xWorld = 2.0 * CLHEP::m; - G4double yWorld = 2.0 * CLHEP::m; - G4double zWorld = 2.0 * CLHEP::m; + G4double xWorld = 4.0 * CLHEP::m; + G4double yWorld = 4.0 * CLHEP::m; + G4double zWorld = 4.0 * CLHEP::m; + + // SC4 dimensions + G4double scinHeight_ = 1 * 2.54 * CLHEP::cm; + G4double scinWidth_ = 1 * 2.54 * CLHEP::cm; + G4double scinThickness_ = 1 * 2.54 * CLHEP::cm; // Create world volume - G4Box* world_box = new G4Box("solid_world", xWorld, yWorld, zWorld); + G4Box* world_box = new G4Box("solid_world", xWorld*0.5, yWorld*0.5, zWorld*0.5); G4LogicalVolume* worldLog = new G4LogicalVolume(world_box, G4Material::GetMaterial("Air"), "log_world", 0, 0, 0); G4VPhysicalVolume* worldPhys = @@ -56,15 +62,29 @@ G4VPhysicalVolume* G4BeamTestDetectorConstruction::Construct() tank_->InstallTank(worldPhys, origin_); // } + // Define SC4 + G4Box* sc4_box = new G4Box("sc4",scinHeight_*0.5, scinWidth_*0.5, scinThickness_*0.5); + G4LogicalVolume* sc4Log = + new G4LogicalVolume(sc4_box, G4Material::GetMaterial("SC4"), "log_sc4", 0, 0, 0); + G4VPhysicalVolume* sc4Phys = + new G4PVPlacement(0, G4ThreeVector(1.2*CLHEP::m,0,0), sc4Log, "sc4", worldLog, false, 0); + // User limits (energy cutoffs) // Do not create photons or electrons below cherenkov threshold // See also corresponding UserSpecialCuts in Physicslist !!!! G4UserLimits* energyLimit = new G4UserLimits(); /* energyLimit->SetUserMinEkine(280.0 * CLHEP::keV); // Cherenkov threshold of electrons in ice */ - energyLimit->SetUserMinEkine(1.907 * CLHEP::eV); // Lower threshold of PMT - 600nm + energyLimit->SetUserMinEkine(2.26 * CLHEP::eV); // Lower threshold of PMT - 550nm + // energyLimit->SetUserMaxEkine(3.55 * CLHEP::eV); //upper threshold of PMT - 350nm worldLog->SetUserLimits(energyLimit); /* snowLog->SetUserLimits(energyLimit); */ + sc4Log->SetUserLimits(energyLimit); + // G4SDManager* sdManager = G4SDManager::GetSDMpointer(); + // sc4SD_ = new G4BeamTestSC4SD("sc4_SD_", "HitsCollection"); + // sdManager->AddNewDetector(sc4SD_); + // sc4Log->SetSensitiveDetector(sc4SD_); + return worldPhys; } @@ -80,6 +100,7 @@ void G4BeamTestDetectorConstruction::CreateMaterials() /* CreatePerlite(); */ CreateGlassSphere(); CreateEffectiveDOMMaterial(); + CreateSC4(); //if(verboseLevel_>0) G4cout << *G4Material::GetMaterialTable() << G4endl; } @@ -105,7 +126,8 @@ void G4BeamTestDetectorConstruction::CreateWater() // https://docushare.icecube.wisc.edu/dsweb/Get/Document-6637/R7081-02%20data%20sheet.pdf<Paste> // TODO(shivesh): add more properties? const G4int water_bins = 2; - G4double water_ephot[water_bins] = {1.91 * CLHEP::eV, 4.13 * CLHEP::eV}; + // G4double water_ephot[water_bins] = {1.91 * CLHEP::eV, 4.13 * CLHEP::eV}; + G4double water_ephot[water_bins] = {0.1 * CLHEP::eV, 10 * CLHEP::eV}; G4double water_refr[water_bins] = {1.33, 1.33}; G4MaterialPropertiesTable *mpt_water = new G4MaterialPropertiesTable (); mpt_water->AddProperty ("RINDEX", water_ephot, water_refr, water_bins); @@ -144,9 +166,40 @@ void G4BeamTestDetectorConstruction::CreateGlassSphere() // 20 lbs. weight = 9072 g // 6.5" outer radius & 0.5" thickness = 4024 cm3 G4NistManager* nistManager = G4NistManager::Instance(); - G4Material* glass = new G4Material("Glass", 2.254 * CLHEP::g / CLHEP::cm3, 2, kStateSolid); - glass->AddElement(nistManager->FindOrBuildElement("Si"), 1); - glass->AddElement(nistManager->FindOrBuildElement("O"), 2); + // G4Material* glass = new G4Material("Glass", 2.254 * CLHEP::g / CLHEP::cm3, 2, kStateSolid); + // glass->AddElement(nistManager->FindOrBuildElement("Si"), 1); + // glass->AddElement(nistManager->FindOrBuildElement("O"), 2); + + // http://hypernews.slac.stanford.edu/HyperNews/geant4/get/AUX/2013/03/11/12.39-85121-chDetectorConstruction.cc + // Define elements for all materials not found in the NIST database + G4Element* Si = nistManager->FindOrBuildElement("Si"); + G4Element* B = nistManager->FindOrBuildElement("B"); + G4Element* O = nistManager->FindOrBuildElement("O"); + G4Element* Na = nistManager->FindOrBuildElement("Na"); + G4Element* Al = nistManager->FindOrBuildElement("Al"); + G4Element* K = nistManager->FindOrBuildElement("K"); + + G4double density; + G4int ncomponents; + G4double fractionmass; + G4Material* glass = new G4Material("Glass", density= 2.23*CLHEP::g/CLHEP::cm3, ncomponents=6); + glass->AddElement(B, fractionmass=0.040064); + glass->AddElement(O, fractionmass=0.539562); + glass->AddElement(Na, fractionmass=0.028191); + glass->AddElement(Al, fractionmass=0.011644); + glass->AddElement(Si, fractionmass=0.377220); + glass->AddElement(K, fractionmass=0.003321); + + // pmt spectral response 300-650nm + // https://docushare.icecube.wisc.edu/dsweb/Get/Document-6637/R7081-02%20data%20sheet.pdf<Paste> + // TODO(shivesh): add more properties? + const G4int glass_bins = 2; + // G4double glass_ephot[glass_bins] = {1.91 * CLHEP::eV, 4.13 * CLHEP::eV}; + G4double glass_ephot[glass_bins] = {0.1 * CLHEP::eV, 10 * CLHEP::eV}; + G4double glass_refr[glass_bins] = {1.47, 1.47}; + G4MaterialPropertiesTable *mpt_glass = new G4MaterialPropertiesTable (); + mpt_glass->AddProperty ("RINDEX", glass_ephot, glass_refr, glass_bins); + glass->SetMaterialPropertiesTable(mpt_glass); } /*****************************************************************/ @@ -160,4 +213,23 @@ void G4BeamTestDetectorConstruction::CreateEffectiveDOMMaterial() G4Material* glass = new G4Material("effectiveDOM", 0.2 * CLHEP::g / CLHEP::cm3, 2, kStateSolid); glass->AddElement(nistManager->FindOrBuildElement("Si"), 1); glass->AddElement(nistManager->FindOrBuildElement("O"), 2); + + // const G4int glass_bins = 2; + // G4double glass_ephot[glass_bins] = {1.91 * CLHEP::eV, 4.13 * CLHEP::eV}; + // glass + // G4double glass_ephot[glass_bins] = {0.1 * CLHEP::eV, 10 * CLHEP::eV}; + // G4double glass_refr[glass_bins] = {1.47, 1.47}; + // G4MaterialPropertiesTable *mpt_glass = new G4MaterialPropertiesTable (); + // mpt_glass->AddProperty ("RINDEX", glass_ephot, glass_refr, glass_bins); + // glass->SetMaterialPropertiesTable(mpt_glass); +} + +void G4BeamTestDetectorConstruction::CreateSC4() +{ + G4NistManager* nistManager = G4NistManager::Instance(); + // POM + G4Material* plastic = new G4Material("SC4", 1.425 * CLHEP::g / CLHEP::cm3, 3, kStateSolid); + plastic->AddElement(nistManager->FindOrBuildElement("H"), 2); + plastic->AddElement(nistManager->FindOrBuildElement("C"), 1); + plastic->AddElement(nistManager->FindOrBuildElement("O"), 1); } diff --git a/src/G4BeamTestEventAction.cxx b/src/G4BeamTestEventAction.cxx index 176991d..53402c9 100644 --- a/src/G4BeamTestEventAction.cxx +++ b/src/G4BeamTestEventAction.cxx @@ -33,13 +33,14 @@ void G4BeamTestEventAction::BeginOfEventAction(const G4Event* event ) G4SDManager * SDman = G4SDManager::GetSDMpointer(); + SDman->ListTree(); if(SiCollID<0) { G4String colNam; SiCollID = SDman->GetCollectionID(colNam="G4BeamTestSiSDCollection"); } - + } @@ -49,6 +50,9 @@ void G4BeamTestEventAction::EndOfEventAction(const G4Event* event) { G4cout << ">>> Summary of Event " << event->GetEventID() << G4endl; + // testnew << ">>> Summary of Event " << event->GetEventID() << G4endl; + + G4cout << SiCollID << G4endl; if(SiCollID<0) return; @@ -63,17 +67,23 @@ void G4BeamTestEventAction::EndOfEventAction(const G4Event* event) - if(SiHC) { int n_hit = SiHC->entries(); + // testnew << std::flush; G4cout << G4endl; - G4cout << "Si hits " << - "--------------------------------------------------------------" - << G4endl; + // G4cout << "Si hits " << + // "--------------------------------------------------------------" + // << G4endl; G4cout << n_hit << " hits are stored in G4BeamTestSiHitsCollection." << G4endl; /* G4cout << "List of hits in tracker" << G4endl; */ + // testnew << G4endl; + // testnew << "Si hits " << + // "--------------------------------------------------------------" + testnew << n_hit << " hits are stored in G4BeamTestSiHitsCollection." + << G4endl; + /* testnew << "List of hits in tracker" << G4endl; */ for(int i=0;i<n_hit;i++) { /* (*SiHC)[i]->Print(); */ @@ -81,7 +91,9 @@ void G4BeamTestEventAction::EndOfEventAction(const G4Event* event) } // G4cout << "sid + " << SiCollID << G4endl; - + // testnew << "sid + " << SiCollID << G4endl; + testnew << std::flush; + } } diff --git a/src/G4BeamTestSiHit.cxx b/src/G4BeamTestSiHit.cxx index 591274c..0c90623 100644 --- a/src/G4BeamTestSiHit.cxx +++ b/src/G4BeamTestSiHit.cxx @@ -9,7 +9,6 @@ #include "G4VisAttributes.hh" G4ThreadLocal G4Allocator<G4BeamTestSiHit>* G4BeamTestSiHitAllocator=0; -static std::fstream testnew("./testnew.txt", std::ofstream::out); //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... @@ -86,28 +85,30 @@ void G4BeamTestSiHit::Draw() void G4BeamTestSiHit::Print() { - G4cout - << " trackID: " << fTrackID - << "Edep: " - << std::setw(7) << G4BestUnit(fEdep,"Energy") - << " Position: " - << std::setw(7) << G4BestUnit(fPos ,"Length") - << "Time: " - << std::setw(7) << G4BestUnit(fTime,"Time") - << G4endl; + // G4cout + // << " trackID: " << fTrackID + // << "Edep: " + // << std::setw(7) << G4BestUnit(fEdep,"Energy") + // << " Position: " + // << std::setw(7) << G4BestUnit(fPos ,"Length") + // << "Time: " + // << std::setw(7) << G4BestUnit(fTime,"Time") + // << G4endl; } void G4BeamTestSiHit::Dataout() { -testnew - << " trackID: " << fTrackID - << "Edep: " - << std::setw(7) << G4BestUnit(fEdep,"Energy") - << " Position: " - << std::setw(7) << G4BestUnit( fPos ,"Length") - << "Time: " - << std::setw(7) << G4BestUnit( fTime,"Time") - << G4endl; +// testnew << std::flush; +// testnew +// << " trackID: " << fTrackID +// << "Edep: " +// << std::setw(7) << G4BestUnit(fEdep,"Energy") +// << " Position: " +// << std::setw(7) << G4BestUnit( fPos ,"Length") +// << "Time: " +// << std::setw(7) << G4BestUnit( fTime,"Time") +// << G4endl; +// testnew << std::flush; } diff --git a/src/G4BeamTestSiSD.cxx b/src/G4BeamTestSiSD.cxx index b67092b..afa6bf8 100644 --- a/src/G4BeamTestSiSD.cxx +++ b/src/G4BeamTestSiSD.cxx @@ -45,6 +45,7 @@ void G4BeamTestSiSD::Initialize(G4HCofThisEvent* hce) G4int hcID = G4SDManager::GetSDMpointer()->GetCollectionID(collectionName[0]); + G4cout << "hcID " << hcID << G4endl; hce->AddHitsCollection( hcID, fHitsCollection ); } @@ -61,11 +62,25 @@ G4bool G4BeamTestSiSD::ProcessHits(G4Step* aStep, /* G4cout << " Particle_name = " << name << G4endl; */ if (name == "opticalphoton" || name == "gamma") { -/* G4cout << " Particle_name = " << name << G4endl; */ + // G4cout << " Particle_name = " << name << G4endl; +// total energy + G4double etot = aStep->GetTrack()->GetTotalEnergy(); // energy deposit G4double edep = aStep->GetTotalEnergyDeposit(); + if (etot < 2.26 * CLHEP::eV) { // Lower threshold of PMT - 550nm + // if (etot < 2.48 * CLHEP::eV) { // Lower threshold of PMT - 500nm + // G4cout << "particle " << name << " under threshold with energy " << etot << G4endl; + return true; + } + if (etot > 3.55 * CLHEP::eV) { // Upper threshold of PMT - 350nm + // if (etot > 3.10 * CLHEP::eV) { // Upper threshold of PMT - 400nm + // G4cout << "particle " << name << " over threshold with energy " << etot << G4endl; + return true; + } + // G4cout << "inserting particle " << name << " with energy " << etot << " into record" << G4endl; + // if (edep==0.) return false; /* G4cout << " Particle_name_after_edep = " << name << G4endl; */ diff --git a/src/G4BeamTestTank.cxx b/src/G4BeamTestTank.cxx index 56570a4..02f6751 100644 --- a/src/G4BeamTestTank.cxx +++ b/src/G4BeamTestTank.cxx @@ -28,7 +28,8 @@ G4BeamTestTank::G4BeamTestTank() { // Get tank dimensions - tankThickness_ = 0.0*CLHEP::cm; // TODO(shivesh) : check thickness + // tankThickness_ = 0.0*CLHEP::cm; // TODO(shivesh) : check thickness + tankThickness_ = 0.44 * 2.54 *CLHEP::cm; // TODO(shivesh) : check thickness tankHeight_ = 76.83 * 2.54 * CLHEP::cm; innerRadius_ = 32 * 2.54 * CLHEP::cm; outerRadius_ = innerRadius_ + tankThickness_; @@ -56,7 +57,7 @@ G4VPhysicalVolume* G4BeamTestTank::InstallTank(G4VPhysicalVolume* mother, const // See also corresponding UserSpecialCuts in Physicslist !!!! // TODO(shivesh): Maybe do all of this as stepping action ?????? G4UserLimits* energyLimit = new G4UserLimits(); - energyLimit->SetUserMinEkine(1.907 * CLHEP::eV); // Lower threshold of PMT - 600nm + energyLimit->SetUserMinEkine(2.26 * CLHEP::eV); // Lower threshold of PMT - 550nm // std::string tankName=boost::lexical_cast<std::string>(tankKey_); std::string tankName = "BTT"; @@ -140,7 +141,7 @@ G4VPhysicalVolume* G4BeamTestTank::InstallTank(G4VPhysicalVolume* mother, const new G4LogicalVolume(upperglasssphere, glass, ("log_dom_up_" + omName).c_str(), 0, 0, 0); G4LogicalVolume* logLowerGlass = - new G4LogicalVolume(lowerglasssphere, water, + new G4LogicalVolume(lowerglasssphere, glass, ("log_dom_lo_" + omName).c_str(), 0, 0, 0); G4LogicalVolume* logUpperDOM = new G4LogicalVolume(upperdomsphere, effectiveDOM, @@ -172,7 +173,7 @@ G4VPhysicalVolume* G4BeamTestTank::InstallTank(G4VPhysicalVolume* mother, const G4SDManager* sdManager = G4SDManager::GetSDMpointer(); iceSD_ = new G4BeamTestSiSD(("ice_SD_" + tankName).c_str(), "HitsCollection"); sdManager->AddNewDetector(iceSD_); - /* logLowerDOM->SetSensitiveDetector(iceSD_); */ + // logLowerDOM->SetSensitiveDetector(iceSD_); logLowerGlass->SetSensitiveDetector(iceSD_); // Instantiation of a set of visualization attributes with red colour diff --git a/src/G4BeamTestUserTrackingAction.cxx b/src/G4BeamTestUserTrackingAction.cxx index f7c053b..e41ccc3 100644 --- a/src/G4BeamTestUserTrackingAction.cxx +++ b/src/G4BeamTestUserTrackingAction.cxx @@ -15,6 +15,7 @@ void G4BeamTestUserTrackingAction::PostUserTrackingAction(const G4Track* track) G4UserLimits *limit = volume->GetUserLimits(); if(!limit) G4cout << "----> G4LogicalVolume: " << volume->GetName() << " has no defined G4UserLimit" << G4endl; G4double threshold = limit->GetUserMinEkine(*track); + G4double max_threshold = 3.54; G4TrackVector* secondaries = fpTrackingManager->GimmeSecondaries(); if(secondaries) { @@ -29,10 +30,15 @@ void G4BeamTestUserTrackingAction::PostUserTrackingAction(const G4Track* track) { //check if particle energy is below threshold; if true, kill the particle G4double energy = (*secondaries)[i]->GetTotalEnergy(); - if(energy < threshold){ - G4cout << "TrackingAction: killing particle " << particle << " with energy " << energy << " < " << threshold << G4endl; - (*secondaries)[i]->SetTrackStatus(fStopAndKill); - } + // if(energy < threshold){ + // G4cout << "TrackingAction: killing particle " << particle << " with energy " << energy << " < " << threshold << G4endl; + // (*secondaries)[i]->SetTrackStatus(fStopAndKill); + // } + // if (energy > max_threshold * CLHEP::eV){ + // G4cout << "TrackingAction: killing particle " << particle << " with energy " << energy << " > " << max_threshold << G4endl; + // (*secondaries)[i]->SetTrackStatus(fStopAndKill); + + // } } } } |