ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
DetectorConstruction.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file DetectorConstruction.cc
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
28 //
29 //
31 //
32 // TestEm8: Gaseous detector
33 //
34 // Created: 31.08.2010 V.Ivanchenko ob base of V.Grichine code
35 //
36 // Modified:
37 //
39 //
40 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
41 
42 #include "DetectorConstruction.hh"
43 #include "DetectorMessenger.hh"
44 #include "TargetSD.hh"
45 
46 #include "G4Material.hh"
47 #include "G4Tubs.hh"
48 #include "G4LogicalVolume.hh"
49 #include "G4PVPlacement.hh"
50 
51 #include "G4SDManager.hh"
52 #include "G4GeometryManager.hh"
53 #include "G4RunManager.hh"
54 #include "G4NistManager.hh"
55 
56 #include "G4Region.hh"
57 #include "G4RegionStore.hh"
58 #include "G4PhysicalVolumeStore.hh"
59 #include "G4LogicalVolumeStore.hh"
60 #include "G4SolidStore.hh"
61 #include "G4ProductionCuts.hh"
62 
63 #include "G4VisAttributes.hh"
64 #include "G4Colour.hh"
65 
66 #include "G4UnitsTable.hh"
67 #include "G4PhysicalConstants.hh"
68 #include "G4SystemOfUnits.hh"
69 #include "G4ios.hh"
70 #include "TestParameters.hh"
71 #include "G4PionPlus.hh"
72 
73 #include <vector>
74 
75 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
76 
77 DetectorConstruction::DetectorConstruction()
78  : G4VUserDetectorConstruction(),
79  fGasMat(nullptr), fWindowMat(nullptr), fWorldMaterial(nullptr),
80  fSolidWorld(nullptr), fSolidContainer(nullptr), fSolidDetector(nullptr),
81  fPhysWorld(nullptr), fLogicWorld(nullptr), fLogicContainer(nullptr),
82  fLogicDetector(nullptr),fRegGasDet(nullptr)
83 {
84  fGasThickness = 23.0*mm;
85  fGasRadius = 10.*cm;
86  fMaxStep = DBL_MAX;
87 
88  fWindowThick = 51.0*micrometer;
89 
90  DefineMaterials();
91 
92  fDetectorMessenger = new DetectorMessenger(this);
93 
94  G4double cut = 0.7*mm;
95  fGasDetectorCuts = new G4ProductionCuts();
96  fGasDetectorCuts->SetProductionCut(cut,"gamma");
97  fGasDetectorCuts->SetProductionCut(cut,"e-");
98  fGasDetectorCuts->SetProductionCut(cut,"e+");
99  fGasDetectorCuts->SetProductionCut(cut,"proton");
100 }
101 
102 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
103 
104 DetectorConstruction::~DetectorConstruction()
105 {
106  delete fDetectorMessenger;
107 }
108 
109 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
110 
111 void DetectorConstruction::DefineMaterials()
112 {
113  //This function illustrates the possible ways to define materials
114  G4String name, symbol ;
115  G4double density;
116  G4int nel;
117  G4int ncomponents;
118  G4double fractionmass;
119 
120  G4NistManager* manager = G4NistManager::Instance();
121  //
122  // define Elements
123  //
124  G4Element* elH = manager->FindOrBuildElement(1);
125  G4Element* elC = manager->FindOrBuildElement(6);
126  G4Element* elO = manager->FindOrBuildElement(8);
127  G4Element* elF = manager->FindOrBuildElement(9);
128  G4Element* elNe = manager->FindOrBuildElement(10);
129  G4Element* elXe = manager->FindOrBuildElement(54);
130  //
131  // simple gases at STP conditions
132  //
133  G4Material* Argon = manager->FindOrBuildMaterial("G4_Ar");
134  G4Material* Kr = manager->FindOrBuildMaterial("G4_Kr");
135  G4Material* Xe = manager->FindOrBuildMaterial("G4_Xe");
136  //
137  // gases at STP conditions
138  //
139  G4Material* CarbonDioxide =
140  manager->FindOrBuildMaterial("G4_CARBON_DIOXIDE");
141  G4Material* Mylar = manager->FindOrBuildMaterial("G4_MYLAR");
142  G4Material* Methane= manager->FindOrBuildMaterial("G4_METHANE");
143  G4Material* Propane= manager->FindOrBuildMaterial("G4_PROPANE");
144 
145  // propane at 10 atmospheres
146  manager->ConstructNewGasMaterial("Propane10","G4_PROPANE",
147  NTP_Temperature,10.*atmosphere);
148 
149  G4Material* empty = manager->FindOrBuildMaterial("G4_Galactic");
150 
151  // 93% Kr + 7% CH4, STP
152  density = 3.491*mg/cm3 ;
153  G4Material* Kr7CH4 =
154  new G4Material(name="Kr7CH4" , density,
155  ncomponents=2);
156  Kr7CH4->AddMaterial( Kr, fractionmass = 0.986 ) ;
157  Kr7CH4->AddMaterial( Methane, fractionmass = 0.014 ) ;
158 
159  G4double TRT_Xe_density = 5.485*mg/cm3;
160  G4Material* TRT_Xe =
161  new G4Material(name="TRT_Xe", TRT_Xe_density, nel=1,
162  kStateGas,293.15*kelvin,1.*atmosphere);
163  TRT_Xe->AddElement(elXe,1);
164 
165  G4double TRT_CO2_density = 1.842*mg/cm3;
166  G4Material* TRT_CO2 =
167  new G4Material(name="TRT_CO2", TRT_CO2_density, nel=2,
168  kStateGas,293.15*kelvin,1.*atmosphere);
169  TRT_CO2->AddElement(elC,1);
170  TRT_CO2->AddElement(elO,2);
171 
172  // check alternative constructor
173  std::vector<G4String> trtatom = {"C","O"};
174  std::vector<G4int> trtnum = {1, 2};
175  manager->ConstructNewMaterial("TRT_CO2p",trtatom,trtnum,TRT_CO2_density,
176  true,kStateGas,NTP_Temperature,atmosphere);
177 
178  G4double TRT_CF4_density = 3.9*mg/cm3;
179  G4Material* TRT_CF4 =
180  new G4Material(name="TRT_CF4", TRT_CF4_density, nel=2,
181  kStateGas,293.15*kelvin,1.*atmosphere);
182  TRT_CF4->AddElement(elC,1);
183  TRT_CF4->AddElement(elF,4);
184 
185  // ATLAS TRT straw tube gas mixture (20 C, 1 atm)
186  G4double XeCO2CF4_density = 4.76*mg/cm3;
187  G4Material* XeCO2CF4 =
188  new G4Material(name="XeCO2CF4", XeCO2CF4_density,
189  ncomponents=3,
190  kStateGas,293.15*kelvin,1.*atmosphere);
191  XeCO2CF4->AddMaterial(TRT_Xe,0.807);
192  XeCO2CF4->AddMaterial(TRT_CO2,0.039);
193  XeCO2CF4->AddMaterial(TRT_CF4,0.154);
194 
195  // C3H8,20 C, 2 atm
196  density = 3.758*mg/cm3;
197  G4Material* C3H8 = new G4Material(name="C3H8",density,nel=2,
198  kStateGas,293.15*kelvin,2.*atmosphere);
199  C3H8->AddElement(elC,3);
200  C3H8->AddElement(elH,8);
201 
202  // The same material via different constructor
203  std::vector<G4String> elmname = {"C","H"};
204  std::vector<G4int> atomnum = {3, 8};
205  manager->ConstructNewIdealGasMaterial("C3H8p",elmname,atomnum,true,
206  293.15*kelvin,2.*atmosphere);
207 
208  // 87.5% Xe + 7.5% CH4 + 5% C3H8, 20 C, 1. atm
209  density = 4.9196*mg/cm3 ;
210  G4Material* XeCH4C3H8 =
211  new G4Material(name="XeCH4C3H8" ,
212  density, ncomponents=3,
213  kStateGas,NTP_Temperature,1.*atmosphere);
214  XeCH4C3H8->AddMaterial( Xe, fractionmass = 0.971);
215  XeCH4C3H8->AddMaterial( Methane, fractionmass = 0.010);
216  XeCH4C3H8->AddMaterial( Propane, fractionmass = 0.019);
217 
218  // 93% Ar + 7% CH4, STP
219  density = 1.709*mg/cm3 ;
220  G4Material* Ar7CH4 =
221  new G4Material(name="Ar7CH4", density, ncomponents=2,
222  kStateGas,STP_Temperature,STP_Pressure);
223  Ar7CH4->AddMaterial( Argon, fractionmass = 0.971 ) ;
224  Ar7CH4->AddMaterial( Methane, fractionmass = 0.029 ) ;
225 
226  // 80% Ar + 20% CO2, STP
227  density = 1.8223*mg/cm3 ;
228  G4Material* Ar_80CO2_20 =
229  new G4Material(name="ArCO2" , density, ncomponents=2,
230  kStateGas,STP_Temperature,STP_Pressure);
231  Ar_80CO2_20->AddMaterial( Argon, fractionmass = 0.783 ) ;
232  Ar_80CO2_20->AddMaterial( CarbonDioxide, fractionmass = 0.217 ) ;
233 
234  // 80% Xe + 20% CO2, STP
235  density = 5.0818*mg/cm3 ;
236  G4Material* Xe20CO2 =
237  new G4Material(name="Xe20CO2", density, ncomponents=2,
238  kStateGas,STP_Temperature,STP_Pressure);
239  Xe20CO2->AddMaterial( Xe, fractionmass = 0.922 ) ;
240  Xe20CO2->AddMaterial( CarbonDioxide, fractionmass = 0.078 ) ;
241 
242  // 80% Kr + 20% CO2, STP
243  density = 3.601*mg/cm3 ;
244  G4Material* Kr20CO2 =
245  new G4Material(name="Kr20CO2", density, ncomponents=2,
246  kStateGas,STP_Temperature,STP_Pressure);
247  Kr20CO2->AddMaterial( Kr, fractionmass = 0.89 ) ;
248  Kr20CO2->AddMaterial( CarbonDioxide, fractionmass = 0.11 ) ;
249 
250  // ALICE mixture TPC_Ne-CO2-2
251  density = 0.939*mg/cm3 ;
252  G4Material* NeCO2 =
253  new G4Material(name="TPC_Ne-CO2-2", density, ncomponents=3,
254  kStateGas,NTP_Temperature,1.*atmosphere);
255  NeCO2->AddElement( elNe, fractionmass = 0.8039 ) ;
256  NeCO2->AddElement( elO, fractionmass = 0.1426 ) ;
257  NeCO2->AddElement( elC, fractionmass = 0.0535 ) ;
258 
259  // check alternative constructor
260  std::vector<G4String> neatom = {"Ne","O","C"};
261  std::vector<G4double> nefr = {0.8039, 0.1426, 0.0536};
262  manager->ConstructNewMaterial("TPC_Ne-CO2-2p",neatom,nefr,density,true,
263  kStateGas,NTP_Temperature,atmosphere);
264 
265  // ALICE TRD mixure 85% Xe + 15% CO2 NTP
266  density = 4.9389*mg/cm3 ;
267  G4Material* Xe15CO2 =
268  new G4Material(name="Xe15CO2", density, ncomponents=2,
269  kStateGas,NTP_Temperature,atmosphere);
270  Xe15CO2->AddMaterial( Xe, fractionmass = 0.944 );
271  Xe15CO2->AddMaterial( CarbonDioxide, fractionmass = 0.056 );
272 
273  fGasMat = XeCH4C3H8;
274  fWindowMat = Mylar;
275  fWorldMaterial = empty;
276 
277  G4cout << *(G4Material::GetMaterialTable()) << G4endl;
278 
279 }
280 
281 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
282 
283 G4VPhysicalVolume* DetectorConstruction::Construct()
284 {
285  G4double contThick = fWindowThick*2 + fGasThickness;
286  G4double contR = fWindowThick*2 + fGasRadius;
287 
288  G4double worldSizeZ = contThick*1.2;
289  G4double worldSizeR = contR*1.2;
290 
291  TestParameters::GetPointer()->SetPositionZ(-0.55*contThick);
292 
293  // Printout parameters
294  G4cout << "\n The WORLD is made of "
295  << worldSizeZ/mm << "mm of " << fWorldMaterial->GetName() ;
296  G4cout << ", the transverse size (R) of the world is " << worldSizeR/mm
297  << " mm. " << G4endl;
298  G4cout << " The CONTAINER is made of "
299  << fWindowThick/mm << "mm of " << fWindowMat->GetName() << G4endl;
300  G4cout << " The TARGET is made of "
301  << fGasThickness/mm << "mm of " << fGasMat->GetName() ;
302  G4cout << ", the transverse size (R) is " << fGasRadius/mm << " mm. "
303  << G4endl;
304  G4cout << G4endl;
305 
306  // World
307  fSolidWorld =
308  new G4Tubs("World",0.,worldSizeR,worldSizeZ/2.,0.,CLHEP::twopi);
309 
310  fLogicWorld = new G4LogicalVolume(fSolidWorld, fWorldMaterial, "World");
311 
312  fPhysWorld = new G4PVPlacement(0,
313  G4ThreeVector(0.,0.,0.),
314  "World",
315  fLogicWorld,
316  0,
317  false,
318  0);
319 
320  // Window
321  fSolidContainer = new G4Tubs("Absorber",
322  0.,contR,contThick/2.,0.,CLHEP::twopi);
323 
324  fLogicContainer = new G4LogicalVolume(fSolidContainer, fWindowMat, "Window");
325 
326  G4PVPlacement* PhysWind = new G4PVPlacement(0, G4ThreeVector(0.,0.,0.),
327  "Window", fLogicContainer,
328  fPhysWorld, false, 0);
329 
330  // Detector volume
331  fSolidDetector = new G4Tubs("Gas", 0., fGasRadius, fGasThickness/2.,
332  0., CLHEP::twopi);
333 
334  fLogicDetector = new G4LogicalVolume(fSolidDetector, fGasMat, "Gas");
335 
336  new G4PVPlacement(0, G4ThreeVector(0.,0.,0.), "Gas", fLogicDetector,
337  PhysWind, false, 0);
338 
339  // defined gas detector region
340  fRegGasDet = new G4Region("GasDetector");
341  fRegGasDet->SetProductionCuts(fGasDetectorCuts);
342  fRegGasDet->AddRootLogicalVolume(fLogicDetector);
343 
344  // visualisation
345  fLogicWorld->SetVisAttributes(G4VisAttributes::GetInvisible());
346  G4VisAttributes* color1 = new G4VisAttributes(G4Colour(0.3, 0.3, 0.3));
347  fLogicContainer->SetVisAttributes(color1);
348  G4VisAttributes* color2 = new G4VisAttributes(G4Colour(0.0, 0.3, 0.7));
349  fLogicDetector->SetVisAttributes(color2);
350 
351  if(0.0 == fGasMat->GetIonisation()->GetMeanEnergyPerIonPair()) {
352  SetPairEnergy(20*eV);
353  }
354  return fPhysWorld;
355 }
356 
357 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
358 
359 void DetectorConstruction::ConstructSDandField()
360 {
361  auto sd = new TargetSD("GasSD");
362  G4SDManager::GetSDMpointer()->AddNewDetector(sd);
363  SetSensitiveDetector(fLogicDetector, sd);
364 }
365 
366 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
367 
368 void DetectorConstruction::SetGasMaterial(const G4String& name)
369 {
370  // get the pointer to the existing material
371  G4Material* mat = G4Material::GetMaterial(name, false);
372 
373  // create the material by its name
374  if(!mat) { mat = G4NistManager::Instance()->FindOrBuildMaterial(name); }
375 
376  if (mat && mat != fGasMat) {
377  G4cout << "### New target material: " << mat->GetName() << G4endl;
378  fGasMat = mat;
379  if(fLogicDetector) {
380  fLogicDetector->SetMaterial(mat);
381  G4RunManager::GetRunManager()->PhysicsHasBeenModified();
382  }
383  }
384 }
385 
386 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
387 
388 void DetectorConstruction::SetContainerMaterial(const G4String& name)
389 {
390  // get the pointer to the existing material
391  G4Material* mat = G4Material::GetMaterial(name, false);
392 
393  // create the material by its name
394  if(!mat) { mat = G4NistManager::Instance()->FindOrBuildMaterial(name); }
395 
396  if (mat && mat != fWindowMat) {
397  G4cout << "### New material for container: " << mat->GetName() << G4endl;
398  fWindowMat = mat;
399  if(fLogicContainer) {
400  fLogicContainer->SetMaterial(mat);
401  G4RunManager::GetRunManager()->PhysicsHasBeenModified();
402  }
403  }
404 }
405 
406 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
407 
408 void DetectorConstruction::SetWorldMaterial(const G4String& name)
409 {
410  // get the pointer to the existing material
411  G4Material* mat = G4Material::GetMaterial(name, false);
412 
413  // create the material by its name
414  if(!mat) { mat = G4NistManager::Instance()->FindOrBuildMaterial(name); }
415 
416  if (mat && mat != fWorldMaterial) {
417  G4cout << "### New World material: " << mat->GetName() << G4endl;
418  fWorldMaterial = mat;
419  if(fLogicWorld) {
420  fLogicWorld->SetMaterial(mat);
421  G4RunManager::GetRunManager()->PhysicsHasBeenModified();
422  }
423  }
424 }
425 
426 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
427 
428 void DetectorConstruction::SetGasThickness(G4double val)
429 {
430  fGasThickness = val;
431  if(fPhysWorld) { ChangeGeometry(); }
432 }
433 
434 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
435 
436 void DetectorConstruction::SetGasRadius(G4double val)
437 {
438  fGasRadius = val;
439  if(fPhysWorld) { ChangeGeometry(); }
440 }
441 
442 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
443 
444 void DetectorConstruction::SetContainerThickness(G4double val)
445 {
446  fWindowThick = val;
447  if(fPhysWorld) { ChangeGeometry(); }
448 }
449 
450 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
451 
452 void DetectorConstruction::SetPairEnergy(G4double val)
453 {
454  if(val > 0.0) {
455  fGasMat->GetIonisation()->SetMeanEnergyPerIonPair(val);
456  }
457 }
458 
459 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
460 
461 void DetectorConstruction::ChangeGeometry()
462 {
463  G4double contThick = fWindowThick*2 + fGasThickness;
464  G4double contR = fWindowThick*2 + fGasRadius;
465 
466  G4double worldSizeZ = contThick*1.2;
467  G4double worldSizeR = contR*1.2;
468 
469  TestParameters::GetPointer()->SetPositionZ(-0.55*contThick);
470 
471  fSolidWorld->SetOuterRadius(worldSizeR);
472  fSolidWorld->SetZHalfLength(worldSizeZ*0.5);
473 
474  fSolidContainer->SetOuterRadius(contR);
475  fSolidContainer->SetZHalfLength(contThick*0.5);
476 
477  fSolidDetector->SetOuterRadius(fGasRadius);
478  fSolidDetector->SetZHalfLength(fGasThickness*0.5);
479 }
480 
481 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......