source: JEM-EUSO/esaf_lal/branches/camille/packages/simulation/detector/G4Detector/optics/src/G4TusGeometry.cc @ 184

Last change on this file since 184 was 184, checked in by moretto, 11 years ago

changes from biktem r:3032
File size: 9.8 KB
Line 
1#include "G4TusGeometry.hh"
2#include <G4Sphere.hh>
3#include <G4Material.hh>
4#include <G4MaterialTable.hh>
5#include <G4Element.hh>
6#include <G4ElementTable.hh>
7#include <G4LogicalBorderSurface.hh>
8#include <G4LogicalSkinSurface.hh>
9#include <G4OpticalSurface.hh>
10#include <G4Box.hh>
11#include <G4LogicalVolume.hh>
12#include <G4RotationMatrix.hh>
13#include <G4ThreeVector.hh>
14#include <G4Transform3D.hh>
15#include <G4PVPlacement.hh>
16#include <G4OpBoundaryProcess.hh>
17#include <G4VisAttributes.hh>
18#include <G4Cons.hh>
19#include <G4SubtractionSolid.hh>
20#include <G4IntersectionSolid.hh>
21#include <G4UnionSolid.hh>
22#include <G4Paraboloid.hh>
23#include <G4Tubs.hh>
24#include <G4Polyhedra.hh>
25
26//#include <Geant4GM/volumes/Factory.h>
27//#include <RootGM/volumes/Factory.h>
28
29
30#include <TGeoManager.h>
31#include <TMath.h>
32#include <TEnv.h>
33
34using namespace TMath;
35
36G4TusGeometry::G4TusGeometry(){
37    expHall_z = 5000.*mm;
38    expHall_x = expHall_y = 1600.*mm;
39    cout<<"OPTDETECTOR TUS"<<endl;
40}
41
42//_______________________________________________________________________________
43G4TusGeometry::~G4TusGeometry(){
44
45}
46
47//_______________________________________________________________________________
48G4VPhysicalVolume* G4TusGeometry::Construct(){
49    G4double  a, z, density;
50    G4int nelements;
51
52    G4Element* Ag = new G4Element("Argentum", "Ag", z=47 , a=107.868*g/mole);
53
54    G4Material* Vacuum =
55    new G4Material("Galactic", z=1., a=1.01*g/mole,density= universe_mean_density, kStateGas,2.73*kelvin, 3.e-18*pascal);
56
57    G4Material* metal= new G4Material("metal", density=10.49*g/cm3, nelements=1);
58    metal->AddElement(Ag, 1);
59
60    //metal
61    const G4int nEntries = 2;
62    G4double PhotonEnergy[nEntries] =
63                {0.001*eV,50.*eV};/*0.2034*eV, 0.2068*eV, 0.2103*eV, 0.2139*eV,
64                0.2177*eV, 0.2216*eV, 0.2256*eV, 0.2298*eV,
65                0.2341*eV, 0.2386*eV, 0.2433*eV, 0.2481*eV,
66                0.2532*eV, 0.2585*eV, 0.2640*eV, 0.2697*eV,
67                0.2757*eV, 0.2820*eV, 0.2885*eV, 0.2954*eV,
68                0.3026*eV, 0.3102*eV, 0.3181*eV, 0.3265*eV,
69                0.3353*eV, 0.3446*eV, 0.3545*eV, 0.3649*eV,
70                0.3760*eV, 0.3877*eV, 0.4002*eV, 0.9136*eV };*/
71
72    G4double RefractiveIndex3[nEntries] =
73            { 0., 0.};//, 0., 0.,0., 0., 0.,
74            /*0., 0., 0., 0.,0., 0., 0.,
75            0., 0., 0., 0.,0., 0., 0.,
76            0., 0., 0., 0.,0., 0., 0.,
77            0., 0., 0., 0. };*/
78
79    G4double  Reflectivity3[nEntries] =
80            { 0.99, 0.99};//, 0.99, 0.99,0.99, 0.99, 0.99,
81            /*0.99, 0.99, 0.99, 0.99,0.99, 0.99, 0.99,
82            0.99, 0.99, 0.99, 0.99,0.99, 0.99, 0.99,
83            0.99, 0.99, 0.99, 0.99,0.99, 0.99, 0.99,
84            0.99, 0.99, 0.99, 0.99 };*/
85
86    G4MaterialPropertiesTable* myMPT3 = new G4MaterialPropertiesTable();
87    myMPT3->AddProperty("RINDEX", PhotonEnergy, RefractiveIndex3, nEntries);
88    myMPT3->AddProperty("REFLECTIVITY",PhotonEnergy,Reflectivity3,nEntries);
89    metal->SetMaterialPropertiesTable(myMPT3);
90
91    //Vacuum
92    G4double RefractiveIndex1[nEntries] =
93            { 1.00, 1.00 };/*, 1.00, 1.00, 1.00, 1.00, 1.00,
94            1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00,
95            1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00,
96            1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00,
97            1.00, 1.00, 1.00, 1.00 };*/
98
99    G4double  Reflectivity1[nEntries] =
100            { 0., 0.};/*, 0., 0.,0., 0., 0.,
101            0., 0., 0., 0.,0., 0., 0.,
102            0., 0., 0., 0.,0., 0., 0.,
103            0., 0., 0., 0.,0., 0., 0.,
104            0., 0., 0., 0. };*/
105
106    G4MaterialPropertiesTable* myMPT1 = new G4MaterialPropertiesTable();
107    myMPT1->AddProperty("RINDEX", PhotonEnergy, RefractiveIndex1, nEntries);
108    myMPT1->AddProperty("REFLECTIVITY", PhotonEnergy,Reflectivity1 , nEntries);
109    Vacuum->SetMaterialPropertiesTable(myMPT1);
110
111    //world
112    G4Box* expHall = new G4Box("World",expHall_x,expHall_y,expHall_z);
113    G4LogicalVolume* expHall_log = new G4LogicalVolume(expHall,Vacuum,"World",0,0,0);
114    G4VPhysicalVolume* expHall_phys = new G4PVPlacement(0,G4ThreeVector(0,0,0),expHall_log,"World",0,false,0);
115    expHall_log->SetVisAttributes (G4VisAttributes::Invisible);
116
117    // tus mirror
118    const int ii = 11;
119    double PRmin[ii] = {245.764,348.712,428.486,496.387,556.776,611.882,663.023,711.056,756.571,800.000,841.665};
120    double PRmax[ii] = {347.563,427.083,494.773,554.977,609.918,660.908,708.802,754.188,797.496,839.047,879.090};
121
122
123    //1-10 rings
124    double dz1;
125    double dz2;
126    double dz;
127    G4UnionSolid* segment1[ii-1];
128    G4UnionSolid* segment2[ii];
129
130    G4SubtractionSolid* segment11;
131    for(int i=0; i<ii-1; i++){
132        dz1 = ((PRmax[i]*PRmax[i])/(4.*(1500.+(10.*(i+1.)))))-(10.*(i+1.));
133        dz2 = ((PRmin[i]*PRmin[i])/(4.*(1500.+(10.*(i+1.)))))-(10.*(i+1.));
134        dz=(dz1-dz2);
135        segment1[i] = OneSegment(PRmin[i],PRmax[i],dz*0.5,PRmin[i+1]);
136    }
137
138    //11 ring
139    dz1 = ((PRmax[10]*PRmax[10])/(4.*(1500.+(10.*(10.+1.)))))-(10.*(10.+1.));
140    dz2 = ((PRmin[10]*PRmin[10])/(4.*(1500.+(10.*(10.+1.)))))-(10.*(10.+1.));
141
142    dz=(dz1-dz2);
143    G4Paraboloid* par1 = new G4Paraboloid("par1",dz*0.5*mm,PRmin[10]*mm,PRmax[10]*mm);
144    G4Tubs* tubs = new G4Tubs("tubs1",PRmin[10]*mm,PRmax[10]*mm,dz*0.5*mm,0.*deg,360.*deg);
145    segment11 = new G4SubtractionSolid("kink",tubs,par1);
146
147
148    segment2[0] = new G4UnionSolid("ring", segment1[0],segment1[1]);
149    for(int j=1;j<9;j++)segment2[j] = new G4UnionSolid("ring", segment2[j-1],segment1[j+1]);
150    segment2[9] = new G4UnionSolid("ring", segment2[8],segment11);
151
152
153    dz=10.*mm;
154    G4RotationMatrix*  pRot= new G4RotationMatrix();
155    pRot->rotateX(180.*DegToRad());
156
157    G4Paraboloid* seg0 = new G4Paraboloid("seg0",dz*0.5*mm,0.*mm,PRmin[0]*mm);
158    segment2[10] = new G4UnionSolid("ring", segment2[9],seg0,pRot,G4ThreeVector(0,0,0));
159
160    double p1[2]={-500.*mm,500.*mm};
161    double p2[2]={0.*mm,0.*mm};
162    double iside=sqrt(330.*330.-165.*165.)*mm;
163    double oside=330.*mm;
164    double p3[2]={iside,iside};
165
166    G4Polyhedra* p = new G4Polyhedra("polly",0.,360.*deg,6,2,p1,p2,p3);
167
168//     EInside ins1 = p->Inside(G4ThreeVector(290.*mm,0.,0.));
169//     EInside ins2 = p->Inside(G4ThreeVector(0.*mm,290.*mm,0.));
170
171    double sch = 0.1*sqrt(2.)*mm;
172    G4IntersectionSolid* s[7];
173     s[0] = new G4IntersectionSolid("s_0",segment2[10], p);
174     s[1] = new G4IntersectionSolid("s_1",segment2[10], p,0,G4ThreeVector(0.,2.*iside+0.1,0.));
175     s[2] = new G4IntersectionSolid("s_2",segment2[10], p,0,G4ThreeVector(0.,-2.*iside-0.1,0.));
176     s[3] = new G4IntersectionSolid("s_3",segment2[10], p,0,G4ThreeVector(1.5*oside+sch,iside+sch,0.));
177     s[4] = new G4IntersectionSolid("s_4",segment2[10], p,0,G4ThreeVector(1.5*oside+sch,-iside-sch,0.));
178     s[5] = new G4IntersectionSolid("s_5",segment2[10], p,0,G4ThreeVector(-1.5*oside-sch,iside+sch,0.));
179     s[6] = new G4IntersectionSolid("s_6",segment2[10], p,0,G4ThreeVector(-1.5*oside-sch,-iside-sch,0.));
180
181    G4LogicalVolume* s_log[7];
182    for(int i(0);i<7;i++)
183    s_log[i] = new G4LogicalVolume(s[i],metal,Form("s%d_log",i),0,0,0);
184
185    G4VPhysicalVolume* s_phys[7];
186
187    for(int i(0);i<7;i++)
188    s_phys[i] =new G4PVPlacement(0,G4ThreeVector(0.,0.,-1510.*mm),s_log[i],Form("mirror_phys%d",i),expHall_log,false,0);
189
190    const G4int n=2;
191    G4double PhotonEnergyS[n]    = { 0.01034*eV, 5.9136*eV};
192    G4double RefractiveIndexS[n] = { 0.0, 0.0};
193    G4double  ReflectivityS[n]   = { 0.99, 0.99};
194    G4double efficiency[n]       = {0.99, 0.99};
195
196    G4MaterialPropertiesTable* myST2 = new G4MaterialPropertiesTable();
197    myST2->AddProperty("RINDEX",  PhotonEnergyS, RefractiveIndexS,n);
198    myST2->AddProperty("REFLECTIVITY",PhotonEnergyS,ReflectivityS,n);
199    myST2 -> AddProperty("EFFICIENCY",PhotonEnergyS,efficiency,n);
200
201
202    //Optical Surface
203    G4OpticalSurface* OpSurface[7];
204    G4LogicalSkinSurface* Surface[7];
205    for(int i(0);i<7;i++){
206        OpSurface[i]= new G4OpticalSurface(Form("GlassSurface%d",i));
207        Surface[i] = new G4LogicalSkinSurface(Form("GlassSurface%d",i),s_log[i],OpSurface[i]);
208    //
209        OpSurface[i]->SetType(dielectric_metal);
210        OpSurface[i]->SetFinish(polished);
211        OpSurface[i]->SetModel(glisur);
212
213        OpSurface[i]->SetMaterialPropertiesTable(myST2);
214    }
215
216
217    G4LogicalVolume* absorber_log ;
218    G4VPhysicalVolume* absorber_phys;
219    double x= 120.*mm;
220    G4Box* absorber = new G4Box("absorber",x,x,1.*mm);
221    absorber_log = new G4LogicalVolume(absorber,Vacuum,"absorber",0,0,0);
222    absorber_phys = new G4PVPlacement(0,G4ThreeVector(0.,0.,-14.*mm),absorber_log,"absorber",expHall_log,false,0);
223
224
225//
226//     G4VisAttributes* visatr= new G4VisAttributes(G4Colour(1.0,1.0,1.0));
227//     visatr->SetForceAuxEdgeVisible (1);
228//     visatr->SetVisibility(true);
229//     visatr->SetColour(0.99,0.8,0.8,.8);
230//     FMirror_log->SetVisAttributes(visatr);
231//     absorber_log->SetVisAttributes(visatr);
232
233//    // Export VGM geometry to Root
234//    //
235//    Geant4GM::Factory g4Factory;
236//    g4Factory.SetDebug(1);
237//    g4Factory.Import(expHall_phys);
238//
239//    RootGM::Factory rtFactory;
240//    rtFactory.SetDebug(1);
241//    g4Factory.Export(&rtFactory);
242//    gGeoManager->CloseGeometry();
243//    gGeoManager->Export("SphericalLens.root");
244    return expHall_phys;
245}
246
247//_______________________________________________________________________________
248G4UnionSolid* G4TusGeometry::OneSegment(double  PRmin,double PRmax,double  PDz,double CRmax){
249    G4Paraboloid* par1 = new G4Paraboloid("par1",PDz*mm,PRmin*mm,PRmax*mm);
250    G4Tubs* tub1 = new G4Tubs("tubs1",PRmin*mm,PRmax*mm,PDz*mm,0.*deg,360.*deg);
251    G4SubtractionSolid* sub_zub = new G4SubtractionSolid("kink",tub1,par1);
252
253    G4Cons* cons = new G4Cons("cons", (PRmax-0.5)*mm,CRmax*mm,(PRmax-0.5)*mm,PRmax*mm, PDz*mm,0.*deg,360.*deg);
254    G4Tubs* tub2 = new G4Tubs("tubs2",(PRmax-0.5)*mm,PRmax*mm,PDz*mm,0.*deg,360.*deg);
255    G4SubtractionSolid* conus = new G4SubtractionSolid("conus",cons,tub2);
256    G4UnionSolid* segment =  new G4UnionSolid("ring", sub_zub, conus);
257    return segment;
258}
Note: See TracBrowser for help on using the repository browser.