source: trunk/examples/extended/electromagnetic/TestEm8/src/Em8DetectorConstruction.cc @ 1309

Last change on this file since 1309 was 1230, checked in by garnier, 14 years ago

update to geant4.9.3
File size: 17.2 KB
Line 
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//
26//
27// $Id: Em8DetectorConstruction.cc,v 1.21 2007/10/02 10:12:47 vnivanch Exp $
28// GEANT4 tag $Name: geant4-09-03-cand-01 $
29//
30//
31
32#include "Em8DetectorConstruction.hh"
33#include "Em8DetectorMessenger.hh"
34#include "Em8CalorimeterSD.hh"
35
36
37#include "G4Material.hh"
38#include "G4Tubs.hh"
39#include "G4LogicalVolume.hh"
40#include "G4PVPlacement.hh"
41
42#include "G4FieldManager.hh"
43#include "G4TransportationManager.hh"
44#include "G4SDManager.hh"
45#include "G4GeometryManager.hh"
46#include "G4RunManager.hh"
47
48#include "G4Region.hh"
49#include "G4RegionStore.hh"
50#include "G4PhysicalVolumeStore.hh"
51#include "G4LogicalVolumeStore.hh"
52#include "G4SolidStore.hh"
53#include "G4ProductionCuts.hh"
54
55#include "G4VisAttributes.hh"
56#include "G4Colour.hh"
57
58#include "G4UnitsTable.hh"
59#include "G4ios.hh"
60
61/////////////////////////////////////////////////////////////////////////////
62//
63//
64
65Em8DetectorConstruction::Em8DetectorConstruction()
66:
67fWorldChanged(false),
68fWorldMaterial(NULL),fSolidWorld(NULL),fLogicWorld(NULL),fPhysicsWorld(NULL),
69fAbsorberMaterial(NULL),fSolidAbsorber(NULL),fLogicAbsorber(NULL),
70fPhysicsAbsorber(NULL),fDetectorMessenger(NULL),
71fCalorimeterSD(NULL),fRegGasDet(NULL)
72{
73  fDelta = 0.0001*mm;
74
75  fAbsorberThickness = 23.0*mm;
76
77  fAbsorberRadius    = 10.*cm;
78  fAbsorberZ         = 0.*cm ;
79
80  fWindowThick       = 51.0*micrometer ;
81
82  fGammaCut    = 23*mm; 
83  fElectronCut = 23*mm; 
84  fPositronCut = 23*mm; 
85
86  fDetectorMessenger = new Em8DetectorMessenger(this);
87}
88
89//////////////////////////////////////////////////////////////////////////
90//
91//
92
93Em8DetectorConstruction::~Em8DetectorConstruction()
94{ 
95  delete fDetectorMessenger;
96}
97
98//////////////////////////////////////////////////////////////////////////
99//
100//
101
102G4VPhysicalVolume* Em8DetectorConstruction::Construct()
103{
104  DefineMaterials();
105  return ConstructCalorimeter();
106}
107
108//////////////////////////////////////////////////////////////////////////////
109//
110//
111
112void Em8DetectorConstruction::DefineMaterials()
113{ 
114 //This function illustrates the possible ways to define materials
115 
116  G4String name, symbol ;             //a=mass of a mole;
117  G4double a, z, density ;            //z=mean number of protons; 
118
119  G4int nel ;                       //iz=number of protons  in an isotope;
120                                   // n=number of nucleons in an isotope;
121
122  G4int ncomponents; 
123  G4double fractionmass;
124
125  //
126  // define Elements
127  //
128
129  a = 1.01*g/mole;
130  G4Element* elH  = new G4Element(name="Hydrogen",symbol="H" , z= 1., a);
131
132  a = 12.01*g/mole;
133  G4Element* elC = new G4Element(name="Carbon", symbol="C", z=6., a);
134
135  a = 14.01*g/mole;
136  G4Element* elN  = new G4Element(name="Nitrogen",symbol="N" , z= 7., a);
137
138  a = 16.00*g/mole;
139  G4Element* elO  = new G4Element(name="Oxygen"  ,symbol="O" , z= 8., a);
140
141  a = 39.948*g/mole;
142  G4Element* elAr = new G4Element(name="Argon", symbol="Ar", z=18., a);
143
144  //
145  // define simple materials
146  //
147
148  // Aluminium
149  a = 26.98*g/mole;
150  density = 2.7*g/cm3;
151  G4Material* Al = new G4Material(name="Aluminium", z=13., a, density);
152  if(Al);
153
154  // Mylar
155  density = 1.39*g/cm3;
156  G4Material* Mylar = new G4Material(name="Mylar", density, nel=3);
157  Mylar->AddElement(elO,2);
158  Mylar->AddElement(elC,5);
159  Mylar->AddElement(elH,4);
160
161  // Silicon as detector material
162  density = 2.330*g/cm3;
163  a = 28.09*g/mole;
164  G4Material* Si = new G4Material(name="Silicon", z=14., a, density);
165  if(Si);
166
167  // Krypton as detector gas, STP
168  density = 3.700*mg/cm3 ;
169  a = 83.80*g/mole ;
170  G4Material* Kr  = new G4Material(name="Kr",z=36., a, density );
171
172
173  // Dry air (average composition)
174  density = 1.7836*mg/cm3 ;       // STP
175  G4Material* Argon = new G4Material(name="Argon"  , density, ncomponents=1);
176  Argon->AddElement(elAr, 1);
177
178  density = 1.25053*mg/cm3 ;       // STP
179  G4Material* Nitrogen = new G4Material(name="N2"  , density, ncomponents=1);
180  Nitrogen->AddElement(elN, 2);
181
182  density = 1.4289*mg/cm3 ;       // STP
183  G4Material* Oxygen = new G4Material(name="O2"  , density, ncomponents=1);
184  Oxygen->AddElement(elO, 2);
185
186  density = 1.2928*mg/cm3 ;       // STP
187  G4Material* Air = new G4Material(name="Air"  , density, ncomponents=3);
188  Air->AddMaterial( Nitrogen, fractionmass = 0.7557 ) ;
189  Air->AddMaterial( Oxygen,   fractionmass = 0.2315 ) ;
190  Air->AddMaterial( Argon,    fractionmass = 0.0128 ) ;
191
192  /* **************
193
194  // 93% Kr + 7% CH4, STP
195  density = 3.491*mg/cm3 ;     
196  G4Material* Kr7CH4 = new G4Material(name="Kr7CH4"  , density,
197                                      ncomponents=2);
198  Kr7CH4->AddMaterial( Kr,       fractionmass = 0.986 ) ;
199  Kr7CH4->AddMaterial( metane,   fractionmass = 0.014 ) ;
200
201  G4double TRT_Xe_density = 5.485*mg/cm3;
202  G4Material* TRT_Xe = new G4Material(name="TRT_Xe", TRT_Xe_density, nel=1,
203                                      kStateGas,293.15*kelvin,1.*atmosphere);
204  TRT_Xe->AddElement(elXe,1);
205
206  G4double TRT_CO2_density = 1.842*mg/cm3;
207  G4Material* TRT_CO2 = new G4Material(name="TRT_CO2", TRT_CO2_density, nel=2,
208                                       kStateGas,293.15*kelvin,1.*atmosphere);
209  TRT_CO2->AddElement(elC,1);
210  TRT_CO2->AddElement(elO,2);
211
212  G4double TRT_CF4_density = 3.9*mg/cm3;
213  G4Material* TRT_CF4 = new G4Material(name="TRT_CF4", TRT_CF4_density, nel=2,
214                                       kStateGas,293.15*kelvin,1.*atmosphere);
215  TRT_CF4->AddElement(elC,1);
216  TRT_CF4->AddElement(elF,4);
217
218  // ATLAS TRT straw tube gas mixture (20 C, 1 atm)
219
220  G4double XeCO2CF4_density = 4.76*mg/cm3;
221  G4Material* XeCO2CF4 = new G4Material(name="XeCO2CF4", XeCO2CF4_density,
222                                        ncomponents=3,
223                                        kStateGas,293.15*kelvin,1.*atmosphere);
224  XeCO2CF4->AddMaterial(TRT_Xe,0.807);
225  XeCO2CF4->AddMaterial(TRT_CO2,0.039);
226  XeCO2CF4->AddMaterial(TRT_CF4,0.154);
227
228  *********** */
229
230  // Xenon as detector gas, STP
231  density = 5.858*mg/cm3 ;
232  a = 131.29*g/mole ;
233  G4Material* Xe  = new G4Material(name="Xenon",z=54., a, density );
234
235  // Metane, STP
236  density = 0.7174*mg/cm3 ;
237  G4Material* metane = new G4Material(name="CH4",density,nel=2) ;
238  metane->AddElement(elC,1) ;
239  metane->AddElement(elH,4) ;
240
241  // C3H8,20 C, 2 atm
242  density = 3.758*mg/cm3 ;
243  G4Material* C3H8 = new G4Material(name="C3H8",density,nel=2) ;
244  C3H8->AddElement(elC,3) ;
245  C3H8->AddElement(elH,8) ;
246
247  // Propane, STP
248  density = 2.005*mg/cm3 ;
249  G4Material* propane = new G4Material(name="propane",density,nel=2) ;
250  propane->AddElement(elC,3) ;
251  propane->AddElement(elH,8) ;
252
253  // 87.5% Xe + 7.5% CH4 + 5% C3H8, 20 C, 1 atm
254  density = 4.9196*mg/cm3 ;
255  G4Material* XeCH4C3H8 = new G4Material(name="XeCH4C3H8"  , 
256                                  density,  ncomponents=3);
257  XeCH4C3H8->AddMaterial( Xe,       fractionmass = 0.971 ) ;
258  XeCH4C3H8->AddMaterial( metane,   fractionmass = 0.010 ) ;
259  XeCH4C3H8->AddMaterial( propane,  fractionmass = 0.019 ) ;
260
261  // 93% Ar + 7% CH4, STP
262  density = 1.709*mg/cm3 ;     
263  G4Material* Ar7CH4 = new G4Material(name="Ar7CH4", density, ncomponents=2);
264  Ar7CH4->AddMaterial( Argon,    fractionmass = 0.971 ) ;
265  Ar7CH4->AddMaterial( metane,   fractionmass = 0.029 ) ;
266
267  // Carbon dioxide, STP
268  density = 1.977*mg/cm3;
269  G4Material* CarbonDioxide = new G4Material(name="CO2", density, nel=2);
270  CarbonDioxide->AddElement(elC,1);
271  CarbonDioxide->AddElement(elO,2);
272
273  // 80% Ar + 20% CO2, STP
274  density = 1.8223*mg/cm3 ;     
275  G4Material* Ar_80CO2_20 = new G4Material(name="ArCO2"  , density, 
276                                           ncomponents=2);
277  Ar_80CO2_20->AddMaterial( Argon,           fractionmass = 0.783 ) ;
278  Ar_80CO2_20->AddMaterial( CarbonDioxide,   fractionmass = 0.217 ) ;
279
280  // 80% Xe + 20% CO2, STP
281  density = 5.0818*mg/cm3 ;     
282  G4Material* Xe20CO2 = new G4Material(name="Xe20CO2", density, 
283                                       ncomponents=2);
284  Xe20CO2->AddMaterial( Xe,              fractionmass = 0.922 ) ;
285  Xe20CO2->AddMaterial( CarbonDioxide,   fractionmass = 0.078 ) ;
286
287  // 80% Kr + 20% CO2, STP
288  density = 3.601*mg/cm3 ;     
289  G4Material* Kr20CO2 = new G4Material(name="Kr20CO2"  , density, 
290                                        ncomponents=2);
291  Kr20CO2->AddMaterial( Kr,              fractionmass = 0.89 ) ;
292  Kr20CO2->AddMaterial( CarbonDioxide,   fractionmass = 0.11 ) ;
293
294  // G4cout << *(G4Material::GetMaterialTable()) << G4endl;
295 
296  // fWindowMat = Mylar ;
297 
298  fAbsorberMaterial = XeCH4C3H8;
299  // Al; // Si; // Xe; // Ar7CH4; // C3H8; // XeCH4C3H8;
300
301  fWorldMaterial    = Mylar; // Air ;
302}
303
304/////////////////////////////////////////////////////////////////////////
305//
306//
307 
308G4VPhysicalVolume* Em8DetectorConstruction::ConstructCalorimeter()
309{
310  // Cleanup old geometry
311
312  G4GeometryManager::GetInstance()->OpenGeometry();
313  G4PhysicalVolumeStore::GetInstance()->Clean();
314  G4LogicalVolumeStore::GetInstance()->Clean();
315  G4SolidStore::GetInstance()->Clean();
316
317  //  G4RegionStore::GetInstance()->Clean();
318
319  // complete the Calor parameters definition and print
320
321  ComputeCalorParameters();
322  PrintCalorParameters();
323     
324  // World
325  fSolidWorld = new G4Tubs("World",                             //its name
326                   0.,fWorldSizeR,fWorldSizeZ/2.,0.,twopi)       ;//its size
327                         
328  fLogicWorld = new G4LogicalVolume(fSolidWorld,                //its solid
329                                   fWorldMaterial,      //its material
330                                   "World");            //its name
331                                   
332  fPhysicsWorld = new G4PVPlacement(0,                  //no rotation
333                                 G4ThreeVector(),       //at (0,0,0)
334                                 "World",               //its name
335                                 fLogicWorld,           //its logical volume
336                                 NULL,                  //its mother  volume
337                                 false,                 //no boolean operation
338                                 0);                    //copy number
339
340  // Absorber
341
342  if (fAbsorberThickness > 0.) 
343  { 
344
345    fSolidAbsorber = new G4Tubs("Absorber",             
346                          0.,fAbsorberRadius,fAbsorberThickness/2.,0.,twopi); 
347                         
348    fLogicAbsorber = new G4LogicalVolume(fSolidAbsorber,   
349                                          fAbsorberMaterial, 
350                                          "Absorber");     
351                                         
352    fPhysicsAbsorber = new G4PVPlacement(0,               
353                    G4ThreeVector(0.,0.,fAbsorberZ),       
354                                        "Absorber",       
355                                        fLogicAbsorber,     
356                                        fPhysicsWorld,       
357                                        false,             
358                                        0);               
359                                       
360  }
361  if( fRegGasDet != 0 )  // remove obsolete root logical volume
362  {
363    fRegGasDet->RemoveRootLogicalVolume(fLogicAbsorber);
364  }
365  G4ProductionCuts* cuts = 0;
366
367  if( fRegGasDet == 0 ) // First time - instantiate a region and a cut objects
368  {   
369    fRegGasDet = new G4Region("VertexDetector");
370    cuts = new G4ProductionCuts();
371    fRegGasDet->SetProductionCuts(cuts);
372  }
373  else  // Second time - get a cut object from region
374  {   
375    cuts = fRegGasDet->GetProductionCuts();
376  }
377  fRegGasDet->AddRootLogicalVolume(fLogicAbsorber);                               
378
379  cuts->SetProductionCut(fGammaCut,"gamma");
380  cuts->SetProductionCut(fElectronCut,"e-");
381  cuts->SetProductionCut(fPositronCut,"e+");
382
383  // Sensitive Detectors: Absorber
384 
385  G4SDManager* SDman = G4SDManager::GetSDMpointer();
386
387  if(!fCalorimeterSD)
388  {
389    fCalorimeterSD = new Em8CalorimeterSD("CalorSD",this);
390    SDman->AddNewDetector( fCalorimeterSD );
391  }
392  if (fLogicAbsorber)  fLogicAbsorber->SetSensitiveDetector(fCalorimeterSD);
393
394  // Parameterisation
395
396  //   G4VXrayTRmodel* pTRModel = new G4IrregularXrayTRmodel(logicRadiator,
397  //        fRadThickness,fGasGap);
398
399  //  G4VXrayTRmodel* pTRModel = new G4FoamXrayTRmodel(logicRadiator,
400  //                                    fRadThickness,fGasGap);
401
402  //   G4VXrayTRmodel* pTRModel = new G4RegularXrayTRmodel(logicRadiator,
403  //                                                       fRadThickness,fGasGap);
404
405  //  G4double alphaPlate = 160.0 ;
406  //  G4double alphaGas   = 160.0 ;
407
408  //  G4VXrayTRmodel* pTRModel = new G4GamDistrXrayTRmodel(logicRadiator,
409  //                                           fRadThickness,alphaPlate,
410  //                                              fGasGap,alphaGas);
411
412  //  G4VXrayTRmodel* pTRModel = new G4PlateIrrGasXrayTRmodel(logicRadiator,
413  //       fRadThickness,fGasGap);
414
415  //  pTRModel->GetPlateZmuProduct() ;
416  //  pTRModel->GetGasZmuProduct() ;
417
418  //  pTRModel->GetNumberOfPhotons() ; 
419 
420  // always return physics world
421
422  return fPhysicsWorld;
423}
424
425////////////////////////////////////////////////////////////////////////////
426//
427//
428
429void Em8DetectorConstruction::PrintCalorParameters()
430{
431  G4cout << "\n The  WORLD   is made of " 
432         << fWorldSizeZ/mm << "mm of " << fWorldMaterial->GetName() ;
433  G4cout << ", the transverse size (R) of the world is " << fWorldSizeR/mm
434         << " mm. " << G4endl;
435  G4cout << " The ABSORBER is made of " 
436         << fAbsorberThickness/mm << "mm of " << fAbsorberMaterial->GetName() ;
437  G4cout << ", the transverse size (R) is " << fAbsorberRadius/mm << " mm. " << G4endl;
438  G4cout << " Z position of the (middle of the) absorber " 
439         << fAbsorberZ/mm << "  mm." << G4endl;
440  G4cout << G4endl;
441}
442
443///////////////////////////////////////////////////////////////////////////
444//
445//
446
447void Em8DetectorConstruction::SetAbsorberMaterial(G4String materialChoice)
448{
449  // get the pointer to the material table
450  const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
451
452  // search the material by its name   
453  G4Material* pttoMaterial;
454
455  for (size_t J = 0 ; J < theMaterialTable->size() ; J++)
456  { 
457    pttoMaterial = (*theMaterialTable)[J];
458     
459    if(pttoMaterial->GetName() == materialChoice)
460    {
461      fAbsorberMaterial = pttoMaterial;
462      fLogicAbsorber->SetMaterial(pttoMaterial); 
463
464        // PrintCalorParameters();
465    }             
466  }
467}
468
469////////////////////////////////////////////////////////////////////////////
470//
471//
472
473void Em8DetectorConstruction::SetWorldMaterial(G4String materialChoice)
474{
475  // get the pointer to the material table
476  const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
477
478  // search the material by its name   
479  G4Material* pttoMaterial;
480  for (size_t J=0 ; J<theMaterialTable->size() ; J++)
481  { 
482    pttoMaterial = (*theMaterialTable)[J]; 
483   
484    if(pttoMaterial->GetName() == materialChoice)
485    {
486      fWorldMaterial = pttoMaterial;
487      fLogicWorld->SetMaterial(pttoMaterial);
488 
489       //  PrintCalorParameters();
490    }             
491  }
492}
493
494///////////////////////////////////////////////////////////////////////////
495//
496//
497
498void Em8DetectorConstruction::SetAbsorberThickness(G4double val)
499{
500  // change Absorber thickness and recompute the calorimeter parameters
501  fAbsorberThickness = val;
502  ComputeCalorParameters();
503} 
504
505/////////////////////////////////////////////////////////////////////////////
506//
507//
508
509void Em8DetectorConstruction::SetAbsorberRadius(G4double val)
510{
511  // change the transverse size and recompute the calorimeter parameters
512  fAbsorberRadius = val;
513  ComputeCalorParameters();
514} 
515
516////////////////////////////////////////////////////////////////////////////
517//
518//
519
520void Em8DetectorConstruction::SetWorldSizeZ(G4double val)
521{
522  fWorldChanged=true;
523  fWorldSizeZ = val;
524  ComputeCalorParameters();
525} 
526
527///////////////////////////////////////////////////////////////////////////
528//
529//
530
531void Em8DetectorConstruction::SetWorldSizeR(G4double val)
532{
533  fWorldChanged=true;
534  fWorldSizeR = val;
535  ComputeCalorParameters();
536} 
537
538//////////////////////////////////////////////////////////////////////////////
539//
540//
541
542void Em8DetectorConstruction::SetAbsorberZpos(G4double val)
543{
544  fAbsorberZ  = val;
545  ComputeCalorParameters();
546} 
547
548
549///////////////////////////////////////////////////////////////////////////////
550//
551//
552 
553void Em8DetectorConstruction::UpdateGeometry()
554{
555  G4RunManager::GetRunManager()->DefineWorldVolume(ConstructCalorimeter());
556}
557
558//
559//
560////////////////////////////////////////////////////////////////////////////
561
Note: See TracBrowser for help on using the repository browser.