electronScattering2/src/ElectronBenchmarkDetector.cc

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 /// \file ElectronBenchmarkDetector.cc
0027 /// \brief Implementation of the ElectronBenchmarkDetector class
0028
0029 #include "ElectronBenchmarkDetector.hh"
0030
0031 #include "ElectronBenchmarkDetectorMessenger.hh"
0032
0033 #include "G4Colour.hh"
0034 #include "G4GeometryManager.hh"
0035 #include "G4LogicalVolume.hh"
0036 #include "G4LogicalVolumeStore.hh"
0037 #include "G4Material.hh"
0038 #include "G4MultiFunctionalDetector.hh"
0039 #include "G4NistManager.hh"
0040 #include "G4PSCellFlux.hh"
0041 #include "G4PSPopulation.hh"
0042 #include "G4PVPlacement.hh"
0043 #include "G4PVReplica.hh"
0044 #include "G4PhysicalVolumeStore.hh"
0045 #include "G4RunManager.hh"
0046 #include "G4SDManager.hh"
0047 #include "G4SDParticleFilter.hh"
0048 #include "G4SolidStore.hh"
0049 #include "G4SystemOfUnits.hh"
0050 #include "G4Tubs.hh"
0051 #include "G4UImanager.hh"
0052 #include "G4VPrimitiveScorer.hh"
0053 #include "G4VisAttributes.hh"
0054
0055 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0056
0057 ElectronBenchmarkDetector::ElectronBenchmarkDetector() : G4VUserDetectorConstruction()
0058 {
0059   // Exit Window
0060   fPosWindow0 = 0.000000 * cm;
0061   fPosWindow1 = 0.004120 * cm;
0062
0063   // Scattering Foil
0064   fPosPrimFoil = 2.650000 * cm;
0065   fHalfThicknessPrimFoil = 0.0 * cm;
0066
0067   // Monitor Chamber
0068   fPosMon0 = 5.000000 * cm;
0069   fPosMon1 = 5.011270 * cm;
0070
0071   // Helium Bag
0072   fPosBag0 = 6.497500 * cm;
0073   fPosHelium0 = 6.500000 * cm;
0074   fPosHelium1 = 116.500000 * cm;
0075   fPosBag1 = 116.502500 * cm;
0076   fThicknessRing = 1.4 * cm;
0077
0078   // Scoring Plane
0079   fPosScorer = 118.200000 * cm;
0080   fThicknessScorer = 0.001 * cm;
0081   fWidthScorerRing = 0.1 * cm;
0082
0083   // Radii
0084   fRadOverall = 23.3 * cm;
0085   fRadRingInner = 20.0 * cm;
0086
0087   // Extra space remaining in world volume around apparatus
0088   fPosDelta = 1. * cm;
0089   fRadDelta = 0.1 * cm;
0090
0091   fMessenger = new ElectronBenchmarkDetectorMessenger(this);
0092   DefineMaterials();
0093 }
0094
0095 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0096
0097 ElectronBenchmarkDetector::~ElectronBenchmarkDetector()
0098 {
0099   delete fMessenger;
0100
0101   delete fWorldVisAtt;
0102   delete fWindowVisAtt;
0103   delete fPrimFoilVisAtt;
0104   delete fMonVisAtt;
0105   delete fBagVisAtt;
0106   delete fHeliumVisAtt;
0107   delete fRingVisAtt;
0108   delete fScorerVisAtt;
0109 }
0110
0111 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0112
0113 G4VPhysicalVolume* ElectronBenchmarkDetector::Construct()
0114 {
0115   return CreateGeometry();
0116 }
0117
0118 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0119
0120 void ElectronBenchmarkDetector::DefineMaterials()
0121 {
0122   // Use NIST database for elements and materials whereever possible.
0123   G4NistManager* man = G4NistManager::Instance();
0124   man->SetVerbose(1);
0125
0126   // Take all elements and materials from NIST
0127   man->FindOrBuildMaterial("G4_He");
0128   man->FindOrBuildMaterial("G4_Be");
0129   man->FindOrBuildMaterial("G4_Al");
0130   man->FindOrBuildMaterial("G4_Ti");
0131   man->FindOrBuildMaterial("G4_Ta");
0132   man->FindOrBuildMaterial("G4_AIR");
0133   man->FindOrBuildMaterial("G4_MYLAR");
0134
0135   G4Element* C = man->FindOrBuildElement("C");
0136   G4Element* Cu = man->FindOrBuildElement("Cu");
0137   G4Element* Au = man->FindOrBuildElement("Au");
0138   G4Element* Ti = man->FindOrBuildElement("Ti");
0139   G4Element* Al = man->FindOrBuildElement("Al");
0140   G4Element* V = man->FindOrBuildElement("V");
0141
0142   // Define materials not in NIST.
0143   // While the NIST database does contain default materials for C, Cu and Au,
0144   // those defaults have different densities than the ones used in the
0145   // benchmark specification.
0146   G4double density;
0147   G4int ncomponents;
0148   G4double fractionmass;
0149
0150   G4Material* G4_C = new G4Material("G4_C", density = 2.18 * g / cm3, ncomponents = 1);
0151   G4_C->AddElement(C, fractionmass = 1.00);
0152
0153   G4Material* G4_Cu = new G4Material("G4_Cu", density = 8.92 * g / cm3, ncomponents = 1);
0154   G4_Cu->AddElement(Cu, fractionmass = 1.00);
0155
0156   G4Material* G4_Au = new G4Material("G4_Au", density = 19.30 * g / cm3, ncomponents = 1);
0157   G4_Au->AddElement(Au, fractionmass = 1.00);
0158
0159   G4Material* TiAlloy = new G4Material("TiAlloy", density = 4.42 * g / cm3, ncomponents = 3);
0160   TiAlloy->AddElement(Ti, fractionmass = 0.90);
0161   TiAlloy->AddElement(Al, fractionmass = 0.06);
0162   TiAlloy->AddElement(V, fractionmass = 0.04);
0163
0164   // Print materials table
0165   G4cout << *(G4Material::GetMaterialTable()) << G4endl;
0166 }
0167
0168 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0169
0170 G4VPhysicalVolume* ElectronBenchmarkDetector::CreateGeometry()
0171 {
0172   if (fPhysiWorld) return fPhysiWorld;
0173
0174   // Instantiate the world
0175   fPhysiWorld = CreateWorld();
0176   fLogWorld = fPhysiWorld->GetLogicalVolume();
0177
0178   // Instantiate the geometry
0179   CreateExitWindow(fLogWorld);
0180   CreatePrimaryFoil(fLogWorld);
0181   CreateMonitor(fLogWorld);
0182   CreateHeliumBag(fLogWorld);
0183
0184   // Create the scorers
0185   CreateScorer(fLogWorld);
0186
0187   return fPhysiWorld;
0188 }
0189
0190 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0191
0192 G4VPhysicalVolume* ElectronBenchmarkDetector::CreateWorld()
0193 {
0194   G4double halfLengthWorld = fPosScorer / 2. + fPosDelta;
0195   G4double radWorld = fRadOverall + fRadDelta;
0196   G4VSolid* worldSolid =
0197     new G4Tubs("WorldSolid", 0. * cm, radWorld, halfLengthWorld, 0. * deg, 360. * deg);
0198   G4LogicalVolume* worldLog =
0199     new G4LogicalVolume(worldSolid, G4Material::GetMaterial("G4_AIR"), "WorldLog");
0200
0201   fWorldVisAtt = new G4VisAttributes(G4Colour(1.0, 1.0, 1.0));
0202   worldLog->SetVisAttributes(fWorldVisAtt);
0203
0204   G4VPhysicalVolume* worldPhys =
0205     new G4PVPlacement(0, G4ThreeVector(0., 0., 0.), worldLog, "World", 0, false, 0);
0206
0207   return worldPhys;
0208 }
0209
0210 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0211
0212 void ElectronBenchmarkDetector::CreateExitWindow(G4LogicalVolume* worldLog)
0213 {
0214   G4double halfLengthWorld = fPosScorer / 2.;
0215   G4double halfThicknessWindow = fPosWindow1 / 2.;
0216   G4VSolid* windowSolid =
0217     new G4Tubs("windowSolid", 0. * cm, fRadOverall, halfThicknessWindow, 0. * deg, 360. * deg);
0218   G4LogicalVolume* windowLog =
0219     new G4LogicalVolume(windowSolid, G4Material::GetMaterial("TiAlloy"), "windowLog");
0220
0221   fWindowVisAtt = new G4VisAttributes(G4Colour(0.5, 1.0, 0.5));
0222   windowLog->SetVisAttributes(fWindowVisAtt);
0223
0224   new G4PVPlacement(0, G4ThreeVector(0., 0., halfThicknessWindow - halfLengthWorld), windowLog,
0225                     "ExitWindow", worldLog, false, 0);
0226 }
0227
0228 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0229
0230 void ElectronBenchmarkDetector::CreatePrimaryFoil(G4LogicalVolume* worldLog)
0231 {
0232   G4double halfLengthWorld = fPosScorer / 2.;
0233
0234   // For some energies, we have no Primary Foil.
0235   if (fHalfThicknessPrimFoil == 0.) return;
0236
0237   fSolidPrimFoil =
0238     new G4Tubs("PrimFoilSolid", 0. * cm, fRadOverall, fHalfThicknessPrimFoil, 0. * deg, 360. * deg);
0239   fLogPrimFoil = new G4LogicalVolume(fSolidPrimFoil, fMaterialPrimFoil, "PrimFoilLog");
0240
0241   fPrimFoilVisAtt = new G4VisAttributes(G4Colour(0.5, 1.0, 0.5));
0242   fLogPrimFoil->SetVisAttributes(fPrimFoilVisAtt);
0243
0244   new G4PVPlacement(0,
0245                     G4ThreeVector(0., 0., fPosPrimFoil + fHalfThicknessPrimFoil - halfLengthWorld),
0246                     fLogPrimFoil, "ScatteringFoil", worldLog, false, 0);
0247 }
0248
0249 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0250
0251 void ElectronBenchmarkDetector::CreateMonitor(G4LogicalVolume* worldLog)
0252 {
0253   G4double halfLengthWorld = fPosScorer / 2.;
0254   G4double halfThicknessMon = (fPosMon1 - fPosMon0) / 2.;
0255   G4VSolid* monSolid =
0256     new G4Tubs("monSolid", 0. * cm, fRadOverall, halfThicknessMon, 0. * deg, 360. * deg);
0257   G4LogicalVolume* monLog =
0258     new G4LogicalVolume(monSolid, G4Material::GetMaterial("G4_MYLAR"), "monLog");
0259
0260   fMonVisAtt = new G4VisAttributes(G4Colour(0.5, 1.0, 0.5));
0261   monLog->SetVisAttributes(fMonVisAtt);
0262
0263   new G4PVPlacement(0, G4ThreeVector(0., 0., fPosMon0 + halfThicknessMon - halfLengthWorld), monLog,
0264                     "MonitorChamber", worldLog, false, 0);
0265 }
0266
0267 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0268
0269 void ElectronBenchmarkDetector::CreateHeliumBag(G4LogicalVolume* worldLog)
0270 {
0271   G4double halfLengthWorld = fPosScorer / 2.;
0272
0273   // Construct cylinder of Mylar
0274   G4double halfThicknessBag = (fPosBag1 - fPosBag0) / 2.;
0275   G4VSolid* bagSolid =
0276     new G4Tubs("bagSolid", 0. * cm, fRadOverall, halfThicknessBag, 0. * deg, 360. * deg);
0277   G4LogicalVolume* bagLog =
0278     new G4LogicalVolume(bagSolid, G4Material::GetMaterial("G4_MYLAR"), "bagLog");
0279
0280   fBagVisAtt = new G4VisAttributes(G4Colour(0.5, 1.0, 0.5));
0281   bagLog->SetVisAttributes(fBagVisAtt);
0282
0283   new G4PVPlacement(0, G4ThreeVector(0., 0., fPosBag0 + halfThicknessBag - halfLengthWorld), bagLog,
0284                     "HeliumBag", worldLog, false, 0);
0285
0286   // Insert cylinder of Helium into the Cylinder of Mylar
0287   G4double halfThicknessHelium = (fPosHelium1 - fPosHelium0) / 2.;
0288   G4VSolid* heliumSolid =
0289     new G4Tubs("heliumSolid", 0. * cm, fRadOverall, halfThicknessHelium, 0. * deg, 360. * deg);
0290   G4LogicalVolume* heliumLog =
0291     new G4LogicalVolume(heliumSolid, G4Material::GetMaterial("G4_He"), "heliumLog");
0292
0293   fHeliumVisAtt = new G4VisAttributes(G4Colour(0.5, 1.0, 0.5));
0294   heliumLog->SetVisAttributes(fHeliumVisAtt);
0295
0296   new G4PVPlacement(0, G4ThreeVector(0., 0., 0.), heliumLog, "Helium", bagLog, false, 0);
0297
0298   // Insert two rings of Aluminum into the Cylinder of Helium
0299   G4double halfThicknessRing = fThicknessRing / 2.;
0300   G4VSolid* ringSolid =
0301     new G4Tubs("ringSolid", fRadRingInner, fRadOverall, halfThicknessRing, 0. * deg, 360. * deg);
0302   G4LogicalVolume* ring0Log =
0303     new G4LogicalVolume(ringSolid, G4Material::GetMaterial("G4_Al"), "ring0Log");
0304   G4LogicalVolume* ring1Log =
0305     new G4LogicalVolume(ringSolid, G4Material::GetMaterial("G4_Al"), "ring1Log");
0306
0307   fRingVisAtt = new G4VisAttributes(G4Colour(0.5, 1.0, 0.5));
0308   ring0Log->SetVisAttributes(fRingVisAtt);
0309   ring1Log->SetVisAttributes(fRingVisAtt);
0310
0311   new G4PVPlacement(0, G4ThreeVector(0., 0., -halfThicknessHelium + halfThicknessRing), ring0Log,
0312                     "Ring0", heliumLog, false, 0);
0313
0314   new G4PVPlacement(0, G4ThreeVector(0., 0., halfThicknessHelium - halfThicknessRing), ring1Log,
0315                     "Ring1", heliumLog, false, 0);
0316 }
0317
0318 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0319
0320 void ElectronBenchmarkDetector::CreateScorer(G4LogicalVolume* worldLog)
0321 {
0322   G4double halfLengthWorld = fPosScorer / 2.;
0323   G4double halfThicknessScorer = fThicknessScorer / 2.;
0324
0325   G4VSolid* scorerSolid =
0326     new G4Tubs("scorerSolid", 0. * cm, fRadOverall, halfThicknessScorer, 0. * deg, 360. * deg);
0327   G4LogicalVolume* scorerLog =
0328     new G4LogicalVolume(scorerSolid, G4Material::GetMaterial("G4_AIR"), "scorerLog");
0329
0330   fScorerVisAtt = new G4VisAttributes(G4Colour(0.5, 1.0, 0.5));
0331   scorerLog->SetVisAttributes(fScorerVisAtt);
0332   new G4PVPlacement(0, G4ThreeVector(0., 0., halfLengthWorld - halfThicknessScorer), scorerLog,
0333                     "Scorer", worldLog, false, 0);
0334
0335   G4VSolid* scorerRingSolid =
0336     new G4Tubs("scorerRingSolid", 0. * cm, fRadOverall, halfThicknessScorer, 0. * deg, 360. * deg);
0337   fScorerRingLog =
0338     new G4LogicalVolume(scorerRingSolid, G4Material::GetMaterial("G4_AIR"), "scorerRingLog");
0339   new G4PVReplica("ScorerRing", fScorerRingLog, scorerLog, kRho,
0340                   G4int(fRadOverall / fWidthScorerRing), fWidthScorerRing);
0341 }
0342
0343 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0344
0345 // Note that this method is called both at start of job and again after
0346 // any command causes a change to detector geometry
0347 void ElectronBenchmarkDetector::ConstructSDandField()
0348 {
0349   G4SDManager::GetSDMpointer()->SetVerboseLevel(1);
0350
0351   // G4Cache mechanism is necessary for multi-threaded operation
0352   // as it allows us to store separate detector pointer per thread
0353   G4MultiFunctionalDetector*& sensitiveDetector = fSensitiveDetectorCache.Get();
0354
0355   if (!sensitiveDetector) {
0356     sensitiveDetector = new G4MultiFunctionalDetector("MyDetector");
0357
0358     G4VPrimitiveScorer* primitive;
0359
0360     G4SDParticleFilter* electronFilter = new G4SDParticleFilter("electronFilter", "e-");
0361
0362     primitive = new G4PSCellFlux("cell flux");
0363     sensitiveDetector->RegisterPrimitive(primitive);
0364
0365     primitive = new G4PSCellFlux("e cell flux");
0366     primitive->SetFilter(electronFilter);
0367     sensitiveDetector->RegisterPrimitive(primitive);
0368
0369     primitive = new G4PSPopulation("population");
0370     sensitiveDetector->RegisterPrimitive(primitive);
0371
0372     primitive = new G4PSPopulation("e population");
0373     primitive->SetFilter(electronFilter);
0374     sensitiveDetector->RegisterPrimitive(primitive);
0375   }
0376   G4SDManager::GetSDMpointer()->AddNewDetector(sensitiveDetector);
0377   fScorerRingLog->SetSensitiveDetector(sensitiveDetector);
0378 }
0379
0380 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0381
0382 void ElectronBenchmarkDetector::SetPrimFoilMaterial(const G4String& matname)
0383 {
0384   G4Material* material = G4NistManager::Instance()->FindOrBuildMaterial(matname);
0385
0386   if (material && material != fMaterialPrimFoil) {
0387     fMaterialPrimFoil = material;
0388     if (fLogPrimFoil) {
0389       fLogPrimFoil->SetMaterial(fMaterialPrimFoil);
0390     }
0391     G4RunManager::GetRunManager()->PhysicsHasBeenModified();
0392   }
0393 }
0394
0395 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0396
0397 void ElectronBenchmarkDetector::SetPrimFoilThickness(G4double thicknessPrimFoil)
0398 {
0399   fHalfThicknessPrimFoil = thicknessPrimFoil / 2.;
0400 }
0401
0402 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......