source: JEM-EUSO/esaf_lal/branches/camille/packages/simulation/detector/G4Detector/optics/src/G4EusoGeometryShape.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: 17.5 KB
Line 
1#include "G4EusoGeometryShape.hh"
2
3#include "FresnelSurface.h"
4#include "FresnelSurfaceHelper.h"
5#include "ParabolicSurface.h"
6#include "SphericSurface.h"
7#include "FlatSurface.h"
8#include "ShapeOfLens.h"
9#include "MaterialParameters.hh"
10
11#include <G4AssemblyVolume.hh>
12#include <G4Box.hh>
13#include <G4IntersectionSolid.hh>
14#include <G4Cons.hh>
15#include <G4DisplacedSolid.hh>
16#include <G4Element.hh>
17#include <G4ElementTable.hh>
18#include <G4GeometryTolerance.hh>
19#include <G4IntersectionSolid.hh>
20#include <G4LogicalBorderSurface.hh>
21#include <G4LogicalSkinSurface.hh>
22#include <G4LogicalVolume.hh>
23#include <G4Material.hh>
24#include <G4MaterialTable.hh>
25#include <G4OpBoundaryProcess.hh>
26#include <G4PVPlacement.hh>
27#include <G4Paraboloid.hh>
28#include <G4Sphere.hh>
29#include <G4SubtractionSolid.hh>
30#include <G4Torus.hh>
31#include <G4Transform3D.hh>
32#include <G4Tubs.hh>
33
34#include <TEnv.h>
35#include <TGraph.h>
36#include <TGeoManager.h>
37#include <TSystem.h>
38#include <TF1.h>
39#include <TFile.h>
40#include <cmath>
41#include <iostream>
42using namespace std;
43using namespace TMath;
44using namespace G4FresnelLens;
45
46#include "Ntrace_optics.hh"
47using namespace NTraceLens;
48//using std::sqrt;
49//using std::cout;
50//using std::flush;
51
52// #include <Geant4GM/volumes/Factory.h>
53// #include <RootGM/volumes/Factory.h>
54
55#include <stdlib.h>
56
57// #define DISABLE_1_LENSphys
58// #define DISABLE_2_LENSphys
59// #define DISABLE_3_LENSphys
60
61G4EusoGeometryShape::G4EusoGeometryShape()
62{
63    fGeomTolerance=G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
64    expHall_x = expHall_y = 2000.*mm;
65    expHall_z = 4000.*mm;
66    fit = Conf()->GetBool("G4EusoGeometryShape.UseFitTransmittivity");
67    CYTOP=0;
68    PMMA=0;
69    fOpticalSurfaceC=0;
70    fOpticalSurfaceP=0;
71
72
73}
74
75G4EusoGeometryShape::~G4EusoGeometryShape(){
76    if(CYTOP)delete CYTOP;
77    if(PMMA)delete PMMA;
78    if(fOpticalSurfaceC)delete fOpticalSurfaceC;
79    if(fOpticalSurfaceP)delete fOpticalSurfaceP;
80
81}
82
83
84G4VPhysicalVolume* G4EusoGeometryShape::Construct(){
85    lens_dir =fOpticalSystemParam->lens_dir;
86#ifdef DEBUG
87    printf("%p %p\n",fOpticalSystemParam, fOpticalSystemParam->lens_dir);
88    printf("%s\n",lens_dir);
89    printf("%s \n", fOpticalSystemParam->name);
90#endif
91    G4double  a, z, density;
92    G4int nelements;
93    /// elements
94    G4Element* H = new G4Element("Hydrogen", "H", z=1 , a=1.00794*g/mole);
95    G4Element* O = new G4Element("Oxygen"  , "O", z=8 , a=16.00*g/mole);
96    G4Element* C = new G4Element("Carbon", "C", z=6 , a=12.0107*g/mole);
97    G4Element* P = new G4Element("Phosphorus", "P", z=15 , a=30.973761*g/mole);
98
99    G4Material* Tributylphosphine = new G4Material("Tributylphosphine", density=0.81*g/cm3, nelements=3);
100    Tributylphosphine->AddElement(P, 1);
101    Tributylphosphine->AddElement(C, 12);
102    Tributylphosphine->AddElement(H, 27);
103
104    G4Material* Dibutylphosphine = new G4Material("Dibutylphosphine", density=0.81*g/cm3, nelements=3);
105    Dibutylphosphine->AddElement(P, 1);
106    Dibutylphosphine->AddElement(C, 8);
107    Dibutylphosphine->AddElement(H, 19);
108
109    G4Material* TributylphosphineOxide = new G4Material("TributylphosphineOxide", density=0.81*g/cm3, nelements=4);
110    TributylphosphineOxide->AddElement(P, 1);
111    TributylphosphineOxide->AddElement(C, 12);
112    TributylphosphineOxide->AddElement(H, 27);
113    TributylphosphineOxide->AddElement(O, 1);
114    ///material
115    CYTOP = new G4Material("CYTOP", density=2.03*g/cm3, nelements=3);
116    CYTOP->AddMaterial(Tributylphosphine,0.98);
117    CYTOP->AddMaterial(TributylphosphineOxide,0.005);
118    CYTOP->AddMaterial(Dibutylphosphine,0.015);
119
120    PMMA = new G4Material("PMMA", density=1.19*g/cm3, nelements=3);
121    PMMA->AddElement(O,2);
122    PMMA->AddElement(H,8);
123    PMMA->AddElement(C,5);
124    // Vacuum-world; Vacuum1-FS; Vacuum2-diffraction
125    G4Material* Vacuum =
126    new G4Material("Galactic", z=1., a=1.01*g/mole,density= universe_mean_density,kStateGas, 2.73*kelvin, 3.e-18*pascal);
127    G4Material* Vacuum1 =
128    new G4Material("Galactic1", z=1., a=1.01*g/mole,density= 10.*g/cm3);
129    G4Material* Vacuum2 =
130    new G4Material("Galactic", z=1., a=1.01*g/mole,density= universe_mean_density,kStateGas, 2.73*kelvin, 3.e-18*pascal);
131
132    //FIXME
133    //Iris material, now IRIS is just black body, cover also is a black body
134    G4Material* fe_mat= new G4Material("fe_mat", z=26. , a=55.845*g/mole, density=7.87*g/cm3);
135
136    ///Material Properties Table
137
138    const G4int nEntries = 18;
139    const double lambdas1[nEntries]={245.,275.,305., 310.,320., 330., 335., 340.,
140                                    350., 360., 370., 380., 390., 400., 410., 420., 435., 546.} ;
141
142    double cnst=hbarc*twopi/eV;
143    G4double PhotonEnergy1[nEntries];
144    for (int i(0); i<nEntries; i++){
145        PhotonEnergy1[i]=cnst/lambdas1[nEntries-i-1];
146    }
147
148    G4double RefractiveIndex1[nEntries] =
149            { 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00,
150              1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00,
151              1.00, 1.00, 1.00, 1.00};
152
153     G4double  EUSO_DINDEX[nEntries] =
154            { 0., 0., 0., 0.,0., 0., 0.,
155              0., 0., 0., 0.,0., 0., 0.,
156              0., 0., 0., 0.};
157
158    G4MaterialPropertiesTable* myMPT1 = new G4MaterialPropertiesTable();
159    myMPT1->AddProperty("RINDEX", PhotonEnergy1, RefractiveIndex1, nEntries);
160    Vacuum->SetMaterialPropertiesTable(myMPT1);
161    Vacuum1->SetMaterialPropertiesTable(myMPT1);
162
163    fOpticalSurfaceC = new G4OpticalSurface("vacuum_cytop_os");
164    fOpticalSurfaceC->SetType(dielectric_dielectric);
165    fOpticalSurfaceC->SetModel(unified);
166    fOpticalSurfaceC->SetFinish(polished);
167
168    fOpticalSurfaceP = new G4OpticalSurface("vacuum_pmma_os");
169    fOpticalSurfaceP->SetType(dielectric_dielectric);
170    fOpticalSurfaceP->SetModel(unified);
171    //fOpticalSurfaceP->SetFinish(polishedbackpainted);
172    fOpticalSurfaceP->SetFinish(polished);
173
174
175    MaterialParameters* mp1 = new MaterialParameters();
176    MaterialParameters* mp2 = new MaterialParameters();
177    char m1mat[50];
178    char m2mat[50];
179
180    SetMaterial(2,mp1,m1mat);  //2 and 3 number of materials, 1&2-air,vacuum. 2&3-materials of lenses
181    SetMaterial(3,mp2,m2mat);
182
183
184    /// material properties to use with OpESAFDiffraction
185    /// there will be thin vacuum layer after 2-d lens with normal vacuum properties
186    /// but special properties of PMMA
187    /// so G4OpBoundaryProcess will see vacuum and do nothing
188    /// and OpESAFDiffraction will see PMMA and do special diffraction
189
190    G4MaterialPropertiesTable* myMPT3a = new G4MaterialPropertiesTable();
191    myMPT3a->AddProperty("RINDEX", PhotonEnergy1, RefractiveIndex1, nEntries);
192    myMPT3a->AddProperty("EUSO_DINDEX",PhotonEnergy1,EUSO_DINDEX,nEntries);
193    //diffractive vacuum
194    //for the EUSO diffraction it has RINDEX properties as second lens
195    char* matdoe = fOpticalSystemParam->material[fOpticalSystemParam->surface[5].matid].name;
196    G4MaterialPropertiesTable* mprop=0;
197
198    if(!strcmp(matdoe,m1mat)){
199        mprop=mp1->GetTable();
200     }
201    else if(!strcmp(matdoe,m2mat)){
202        mprop=mp2->GetTable();
203    }
204    else cerr<<"Unknown doe material"<<endl;
205
206    G4MaterialPropertyVector *mv = mprop->GetProperty ("RINDEX");
207    myMPT3a->AddProperty ("EUSO_RINDEX",mv);
208    Vacuum2->SetMaterialPropertiesTable(myMPT3a);
209
210    //volumes
211    ///world
212    G4Box* expHall = new G4Box("World",expHall_x,expHall_y,expHall_z);
213    G4LogicalVolume* expHall_log = new G4LogicalVolume(expHall,Vacuum,"World",0,0,0);
214    G4VPhysicalVolume* expHall_phys = new G4PVPlacement(0,G4ThreeVector(),expHall_log,"World",0,false,0);
215
216    //intersection volume for lenses
217    double r_wall=fOpticalSystemParam->r_wall;
218    double cut =fOpticalSystemParam->surface[1].r_cut;
219    G4Box* box = new G4Box("box",(r_wall+10.)*mm,cut*mm,1.*m);
220   
221   
222    /// 1-st lens
223
224    double c_rad, c_z1, c_z2, c_dz;
225    char* name=fOpticalSystemParam->name;
226   
227
228    double surf1=fOpticalSystemParam->surface[1].Sz;
229    double surf2=fOpticalSystemParam->surface[2].Sz;
230    double thin= surf2 - surf1;
231    double Rc1=fOpticalSystemParam->surface[1].Rc;
232    double Rc2=fOpticalSystemParam->surface[2].Rc;
233    char* mat = fOpticalSystemParam->material[fOpticalSystemParam->surface[2].matid].name;
234 //
235    //const char* lens_file="output/surfaces.root";
236    FresnelSurfaceHelper* shelper=new FresnelSurfaceHelper();
237    double filter=0.5;
238    if ( !strcmp(name,"CPC_2007") ){
239        filter=DBL_MAX;
240        //printf("Disabling filter for (%s)\n", name);
241    }
242
243//     const char* lens_file="output/surfaces.root";
244#ifndef DISABLE_1_LENS
245    ParametricSurface* uplimit=new SphericSurface(Rc1*mm,0.,r_wall);
246    ParametricSurface* downlimit=new SphericSurface(Rc1*mm,0.,r_wall);
247
248    shelper->SetOptions(0.*mm, -1., -filter);
249    shelper->Read(Form("%s/S1_%s.dat",lens_dir,name));
250    FresnelSurface* downsurf1=shelper->CreateFresnelLens(downlimit, uplimit);
251
252    uplimit=new SphericSurface(Rc2*mm,0.,r_wall);
253    downlimit=new SphericSurface(Rc2*mm,0.,r_wall);
254
255    shelper->SetOptions(thin*mm, -1., filter);
256    shelper->Read(Form("%s/S2_%s.dat",lens_dir,name));
257    FresnelSurface* upsurf1=shelper->CreateFresnelLens(downlimit, uplimit);
258    ShapeOfLens* lens1_solid = new ShapeOfLens("FresnelLens_1", (thin*0.5+0.5)*mm, (r_wall+1.)*mm, upsurf1, downsurf1);
259
260    G4IntersectionSolid* intersec1 = new G4IntersectionSolid("box*FresnelLens_1", box, lens1_solid);
261    G4LogicalVolume* lens1_log = new G4LogicalVolume(intersec1,!strcmp(mat,"CYTOP")?CYTOP:PMMA, "FresnelLens_1",0,0,0);
262   
263    #ifndef DISABLE_1_LENSphys
264    G4VPhysicalVolume* lens1_phys = new G4PVPlacement(0, G4ThreeVector(0.,0.,surf1), lens1_log, "FresnelLens_1", expHall_log,false,0);
265    /*G4LogicalBorderSurface diffr_surf =*/
266    new G4LogicalBorderSurface("diffractor_surf1",lens1_phys,expHall_phys,!strcmp(mat,"CYTOP")?fOpticalSurfaceC:fOpticalSurfaceP);
267    #endif
268#endif
269
270    ///2-d lens
271#ifndef DISABLE_2_LENS
272
273    surf1 = fOpticalSystemParam->surface[4].Sz;
274    surf2 = fOpticalSystemParam->surface[5].Sz;
275    thin = surf2 -surf1;
276    Rc1 = fOpticalSystemParam->surface[4].Rc;
277    Rc2 = fOpticalSystemParam->surface[5].Rc;
278    mat = fOpticalSystemParam->material[fOpticalSystemParam->surface[5].matid].name;
279    cut =fOpticalSystemParam->surface[4].r_cut;
280
281    if(Rc2>1.0e10){
282         uplimit=new FlatSurface();
283         downlimit=new FlatSurface();
284    }
285    else{
286        uplimit=new SphericSurface(Rc1*mm,0.,r_wall);
287        downlimit=new SphericSurface(Rc1*mm,0.,r_wall);
288    }
289    shelper->SetOptions(0., -1., -filter);
290    shelper->Read(Form("%s/S4_%s.dat",lens_dir,name));
291    FresnelSurface* downsurf2=shelper->CreateFresnelLens(downlimit, uplimit);
292    double rsurf5=downsurf2->GetR();
293
294    ParametricSurface* upsurf2;
295
296    if(Rc2>1.0e10){
297         upsurf2 = new FlatSurface((thin)*mm,rsurf5*mm);
298    }
299    else{
300        upsurf2=new SphericSurface(Rc2*mm,(thin+Rc2)*mm,rsurf5*mm);
301    }
302    ShapeOfLens* lens2_solid = new ShapeOfLens("FresnelLens_2", (thin*0.5+0.5)*mm, (rsurf5+1.)*mm, upsurf2, downsurf2);
303    G4IntersectionSolid* intersec2 = new G4IntersectionSolid("box*FresnelLens_2", box, lens2_solid);
304   
305    G4LogicalVolume* lens2_log = new G4LogicalVolume(intersec2,!strcmp(mat,"CYTOP")?CYTOP:PMMA, "FresnelLens_2",0,0,0);
306
307    double diffr_rad=rsurf5; //FIXME: diffractive surface is shorter than surf5
308    double diffr_dz=1.e-8;
309    G4VSolid* euso_diffractor;
310
311    if(Rc2>1.0e10){
312         euso_diffractor=new G4Tubs("euso_diffractor_tubs", 0., diffr_rad, diffr_dz,0.*deg,360.*deg);
313    }
314    else {
315        ParametricSurface* upsurf_diffr = new SphericSurface(Rc2*mm,Rc2*mm+diffr_dz,rsurf5*mm);
316        ParametricSurface* downsurf_diffr = new SphericSurface(Rc2*mm,Rc2*mm-diffr_dz,rsurf5*mm);
317        euso_diffractor=new ShapeOfLens("euso_diffractor_par", diffr_dz, diffr_rad, upsurf_diffr, downsurf_diffr);
318       
319    }
320    G4IntersectionSolid* intersec_diffr = new G4IntersectionSolid("box*euso_diffractor_par", box, euso_diffractor);
321    G4LogicalVolume* euso_diffractor_log = new G4LogicalVolume(intersec_diffr,Vacuum2,"euso_diffractor",0,0,0);
322
323    #ifndef DISABLE_2_LENSphys
324    G4VPhysicalVolume* lens2_phys = new G4PVPlacement(0, G4ThreeVector(0.,0.,surf1), lens2_log, "FresnelLens_2", expHall_log,false,0);
325    //G4VPhysicalVolume* euso_diffractor_phys =
326            new G4PVPlacement(0,G4ThreeVector(0,0,thin+diffr_dz+fGeomTolerance+surf1),
327                                              euso_diffractor_log,"euso_diffractor",expHall_log,false,0);
328
329    /*G4LogicalBorderSurface diffr_surf =*/
330    new G4LogicalBorderSurface("diffractor_surf21",lens2_phys,expHall_phys,!strcmp(mat,"CYTOP")?fOpticalSurfaceC:fOpticalSurfaceP);
331    #endif
332#endif
333
334    /// 3-d Lens
335#ifndef DISABLE_3_LENS
336    surf1 = fOpticalSystemParam->surface[6].Sz;
337    surf2 = fOpticalSystemParam->surface[7].Sz;
338    thin = surf2 -surf1;
339    Rc1 = fOpticalSystemParam->surface[6].Rc;
340    Rc2 = fOpticalSystemParam->surface[7].Rc;
341    mat = fOpticalSystemParam->material[fOpticalSystemParam->surface[7].matid].name;
342    cut =fOpticalSystemParam->surface[6].r_cut;
343
344     uplimit=new SphericSurface(Rc1,0.,r_wall);
345     downlimit=new SphericSurface(Rc1,0.,r_wall);
346     shelper->SetOptions(0., -1., -filter);
347     shelper->Read(Form("%s/S6_%s.dat",lens_dir,name));
348     FresnelSurface* downsurf3=shelper->CreateFresnelLens(downlimit, uplimit);
349
350     uplimit=new SphericSurface(Rc2,0.,r_wall);
351     downlimit=new SphericSurface(Rc2,0.,r_wall);
352     shelper->SetOptions(thin*mm, -1., -filter);
353     shelper->Read(Form("%s/S7_%s.dat",lens_dir,name));
354     FresnelSurface* upsurf3=shelper->CreateFresnelLens(downlimit, uplimit);
355
356     ShapeOfLens* lens3_solid = new ShapeOfLens("FresnelLens_3",(thin*0.5+0.5)*mm, (r_wall+1.)*mm, upsurf3, downsurf3);
357     G4IntersectionSolid* intersec3 = new G4IntersectionSolid("box*FresnelLens_3", box, lens3_solid);
358
359     G4LogicalVolume* lens3_log = new G4LogicalVolume(intersec3,!strcmp(mat,"CYTOP")?CYTOP:PMMA, "FresnelLens_3",0,0,0);
360
361     lens3_solid->GetCylinderLimits(c_rad, c_z1, c_z2);
362     c_dz=(fabs(c_z2)+fabs(c_z1));
363
364
365    #ifndef DISABLE_3_LENSphys
366    G4VPhysicalVolume* lens3_phys = new G4PVPlacement(0, G4ThreeVector(0.,0.,surf1), lens3_log, "FresnelLens_3", expHall_log,false,0);
367    /*G4LogicalBorderSurface diffr_surf =*/
368    new G4LogicalBorderSurface("diffractor_surf3",lens3_phys,expHall_phys,!strcmp(mat,"CYTOP")?fOpticalSurfaceC:fOpticalSurfaceP);
369    #endif
370#endif
371    delete shelper;
372
373    /// fs
374    double fs_dz=10.*mm;
375    double fs_pos= ((fOpticalSystemParam->FS_OFFSET-surf2)*0.5+surf2)*mm;//surf1+c_dz+15.*mm+fs_dz;
376    double fs_rad=fOpticalSystemParam->r_fs*mm;
377
378    G4Tubs* absorber_cyl=new G4Tubs("absorber",0.,fs_rad,fs_dz,0.*deg,360.*deg);
379    G4IntersectionSolid* intersec_fs = new G4IntersectionSolid("box*absorber", box, absorber_cyl);
380    G4LogicalVolume* absorber_cyl_log = new G4LogicalVolume(intersec_fs,Vacuum1,"absorber",0,0,0);
381    //G4VPhysicalVolume* absorber_cyl_phys =
382            new G4PVPlacement(0,G4ThreeVector(0,0,fs_pos),absorber_cyl_log,"absorber",expHall_log,false,0);
383
384
385    ///IRIS
386    bool useiris=Conf()->GetBool("G4EusoGeometryShape.UseIRIS");
387    if ( useiris ){
388        double iris_r_inner=fOpticalSystemParam->surface[3].r_lens;
389        double iris_r_outer=fOpticalSystemParam->r_wall;
390        double iris_dz=fOpticalSystemParam->surface[3].Sz;
391        //double iris_cut=fOpticalSystemParam->surface[3].r_cut;
392        double iris_thickness=1.*mm;
393        G4Tubs* iris_cyl=new G4Tubs("iris",iris_r_inner,iris_r_outer,iris_thickness,0.*deg,360.*deg);
394        G4IntersectionSolid* intersec_iris = new G4IntersectionSolid("box*iris", box, iris_cyl);
395        G4LogicalVolume* iris_log= new G4LogicalVolume(intersec_iris,fe_mat,"iris",0,0,0);
396
397        /*G4VPhysicalVolume* iris_phys = */
398        new G4PVPlacement(0,G4ThreeVector(0,0,iris_dz+iris_thickness),iris_log,"iris",expHall_log,false,0);
399    }
400    else Msg(EsafMsg::Warning) << "IRIS is turned off." << MsgDispatch;
401
402    if(mp1)delete mp1;
403    if(mp2)delete mp2;
404
405// //     // Export VGM geometry to Root
406//     Geant4GM::Factory g4Factory;
407// //     g4Factory.SetDebug(1);
408//     g4Factory.SetIgnore(1);
409//     g4Factory.Import(expHall_phys);
410//
411//     RootGM::Factory rtFactory;
412// //     rtFactory.SetDebug(1);
413//     g4Factory.Export(&rtFactory);
414//     gGeoManager->CloseGeometry();
415//     gGeoManager->Export("output/G4lens.root");
416
417    return expHall_phys;
418}
419
420void G4EusoGeometryShape::SetMaterial(int i, MaterialParameters* mp, char* material){
421
422    char* matname = fOpticalSystemParam->material[i].name;
423    const char* fmatname;
424    if(!matname){
425        matname = fOpticalSystemParam->material[i-1].name;
426        fmatname = Form("%s/%s",lens_dir,fOpticalSystemParam->material[i-1].filename);
427    }
428    else fmatname = Form("%s/%s",lens_dir,fOpticalSystemParam->material[i].filename);
429#ifdef DEBUG
430    printf("%s %s %s %s\n",fOpticalSystemParam->material[0].name,fOpticalSystemParam->material[1].name,
431                           fOpticalSystemParam->material[2].name,fOpticalSystemParam->material[3].name);
432    printf("%s %s\n",fOpticalSystemParam->material[2].filename,fOpticalSystemParam->material[3].filename);
433#endif
434
435
436    if(!strcmp(matname,"CYTOP")){
437        mp->setRI(fmatname);
438        mp->SetData(CYTOP);
439        mp->SetData(fOpticalSurfaceC);
440        strcpy (material,"CYTOP");
441    }
442    else if(!strcmp(matname,"PMMA")){
443           if(fit){
444            TGraph* ph = new TGraph(fmatname);
445            TF1* f=new TF1("ri","1./(x+[0])*[1]+1./(x*x+[2])*[3]+[4]",200.,600.);
446            ph->Fit(f,"QNO","",250.,500.);//NOS
447            mp->setRealRI(f,250.,500.);
448            mp->setTrans(Form("%s/pmma_transmittivity.dat",lens_dir));
449            delete ph;
450            delete f;
451
452        }
453        else{
454            mp->setRI(fmatname);
455        }
456        mp->SetData(PMMA);
457        mp->SetData(fOpticalSurfaceP);
458        strcpy (material,"PMMA");
459
460    }
461    else cerr<<"mat table: material not found"<<endl;
462
463
464}
Note: See TracBrowser for help on using the repository browser.