source: JEM-EUSO/esaf_cc_at_lal/packages/simulation/detector/G4Detector/optics/src/G4EusoGeometryShape.cc @ 114

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

actual version of ESAF at CCin2p3
File size: 19.2 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=true;
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//     downsurf1->draw(lens_file, 0, surf1, "S1", "RECREATE");
261//     upsurf1->draw(lens_file, 0, surf1, "S2", "UPDATE");
262    G4IntersectionSolid* intersec1 = new G4IntersectionSolid("box*FresnelLens_1", box, lens1_solid);
263    G4LogicalVolume* lens1_log = new G4LogicalVolume(intersec1,!strcmp(mat,"CYTOP")?CYTOP:PMMA, "FresnelLens_1",0,0,0);
264   
265//     lens1_solid->GetCylinderLimits(c_rad, c_z1, c_z2);
266//     c_dz=(fabs(c_z2)+fabs(c_z1));
267//     G4Tubs* lens1_mother_solid=new G4Tubs("lens1_mother",0., c_rad+10.*mm, c_dz+10.*mm, 0.*deg, 360.*deg);
268//     G4LogicalVolume* lens1_mother_log = new G4LogicalVolume(lens1_mother_solid, Vacuum, "FresnelLens_1_mother",0,0,0);
269//     G4VPhysicalVolume* lens1_mother_phys = new G4PVPlacement(0, G4ThreeVector(0.,0.,surf1), lens1_mother_log, "FresnelLens_1_mother", expHall_log,false,0);
270
271    #ifndef DISABLE_1_LENSphys
272    G4VPhysicalVolume* lens1_phys = new G4PVPlacement(0, G4ThreeVector(0.,0.,surf1), lens1_log, "FresnelLens_1", expHall_log,false,0);
273    /*G4LogicalBorderSurface diffr_surf =*/
274    new G4LogicalBorderSurface("diffractor_surf1",lens1_phys,expHall_phys,!strcmp(mat,"CYTOP")?fOpticalSurfaceC:fOpticalSurfaceP);
275    #endif
276#endif
277
278    ///2-d lens
279#ifndef DISABLE_2_LENS
280
281    surf1 = fOpticalSystemParam->surface[4].Sz;
282    surf2 = fOpticalSystemParam->surface[5].Sz;
283    thin = surf2 -surf1;
284    Rc1 = fOpticalSystemParam->surface[4].Rc;
285    Rc2 = fOpticalSystemParam->surface[5].Rc;
286    mat = fOpticalSystemParam->material[fOpticalSystemParam->surface[5].matid].name;
287    cut =fOpticalSystemParam->surface[4].r_cut;
288
289    if(Rc2>1.0e10){
290         uplimit=new FlatSurface();
291         downlimit=new FlatSurface();
292    }
293    else{
294        uplimit=new SphericSurface(Rc1*mm,0.,r_wall);
295        downlimit=new SphericSurface(Rc1*mm,0.,r_wall);
296    }
297    shelper->SetOptions(0., -1., -filter);
298    shelper->Read(Form("%s/S4_%s.dat",lens_dir,name));
299    FresnelSurface* downsurf2=shelper->CreateFresnelLens(downlimit, uplimit);
300    double rsurf5=downsurf2->GetR();
301
302    ParametricSurface* upsurf2;
303
304    if(Rc2>1.0e10){
305         upsurf2 = new FlatSurface((thin)*mm,rsurf5*mm);
306    }
307    else{
308        upsurf2=new SphericSurface(Rc2*mm,(thin+Rc2)*mm,rsurf5*mm);
309    }
310    ShapeOfLens* lens2_solid = new ShapeOfLens("FresnelLens_2", (thin*0.5+0.5)*mm, (rsurf5+1.)*mm, upsurf2, downsurf2);
311    G4IntersectionSolid* intersec2 = new G4IntersectionSolid("box*FresnelLens_2", box, lens2_solid);
312   
313    G4LogicalVolume* lens2_log = new G4LogicalVolume(intersec2,!strcmp(mat,"CYTOP")?CYTOP:PMMA, "FresnelLens_2",0,0,0);
314
315//     downsurf2->draw(lens_file, 0, surf1, "S4", "UPDATE");
316//     upsurf2->draw(lens_file, 0, surf1, "S5", "UPDATE");
317
318//     lens2_solid->GetCylinderLimits(c_rad, c_z1, c_z2);
319//     c_dz=(fabs(c_z2)+fabs(c_z1));
320//     G4Tubs* lens2_mother_solid=new G4Tubs("lens2_mother",0., c_rad+10.*mm, c_dz+10.*mm, 0.*deg, 360.*deg);
321//     G4LogicalVolume* lens2_mother_log = new G4LogicalVolume(lens2_mother_solid, Vacuum, "FresnelLens_2_mother",0,0,0);
322//     G4VPhysicalVolume* lens2_mother_phys = new G4PVPlacement(0, G4ThreeVector(0.,0.,surf1), lens2_mother_log, "FresnelLens_2_mother", expHall_log,false,0);
323
324    double diffr_rad=rsurf5; //FIXME
325    double diffr_dz=1.e-8;
326    G4VSolid* euso_diffractor;
327
328    if(Rc2>1.0e10){
329         euso_diffractor=new G4Tubs("euso_diffractor_tubs", 0., diffr_rad, diffr_dz,0.*deg,360.*deg);
330    }
331    else {
332        ParametricSurface* upsurf_diffr = new SphericSurface(Rc2*mm,Rc2*mm+diffr_dz,rsurf5*mm);
333        ParametricSurface* downsurf_diffr = new SphericSurface(Rc2*mm,Rc2*mm-diffr_dz,rsurf5*mm);
334//         downsurf_diffr->draw(lens_file, 0, surf1+thin+diffr_dz+fGeomTolerance, "D1", "UPDATE");
335//         upsurf_diffr->draw(lens_file, 0, surf1+thin+diffr_dz+fGeomTolerance, "D2", "UPDATE");
336        euso_diffractor=new ShapeOfLens("euso_diffractor_par", diffr_dz, diffr_rad, upsurf_diffr, downsurf_diffr);
337       
338    }
339    G4IntersectionSolid* intersec_diffr = new G4IntersectionSolid("box*euso_diffractor_par", box, euso_diffractor);
340    G4LogicalVolume* euso_diffractor_log = new G4LogicalVolume(intersec_diffr,Vacuum2,"euso_diffractor",0,0,0);
341
342    #ifndef DISABLE_2_LENSphys
343    G4VPhysicalVolume* lens2_phys = new G4PVPlacement(0, G4ThreeVector(0.,0.,surf1), lens2_log, "FresnelLens_2", expHall_log,false,0);
344    //G4VPhysicalVolume* euso_diffractor_phys =
345            new G4PVPlacement(0,G4ThreeVector(0,0,thin+diffr_dz+fGeomTolerance+surf1),
346                                              euso_diffractor_log,"euso_diffractor",expHall_log,false,0);
347
348    /*G4LogicalBorderSurface diffr_surf =*/
349    new G4LogicalBorderSurface("diffractor_surf21",lens2_phys,expHall_phys,!strcmp(mat,"CYTOP")?fOpticalSurfaceC:fOpticalSurfaceP);
350    #endif
351#endif
352
353    /// 3-d Lens
354#ifndef DISABLE_3_LENS
355    surf1 = fOpticalSystemParam->surface[6].Sz;
356    surf2 = fOpticalSystemParam->surface[7].Sz;
357    thin = surf2 -surf1;
358    Rc1 = fOpticalSystemParam->surface[6].Rc;
359    Rc2 = fOpticalSystemParam->surface[7].Rc;
360    mat = fOpticalSystemParam->material[fOpticalSystemParam->surface[7].matid].name;
361    cut =fOpticalSystemParam->surface[6].r_cut;
362
363     uplimit=new SphericSurface(Rc1,0.,r_wall);
364     downlimit=new SphericSurface(Rc1,0.,r_wall);
365     shelper->SetOptions(0., -1., -filter);
366     shelper->Read(Form("%s/S6_%s.dat",lens_dir,name));
367     FresnelSurface* downsurf3=shelper->CreateFresnelLens(downlimit, uplimit);
368
369     uplimit=new SphericSurface(Rc2,0.,r_wall);
370     downlimit=new SphericSurface(Rc2,0.,r_wall);
371     shelper->SetOptions(thin*mm, -1., -filter);
372     shelper->Read(Form("%s/S7_%s.dat",lens_dir,name));
373     FresnelSurface* upsurf3=shelper->CreateFresnelLens(downlimit, uplimit);
374
375     ShapeOfLens* lens3_solid = new ShapeOfLens("FresnelLens_3",(thin*0.5+0.5)*mm, (r_wall+1.)*mm, upsurf3, downsurf3);
376     G4IntersectionSolid* intersec3 = new G4IntersectionSolid("box*FresnelLens_3", box, lens3_solid);
377//      downsurf3->draw(lens_file, 0, surf1, "S6", "UPDATE");
378//      upsurf3->draw(lens_file, 0, surf1, "S7", "UPDATE");
379
380     G4LogicalVolume* lens3_log = new G4LogicalVolume(intersec3,!strcmp(mat,"CYTOP")?CYTOP:PMMA, "FresnelLens_3",0,0,0);
381
382     lens3_solid->GetCylinderLimits(c_rad, c_z1, c_z2);
383     c_dz=(fabs(c_z2)+fabs(c_z1));
384
385
386    #ifndef DISABLE_3_LENSphys
387    G4VPhysicalVolume* lens3_phys = new G4PVPlacement(0, G4ThreeVector(0.,0.,surf1), lens3_log, "FresnelLens_3", expHall_log,false,0);
388    /*G4LogicalBorderSurface diffr_surf =*/
389    new G4LogicalBorderSurface("diffractor_surf3",lens3_phys,expHall_phys,!strcmp(mat,"CYTOP")?fOpticalSurfaceC:fOpticalSurfaceP);
390    #endif
391#endif
392    delete shelper;
393
394    /// fs
395    double fs_dz=10.*mm;
396    double fs_pos= ((fOpticalSystemParam->FS_OFFSET-surf2)*0.5+surf2)*mm;//surf1+c_dz+15.*mm+fs_dz;
397   // printf("fpos************************ %g",fs_pos);
398    double fs_rad=fOpticalSystemParam->r_fs*mm;
399
400    G4Tubs* absorber_cyl=new G4Tubs("absorber",0.,fs_rad,fs_dz,0.*deg,360.*deg);
401    G4IntersectionSolid* intersec_fs = new G4IntersectionSolid("box*absorber", box, absorber_cyl);
402    G4LogicalVolume* absorber_cyl_log = new G4LogicalVolume(intersec_fs,Vacuum1,"absorber",0,0,0);
403    //G4VPhysicalVolume* absorber_cyl_phys =
404            new G4PVPlacement(0,G4ThreeVector(0,0,fs_pos),absorber_cyl_log,"absorber",expHall_log,false,0);
405
406
407    ///IRIS
408    bool useiris=Conf()->GetBool("G4EusoGeometryShape.UseIRIS");
409    if ( useiris ){
410        double iris_r_inner=fOpticalSystemParam->surface[3].r_lens;
411        double iris_r_outer=fOpticalSystemParam->r_wall;
412        double iris_dz=fOpticalSystemParam->surface[3].Sz;
413        //double iris_cut=fOpticalSystemParam->surface[3].r_cut;
414        double iris_thickness=1.*mm;
415        G4Tubs* iris_cyl=new G4Tubs("iris",iris_r_inner,iris_r_outer,iris_thickness,0.*deg,360.*deg);
416        G4IntersectionSolid* intersec_iris = new G4IntersectionSolid("box*iris", box, iris_cyl);
417        G4LogicalVolume* iris_log= new G4LogicalVolume(intersec_iris,fe_mat,"iris",0,0,0);
418
419        /*G4VPhysicalVolume* iris_phys = */
420        new G4PVPlacement(0,G4ThreeVector(0,0,iris_dz+iris_thickness),iris_log,"irisr",expHall_log,false,0);
421    }
422    else Msg(EsafMsg::Warning) << "IRIS is turned off." << MsgDispatch;
423
424    if(mp1)delete mp1;
425    if(mp2)delete mp2;
426
427// //     // Export VGM geometry to Root
428//     Geant4GM::Factory g4Factory;
429// //     g4Factory.SetDebug(1);
430//     g4Factory.SetIgnore(1);
431//     g4Factory.Import(expHall_phys);
432//
433//     RootGM::Factory rtFactory;
434// //     rtFactory.SetDebug(1);
435//     g4Factory.Export(&rtFactory);
436//     gGeoManager->CloseGeometry();
437//     gGeoManager->Export("output/G4lens.root");
438
439    return expHall_phys;
440}
441
442void G4EusoGeometryShape::SetMaterial(int i, MaterialParameters* mp, char* material){
443
444    char* matname = fOpticalSystemParam->material[i].name;
445    const char* fmatname;
446    if(!matname){
447        matname = fOpticalSystemParam->material[i-1].name;
448        fmatname = Form("%s/%s",lens_dir,fOpticalSystemParam->material[i-1].filename);
449    }
450    else fmatname = Form("%s/%s",lens_dir,fOpticalSystemParam->material[i].filename);
451#ifdef DEBUG
452    printf("%s %s %s %s\n",fOpticalSystemParam->material[0].name,fOpticalSystemParam->material[1].name,
453                           fOpticalSystemParam->material[2].name,fOpticalSystemParam->material[3].name);
454    printf("%s %s\n",fOpticalSystemParam->material[2].filename,fOpticalSystemParam->material[3].filename);
455#endif
456
457
458    if(!strcmp(matname,"CYTOP")){
459        mp->setRI(fmatname);
460        mp->SetData(CYTOP);
461        mp->SetData(fOpticalSurfaceC);
462        strcpy (material,"CYTOP");
463    }
464    else if(!strcmp(matname,"PMMA")){
465           if(fit){
466            TGraph* ph = new TGraph(fmatname);
467            TF1* f=new TF1("ri","1./(x+[0])*[1]+1./(x*x+[2])*[3]+[4]",200.,600.);
468            //printf("mp2 fit\n");
469            ph->Fit(f,"QNO","",250.,500.);//NOS
470            //TFile ff("fit.root","RECREATE");
471            //ph->Draw("A*");
472            //ph->Write("fit");
473            //ff.Close();
474            //printf("mp2 fit\n");
475            mp->setRealRI(f,250.,500.);
476            mp->setTrans(Form("%s/pmma_transmittivity.dat",lens_dir));
477            delete ph;
478            delete f;
479
480        }
481        else{
482            mp->setRI(fmatname);
483        }
484        mp->SetData(PMMA);
485        mp->SetData(fOpticalSurfaceP);
486        strcpy (material,"PMMA");
487
488    }
489    else cerr<<"mat table: material not found"<<endl;
490
491
492}
Note: See TracBrowser for help on using the repository browser.