Changeset 184 in JEM-EUSO


Ignore:
Timestamp:
Jul 15, 2013, 2:08:33 PM (11 years ago)
Author:
moretto
Message:

changes from biktem r:3032

Location:
esaf_lal/branches/camille/packages
Files:
21 edited

Legend:

Unmodified
Added
Removed

  • esaf_lal/branches/camille/packages/g4libs.gmk

    r16 r184  
    11#this config file includes Geant4 libs
     2
     3ifneq ($(shell which geant4-config 2>/dev/null),)
     4    G4LIBS := $(shell geant4-config --libs)
     5    G4INCLUDES := $(shell geant4-config --cflags)
     6    USE_G4 = 1
     7    export USE_G4 G4LIBS G4INCLUDES
     8else
    29ifdef G4INSTALL
    3 GEANT_LIBDIR := $(G4INSTALL)/lib/Linux-g++
    4 G4LIBS   := -L${GEANT_LIBDIR} $(shell cat ${GEANT_LIBDIR}/libname.map | ${GEANT_LIBDIR}/liblist -m ${GEANT_LIBDIR})
    5 G4LIBS   += -L${CLHEP_LIB_DIR} -l${CLHEP_LIB}
    6 #-L/usr/lib -lxerces-c-3.1
    7 export G4LIBS
     10    GEANT_LIBDIR := $(G4INSTALL)/lib/Linux-g++
     11    G4LIBS   := -L${GEANT_LIBDIR} $(shell cat ${GEANT_LIBDIR}/libname.map | ${GEANT_LIBDIR}/liblist -m ${GEANT_LIBDIR})
     12    G4LIBS   += -L${CLHEP_LIB_DIR} -l${CLHEP_LIB}
     13
     14    G4INCLUDES += -I$(G4INSTALL)/include
     15    G4INCLUDES += -I$(CLHEP_INCLUDE_DIR)
     16
     17    USE_G4 = 1
     18
     19ifdef VGM_INSTALL
     20    VGMINCLUDES := -I$(VGM_INSTALL)/packages/VGM/include
     21    VGMINCLUDES += -I$(VGM_INSTALL)/packages/BaseVGM/include
     22    VGMINCLUDES += -I$(VGM_INSTALL)/packages/Geant4GM/include
     23    VGMINCLUDES += -I$(VGM_INSTALL)/packages/RootGM/include
     24    export VGMINCLUDES
    825endif
     26
     27    export USE_G4 G4LIBS G4INCLUDES
     28endif
     29endif

  • esaf_lal/branches/camille/packages/simulation/detector/G4Detector/G4fresnellens/GNUmakefile

    r16 r184  
    1717
    1818include $(ESAFINSTALL)/packages/config.gmk
    19 #include $(G4INSTALL)/config/G4VIS_USE.gmk
    2019
    2120# Include paths needed for rootcint dictionary generation
    … …  
    2423INCLUDES += -I$(ESAFPACKAGES)/common/root/include
    2524INCLUDES += -I$(ESAFPACKAGES)/simulation/tools/include
    26 INCLUDES += -I$(G4INSTALL)/include
    27 INCLUDES += -I$(CLHEP_INCLUDE_DIR)
     25INCLUDES += $(G4INCLUDES)
    2826
    2927# headers and settings needed for the local library compilation

  • esaf_lal/branches/camille/packages/simulation/detector/G4Detector/G4fresnellens/include/FresnelSurface.h

    r16 r184  
    2121
    2222    void setSmoothedEdges(double rad);
    23     double PointOnSurf(double x) { AbstractMethod(); return kInfinity; }
     23    double PointOnSurf(double x);
    2424    double DistanceToSurf(const G4ThreeVector& p, const G4ThreeVector& v){ AbstractMethod(); return kInfinity; }
    2525    double DistanceToSurf(const G4ThreeVector& p, const G4ThreeVector& v, FSIntersectedSegment*){ AbstractMethod(); return kInfinity; }
    … …  
    4747    int     GetNTeeth()     { return theToothN-1; }
    4848    std::vector<SSurface*>& GetTeeth() { return fSurfaces; }
     49   
    4950protected:
    5051    #ifdef USE_ROOT
    … …  
    5354    void FillArray(int n, double* rho, double* z);
    5455    void FillArray(const std::vector<SSurface*>& surf);
     56    int FindTooth( double rho );
    5557
    5658    /// Surfaces

  • esaf_lal/branches/camille/packages/simulation/detector/G4Detector/G4fresnellens/include/SSurface.h

    r16 r184  
    99
    1010#include <geomdefs.hh>
     11#include <assert.h>
    1112#include <G4ThreeVector.hh>
    1213class G4PhysicalVolume;
    … …  
    5859
    5960    // plug for temporarily unimplemented methods
    60     void AbstractMethod() { fprintf(stderr, "Error, calling the unimplemented method"); exit(1); }
     61    void AbstractMethod() { fprintf(stderr, "Error, calling the unimplemented method. The method was not supposed to be used."); assert(false); }
    6162protected:
    6263    double fR1;

  • esaf_lal/branches/camille/packages/simulation/detector/G4Detector/G4fresnellens/include/ShapeOfLens.h

    r16 r184  
    6363
    6464    void DescribeYourselfTo(G4VGraphicsScene& scene) const;
    65 /*    G4Polyhedron* CreatePolyhedron() const;
     65    G4Polyhedron* CreatePolyhedron() const;
    6666    G4Polyhedron* GetPolyhedron () const;
    67     G4NURBS*      CreateNURBS() const;*/
     67    G4NURBS*      CreateNURBS() const;
    6868
    6969    void GetCylinderLimits(double& r, double& z1, double& z2) { r=fR; z1=fZ1; z2=fZ2; }

  • esaf_lal/branches/camille/packages/simulation/detector/G4Detector/G4fresnellens/src/FresnelSurface.cc

    r16 r184  
    4646    }
    4747    fSurfacesIns.clear();
     48}
     49
     50//__________________________________________________________________________________
     51int FresnelSurface::FindTooth( double rho ) {
     52    return BinarySearch(theToothEdges, rho, 0, theToothN);
     53}
     54
     55//__________________________________________________________________________________
     56double FresnelSurface::PointOnSurf(double rho) {
     57    int n_tooth = FindTooth( rho );
     58    if( n_tooth<0 ) return kInfinity;
     59    return fSurfacesIns[ n_tooth*2 ]->PointOnSurf( rho );   
    4860}
    4961
    … …  
    180192      }
    181193    #endif
    182     int n_tooth = BinarySearch(theToothEdges,rho,0,theToothN);
     194    int n_tooth = FindTooth( rho );
    183195    if( n_tooth<0 ) {
    184196        printf("Warning. Tooth not found. Shouldn't happen (rho=%g of [%g, %g]).",
    … …  
    255267    }// approach to the limiting surface
    256268
    257     int theCurrentToothIndex=BinarySearch(theToothEdges, perp, 0, theToothN);
     269    int theCurrentToothIndex=FindTooth( perp );
    258270    if ( theCurrentToothIndex<0 ) return kInfinity;
    259271    vector<SSurface*>::const_iterator it=fSurfaces.begin()+theCurrentToothIndex*2;

  • esaf_lal/branches/camille/packages/simulation/detector/G4Detector/G4fresnellens/src/ShapeOfLens.cc

    r16 r184  
    463463                                   GetTrack()->GetTrackID();
    464464            printf("Event %i, %s\n", id, GetName().data());
    465             printf("pos(%.10g,%.10g,%.10g)\n", p.x(), p.y(), p.z());
    466465            double rho=p.perp();
    467             printf("dz1 %g, dz2 %g\n", p.z()-fUpSurface->PointOnSurf(rho),p.z()-fDownSurface->PointOnSurf(rho));
     466            printf("pos (%.10g,%.10g,%.10g), rho=%.10g\n", p.x(), p.y(), p.z(), rho);
     467            double z1=fUpSurface->PointOnSurf(rho), z2=fDownSurface->PointOnSurf(rho);
     468            printf("dz1 %g, dz2 %g: z1, z2 (%g, %g)\n", p.z()-z1,p.z()-z2, z1, z2);
    468469            return norm;
    469470        }
    … …  
    518519    top_volume->DescribeYourselfTo(scene);
    519520}
     521   
     522//__________________________________________________________________________________
     523G4Polyhedron* ShapeOfLens::CreatePolyhedron() const {
     524  return top_volume->CreatePolyhedron();
     525}
     526
     527//__________________________________________________________________________________
     528G4Polyhedron* ShapeOfLens::GetPolyhedron () const {
     529  return top_volume->GetPolyhedron();
     530}
     531
     532//__________________________________________________________________________________
     533G4NURBS*      ShapeOfLens::CreateNURBS() const {
     534  return top_volume->CreateNURBS();
     535}
     536
    520537//__________________________________________________________________________________
    521538G4ThreeVector ShapeOfLens::GetPointOnSurface() const{

  • esaf_lal/branches/camille/packages/simulation/detector/G4Detector/GNUmakefile

    r16 r184  
    11
    2 ifdef G4INSTALL
    32SUBDIRS = optics G4fresnellens
    4 endif
    53
    64.PHONY: clean cleandict obj lib tags doc all

  • esaf_lal/branches/camille/packages/simulation/detector/G4Detector/optics/GNUmakefile

    r16 r184  
    3535INCLUDES += -I$(G4INSTALL)/include
    3636INCLUDES += -I$(CLHEP_INCLUDE_DIR)
    37 ifdef VGM_INSTALL
    38 INCLUDES += -I$(VGM_INSTALL)/packages/VGM/include
    39 INCLUDES += -I$(VGM_INSTALL)/packages/BaseVGM/include
    40 INCLUDES += -I$(VGM_INSTALL)/packages/Geant4GM/include
    41 INCLUDES += -I$(VGM_INSTALL)/packages/RootGM/include
     37INCLUDES += $(G4INCLUDES)
     38ifdef VGMINCLUDES
     39    INCLUDES += $(VGMINCLUDES)
    4240endif
     41#
     42#
    4343# headers and settings needed for the local library compilation
    4444CXXFLAGS += $(INCLUDES) $(CPPFLAGS)

  • esaf_lal/branches/camille/packages/simulation/detector/G4Detector/optics/include/OptPrimaryGeneratorAction.hh

    r16 r184  
    2323        void Init();
    2424        void FillCircle(G4Event*);
     25        void FillSphere(G4Event*);
    2526        void OnePoint(G4Event* anEvent);
    2627        void GenerateEnergy();

  • esaf_lal/branches/camille/packages/simulation/detector/G4Detector/optics/include/OptSteppingAction.hh

    r16 r184  
    2626    G4OpBoundaryProcessStatus   theStatus;
    2727    int                         fStatus;
     28    int                         fStatusKilled;
    2829    OptUserStackingAction*      fStacking;
    2930

  • esaf_lal/branches/camille/packages/simulation/detector/G4Detector/optics/include/OptUserStackingAction.hh

    r16 r184  
    2121class G4TrackInformation;
    2222class ListPhotonsOnPupil;
     23
     24class TGraph2D;
     25
    2326class OptUserStackingAction : public G4UserStackingAction, public G4UserTrackingAction, public G4UserEventAction, public EsafMsgSource, public VirtualTelParm  {
    2427public:
    … …  
    3538    void PrepareNewEvent();
    3639    void SetPupil(PhotonsOnPupil*);
     40 
    3741    void SetNmax(int m) { fNmax=m; }
    3842    void Progress(int, const char*);
    … …  
    4448    void SetTrackHistory(int hist){fTrackHistory=hist;}
    4549    OptSteppingAction* stepping;
    46 
     50   
     51    G4Track* NewPhoton(Photon*,int,int n=0);
     52    void     Reset();
    4753private:
    4854
    4955    bool GoToElectronics(Photon*);
    50     void     Reset();
     56   
    5157    EVector  NewInteractionPoint(const Photon *p, double DZ) const;
    52     G4Track* NewPhoton(Photon*,int,int n=0);
     58   
    5359    void CheckTrace(int id,const G4ThreeVector& pos,const G4ThreeVector& mom,double lambda,double time);
    5460
    … …  
    6470    double   fTimeFirstPhoton;
    6571    double   fTimeLastPhoton;
    66 
     72   
     73    int     fNReflElec;
     74   
    6775    int     fNTracked;
    6876    int     fNG4Tracked;
    … …  
    7078    int     fNTotal;
    7179    int     fNG4focalplane;
     80    int     fNiris;
    7281    int     fNfocalplane;
    7382    int     fNlost;
    7483    int     fNoptics;
     84                int     fNwalls;
     85                int     fNearth;
    7586    int     fProgress;
    7687    int     fTrackHistory;
    7788    int     forigin;
    7889    bool    ftus;
     90   
     91    TGraph2D* f_before;
     92    TGraph2D* f_after;
    7993    ListPhotonsOnPupil* fPhOnPupil;
    8094};

  • esaf_lal/branches/camille/packages/simulation/detector/G4Detector/optics/src/G4Detector.cc

    r16 r184  
    1212#include "OpticsFactory.hh"
    1313#include "Pupil.hh"
    14 
     14#include "ParentPhoton.hh"
     15#include "Photon.hh"
    1516
    1617#include <TMath.h>
    17 
     18#include <TFile.h>
    1819using namespace sou;
    1920
    … …  
    101102
    102103    // wrapper to track the photons onto the real pupil
     104   
     105//     TFile f("output/event.root","RECREATE");
     106   
    103107    Pupil pupil( photons, GetGeometry());
     108/*             
     109                TTree* t = new TTree("t","tt");
     110                Photon* photon =0;
     111                ParentPhoton* pp =0;*/
    104112
     113//              TVector3 vpos(0.,0.,0.);
     114//              TVector3 vdir(0.,0.,0.);
     115//             
     116//              t->Branch("position","TVector3",&vpos);
     117//              t->Branch("direction","TVector3",&vdir);
     118//              t->Branch("wl",&wl,"w/D");
     119//              t->Branch("flag",&fl,"flag/I");
     120
     121//              t->Branch("photon","ParentPhoton",&pp);
    105122    if(pupil.GetNphotons()==0) Msg(EsafMsg::Warning)  << "NO PHOTONS ON PUPIL !!!"<< MsgDispatch;
    106     else Msg(EsafMsg::Info)  << "PHOTONS ON PUPIL: "<<pupil.GetNphotons()<< MsgDispatch;
     123    else {Msg(EsafMsg::Info)  << "PHOTONS ON PUPIL: "<<pupil.GetNphotons()<< MsgDispatch;
     124//                      photon = pupil.GetPhoton();
     125//                      pp=photon->parent;
     126//                      printf("www  %g\n",pp->W());
     127//                      int id=0;
     128//                     
     129//                while(id<pupil.GetNphotons()){
     130//                              photon->id = id;
     131//                              t->Fill();
     132//                             
     133//                              photon = pupil.GetPhoton();
     134//                              if(photon){
     135//                                      pp=photon->parent;
     136//                                      printf("www  %g\n",pp->W());
     137//                                     
     138//                              }
     139//                              id++;
     140//      //                      photon->Dump();
     141//                     
     142//                      }
     143//                      f.Write();
     144//                      f.Close();
     145                }
    107146
    108147    // transport photons from pupil to pmts

  • esaf_lal/branches/camille/packages/simulation/detector/G4Detector/optics/src/G4EusoGeometryShape.cc

    r16 r184  
    6464    expHall_x = expHall_y = 2000.*mm;
    6565    expHall_z = 4000.*mm;
    66     fit=true;
     66    fit = Conf()->GetBool("G4EusoGeometryShape.UseFitTransmittivity");
    6767    CYTOP=0;
    6868    PMMA=0;
    … …  
    258258    ShapeOfLens* lens1_solid = new ShapeOfLens("FresnelLens_1", (thin*0.5+0.5)*mm, (r_wall+1.)*mm, upsurf1, downsurf1);
    259259
    260 //     downsurf1->draw(lens_file, 0, surf1, "S1", "RECREATE");
    261 //     upsurf1->draw(lens_file, 0, surf1, "S2", "UPDATE");
    262260    G4IntersectionSolid* intersec1 = new G4IntersectionSolid("box*FresnelLens_1", box, lens1_solid);
    263261    G4LogicalVolume* lens1_log = new G4LogicalVolume(intersec1,!strcmp(mat,"CYTOP")?CYTOP:PMMA, "FresnelLens_1",0,0,0);
    264262   
    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 
    271263    #ifndef DISABLE_1_LENSphys
    272264    G4VPhysicalVolume* lens1_phys = new G4PVPlacement(0, G4ThreeVector(0.,0.,surf1), lens1_log, "FresnelLens_1", expHall_log,false,0);
    … …  
    313305    G4LogicalVolume* lens2_log = new G4LogicalVolume(intersec2,!strcmp(mat,"CYTOP")?CYTOP:PMMA, "FresnelLens_2",0,0,0);
    314306
    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
     307    double diffr_rad=rsurf5; //FIXME: diffractive surface is shorter than surf5
    325308    double diffr_dz=1.e-8;
    326309    G4VSolid* euso_diffractor;
    … …  
    332315        ParametricSurface* upsurf_diffr = new SphericSurface(Rc2*mm,Rc2*mm+diffr_dz,rsurf5*mm);
    333316        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");
    336317        euso_diffractor=new ShapeOfLens("euso_diffractor_par", diffr_dz, diffr_rad, upsurf_diffr, downsurf_diffr);
    337318       
    … …  
    375356     ShapeOfLens* lens3_solid = new ShapeOfLens("FresnelLens_3",(thin*0.5+0.5)*mm, (r_wall+1.)*mm, upsurf3, downsurf3);
    376357     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");
    379358
    380359     G4LogicalVolume* lens3_log = new G4LogicalVolume(intersec3,!strcmp(mat,"CYTOP")?CYTOP:PMMA, "FresnelLens_3",0,0,0);
    … …  
    395374    double fs_dz=10.*mm;
    396375    double fs_pos= ((fOpticalSystemParam->FS_OFFSET-surf2)*0.5+surf2)*mm;//surf1+c_dz+15.*mm+fs_dz;
    397    // printf("fpos************************ %g",fs_pos);
    398376    double fs_rad=fOpticalSystemParam->r_fs*mm;
    399377
    … …  
    418396
    419397        /*G4VPhysicalVolume* iris_phys = */
    420         new G4PVPlacement(0,G4ThreeVector(0,0,iris_dz+iris_thickness),iris_log,"irisr",expHall_log,false,0);
     398        new G4PVPlacement(0,G4ThreeVector(0,0,iris_dz+iris_thickness),iris_log,"iris",expHall_log,false,0);
    421399    }
    422400    else Msg(EsafMsg::Warning) << "IRIS is turned off." << MsgDispatch;
    … …  
    466444            TGraph* ph = new TGraph(fmatname);
    467445            TF1* f=new TF1("ri","1./(x+[0])*[1]+1./(x*x+[2])*[3]+[4]",200.,600.);
    468             //printf("mp2 fit\n");
    469446            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");
    475447            mp->setRealRI(f,250.,500.);
    476448            mp->setTrans(Form("%s/pmma_transmittivity.dat",lens_dir));

  • esaf_lal/branches/camille/packages/simulation/detector/G4Detector/optics/src/G4TusGeometry.cc

    r16 r184  
    2424#include <G4Polyhedra.hh>
    2525
    26 #include <Geant4GM/volumes/Factory.h>
    27 #include <RootGM/volumes/Factory.h>
     26//#include <Geant4GM/volumes/Factory.h>
     27//#include <RootGM/volumes/Factory.h>
    2828
    2929

  • esaf_lal/branches/camille/packages/simulation/detector/G4Detector/optics/src/OpEUSODiffraction.cc

    r16 r184  
    88#include <G4TransportationManager.hh>
    99#include <G4Navigator.hh>
     10#include <G4Version.hh>
    1011
    1112#include <TGraph.h>
    … …  
    9394    }
    9495    if ( Rindex ) {
    95         Rindex1 = Rindex->GetProperty(thePhotonMomentum);
     96                #if G4VERSION_NUMBER<950
     97            Rindex1 = Rindex->GetProperty(thePhotonMomentum);//g9.4
     98                #else
     99                Rindex1 = Rindex->Value(thePhotonMomentum); //g9.5
     100                #endif
    96101    } else {
    97102        theStatus = edNoINDEX;
    … …  
    108113
    109114    if ( Rindex ) {
    110         Rindex2 = Rindex->GetProperty(thePhotonMomentum);
     115                #if G4VERSION_NUMBER<950
     116                Rindex2 = Rindex->GetProperty(thePhotonMomentum);//g9.4
     117                #else
     118            Rindex2 = Rindex->Value(thePhotonMomentum);//g9.5
     119                #endif
    111120    } else {
    112121        theStatus = edNoINDEX;

  • esaf_lal/branches/camille/packages/simulation/detector/G4Detector/optics/src/OptPhysicsList.cc

    r16 r184  
    115115{
    116116  AddTransportation();
    117    ConstructGeneral();
    118    ConstructEM();
     117//    ConstructGeneral();
     118  ConstructEM();
    119119  ConstructOp();
    120120}
    … …  
    136136    if (theDecayProcess->IsApplicable(*particle)) {
    137137        pmanager->AddDiscreteProcess(theDecayProcess);
    138 //      pmanager ->AddProcess(theDecayProcess);
    139 //       // set ordering for PostStepDoIt and AtRestDoIt
    140 //       pmanager ->SetProcessOrdering(theDecayProcess, idxPostStep);
    141 //       pmanager ->SetProcessOrdering(theDecayProcess, idxAtRest);
     138     pmanager ->AddProcess(theDecayProcess);
     139      // set ordering for PostStepDoIt and AtRestDoIt
     140      pmanager ->SetProcessOrdering(theDecayProcess, idxPostStep);
     141      pmanager ->SetProcessOrdering(theDecayProcess, idxAtRest);
    142142    }
    143143  }
    … …  
    234234  if ( use_diffraction ) theEUSODiffractionProcess = new OpEUSODiffraction();
    235235  else Msg(EsafMsg::Warning) << "Diffraction is turned off." << MsgDispatch;
    236     theCerenkovProcess           -> SetVerboseLevel(0);
    237     theAbsorptionProcess         -> SetVerboseLevel(0);
    238     theRayleighScatteringProcess -> SetVerboseLevel(0);
    239     theBoundaryProcess           -> SetVerboseLevel(0);
     236//   theCerenkovProcess           -> SetVerboseLevel(0);
     237//   theAbsorptionProcess         -> SetVerboseLevel(0);
     238//   theRayleighScatteringProcess -> SetVerboseLevel(0);
     239//   theBoundaryProcess           -> SetVerboseLevel(0);
    240240//  theCerenkovProcess->DumpPhysicsTable();
    241241//  theScintillationProcess->DumpPhysicsTable();
    … …  
    243243//  theRayleighScatteringProcess->DumpPhysicsTable();
    244244
    245  // theCerenkovProcess->SetMaxNumPhotonsPerStep(300);
    246   theCerenkovProcess->SetTrackSecondariesFirst(true);
    247 /*
    248   theScintillationProcess->SetScintillationYieldFactor(1.);
    249   theScintillationProcess->SetTrackSecondariesFirst(true);*/
    250 
    251 //   G4EmSaturation* emSaturation = G4LossTableManager::Instance()->EmSaturation();
    252 //   theScintillationProcess->AddSaturation(emSaturation);
     245 //theCerenkovProcess->SetMaxNumPhotonsPerStep(300);
     246//    theCerenkovProcess->SetTrackSecondariesFirst(true);
     247
    253248
    254249  G4OpticalSurfaceModel themodel = unified;
    … …  
    265260      pmanager->SetProcessOrdering(theCerenkovProcess,idxPostStep);
    266261    }
    267 //     if (theScintillationProcess->IsApplicable(*particle)) {
    268 //       pmanager->AddProcess(theScintillationProcess);
    269 //       pmanager->SetProcessOrderingToLast(theScintillationProcess, idxAtRest);
    270 //       pmanager->SetProcessOrderingToLast(theScintillationProcess, idxPostStep);
    271 //     }
     262
    272263    if (particleName == "opticalphoton") {
    273264      G4cout << " AddDiscreteProcess to OpticalPhoton " << G4endl;

  • esaf_lal/branches/camille/packages/simulation/detector/G4Detector/optics/src/OptPrimaryGeneratorAction.cc

    r16 r184  
    7979    else if(!FillPosition.compare ("Point"))
    8080        OnePoint(anEvent);
     81    else if(!FillPosition.compare ("Sphere"))
     82        FillSphere(anEvent);
    8183
    8284}
    … …  
    142144
    143145
     146//-----------------------------------------------------------------------------------------
     147void OptPrimaryGeneratorAction::FillSphere(G4Event* anEvent){
     148    TVector3 p;
     149    double phi = gRandom->Uniform(0.,TwoPi());
     150    double th = gRandom->Uniform(Pi()-asin(1325./3143.57795946),Pi());
     151    p.SetMagThetaPhi(3143.57795946,th,phi);
     152    p.SetZ(p.Z()+3143.57795946+78.);
     153    particleGun->SetParticlePosition(G4ThreeVector(p.X(),p.Y(),p.Z()));
     154    particleGun->GeneratePrimaryVertex(anEvent);
     155}
    144156

  • esaf_lal/branches/camille/packages/simulation/detector/G4Detector/optics/src/OptSteppingAction.cc

    r16 r184  
    5353    kWalls                    = 5,
    5454    kEarth                    = 6,
    55     kG4FocalSurface           = 7
     55    kIris                     = 7,
     56    kG4FocalSurface           = 8
    5657};
    5758
    … …  
    5960{
    6061    //theStatus = Undefined;
     62    //double fGeomTolerance=G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
     63
    6164    fBoundary = NULL;
    6265    fProcessManager = NULL;
    6366    fProcessVector = NULL;
    6467    fStatus = 0;
     68    fStatusKilled = 0;
    6569    fStacking = stack;
    6670    ConfigFileParser* conf = Config::Get()->GetCF("G4","G4Detector");
    6771    string det=conf->GetStr("G4Detector.DetectorType");
    68     //double fGeomTolerance=G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
    6972    if (det.find("Tus")!=string::npos){
    70         fTus=true;
    71         cover=0;
     73      fTus=true;
     74      cover=0;
    7275    }
    7376    else {
    … …  
    8790    }
    8891    SetSender("OptSteppingAction");
     92       
    8993
    9094}
    … …  
    9296OptSteppingAction::~OptSteppingAction()
    9397{
    94 
     98       
    9599}
    96100
    97101void OptSteppingAction::UserSteppingAction(const G4Step * theStep)
    98102{
    99     G4VPhysicalVolume* prevol=theStep->GetPreStepPoint()->GetPhysicalVolume();
    100     G4VPhysicalVolume* postvol=theStep->GetPostStepPoint()->GetPhysicalVolume();
    101     G4String PreStepVolume,PostStepVolume;
    102     if(prevol){
    103         PreStepVolume = prevol->GetName();
    104     }
    105     if(postvol){
    106         PostStepVolume = postvol->GetName();
    107     }
    108    
    109     G4StepPoint *PostStep = 0;
    110     G4StepPoint *PreStep = 0;
    111     G4ThreeVector pre_pos, pre_dir;
    112     G4ThreeVector post_pos, post_dir;
    113    
    114     PreStep = theStep->GetPreStepPoint();
    115     pre_pos = PreStep->GetPosition();
    116     pre_dir = PreStep->GetMomentumDirection();
    117     PostStep = theStep->GetPostStepPoint();
    118     post_pos =  PostStep->GetPosition();
    119     post_dir =  PostStep->GetMomentumDirection();
    120 
    121    
    122     fStepCounter++;
    123     if ( fStepCounter>fMaxIterations ){
    124         theStep->GetTrack()->SetTrackStatus(fStopAndKill);
    125         MsgForm(EsafMsg::Warning,"This photon exceeded the maximum number of iterations (%i), killed.",fMaxIterations);
    126     }
     103  G4VPhysicalVolume* prevol=theStep->GetPreStepPoint()->GetPhysicalVolume();
     104  G4VPhysicalVolume* postvol=theStep->GetPostStepPoint()->GetPhysicalVolume();
     105  G4String PreStepVolume,PostStepVolume;
     106  if(prevol){
     107    PreStepVolume = prevol->GetName();
     108  }
     109  if(postvol){
     110    PostStepVolume = postvol->GetName();
     111  }
     112   
     113  G4StepPoint *PostStep = 0;
     114  G4StepPoint *PreStep = 0;
     115  G4ThreeVector pre_pos, pre_dir;
     116  G4ThreeVector post_pos, post_dir;
     117 
     118  PreStep = theStep->GetPreStepPoint();
     119  pre_pos = PreStep->GetPosition();
     120  pre_dir = PreStep->GetMomentumDirection();
     121  PostStep = theStep->GetPostStepPoint();
     122  post_pos =  PostStep->GetPosition();
     123  post_dir =  PostStep->GetMomentumDirection();
     124
     125   
     126  fStepCounter++;
     127  if ( fStepCounter>fMaxIterations ){
     128    theStep->GetTrack()->SetTrackStatus(fStopAndKill);
     129    MsgForm(EsafMsg::Warning,"This photon exceeded the maximum number of iterations (%i), killed.",fMaxIterations);
     130    fStatus=-3;
     131    fStatusKilled++;
     132  }
    127133
    128134    //find the boundary process only once
    129     if( !fBoundary ){
    130         fProcessManager = theStep->GetTrack()->GetDefinition()->GetProcessManager();
    131         G4int nprocesses = fProcessManager->GetProcessListLength();
    132         fProcessVector = fProcessManager->GetProcessList();
    133 
    134         for(int i=0;i<nprocesses;i++){
    135             if( (*fProcessVector)[i]->GetProcessName()=="OpBoundary" ){
    136                 fBoundary = (G4OpBoundaryProcess*)(*fProcessVector)[i];
    137                 break;
    138             }
    139         }
    140     }
    141 
    142     theStatus = Undefined;
    143     if(fBoundary)theStatus = fBoundary->GetStatus();
    144 
    145     double rho = post_pos.perp();
    146     if(!fTus && (fabs(post_pos.y())>=r_cut || rho>=r_wall)){
    147         if(cover->Inside (pre_pos)==kInside ){
    148             theStep->GetTrack()->SetTrackStatus(fStopAndKill);
    149             fStatus=-4;
    150             double dist = cover->DistanceToOut (pre_pos, pre_dir);
    151             // intersection with virtual wall
    152            
    153             if(dist!=kInfinity ){
    154               //printf("dist %g\n",dist);
    155                 pre_pos+=pre_dir.unit()*dist;
    156                 PostStep->SetPosition(pre_pos);
    157                 PostStep->SetMomentumDirection(pre_dir);
    158                 post_pos=pre_pos;
    159                 post_dir=pre_dir;
    160                 double t = PostStep->GetLocalTime();
    161                 t+=dist/c_light;
    162                 PostStep->SetLocalTime(t);
    163             }
    164         }
    165         else if(cover->Inside (pre_pos)==kSurface){
    166            
    167             //double newdist= cover->DistanceToIn (pre_pos, pre_dir);
    168             //if(newdist==kInfinity ){
    169               PostStep->SetPosition(pre_pos);
    170               //printf("dist %g\n",newdist);
    171             //}
    172 //          if(newdist!=kInfinity ){
    173 //              pre_pos-=pre_dir.unit()*newdist;
    174                 //pre_pos.z()=0.;
    175 //              PostStep->SetPosition(pre_pos);
    176 //              PostStep->SetMomentumDirection(pre_dir);
    177 //              post_pos=pre_pos;
    178 //              post_dir=pre_dir;
    179 //              double t = PostStep->GetLocalTime();
    180 //              t-=newdist/c_light;
    181 //              PostStep->SetLocalTime(t);
    182 //             
    183 //         
    184 //         
    185 //         
    186 //          }
    187 //         
    188 //     
    189         }
    190         else PostStep->SetPosition(pre_pos);
    191        
    192     }
    193 
    194  
    195     if ( PostStep && fStatus!=-4 && PostStep->GetStepStatus()==fPostStepDoItProc &&
    196          theStep->GetTrack()->GetTrackStatus()==fStopAndKill ) fStatus=-3;
    197 
    198     if(postvol){
     135  if( !fBoundary ){
     136    fProcessManager = theStep->GetTrack()->GetDefinition()->GetProcessManager();
     137    G4int nprocesses = fProcessManager->GetProcessListLength();
     138    fProcessVector = fProcessManager->GetProcessList();
     139
     140    for(int i=0;i<nprocesses;i++){
     141      if( (*fProcessVector)[i]->GetProcessName()=="OpBoundary" ){
     142        fBoundary = (G4OpBoundaryProcess*)(*fProcessVector)[i];
     143        break;
     144      }
     145    }
     146  }
     147
     148  theStatus = Undefined;
     149  if(fBoundary)theStatus = fBoundary->GetStatus();
     150
     151  double rho = post_pos.perp();
     152  if(!fTus && (fabs(post_pos.y())>=r_cut || rho>=r_wall)){
     153    if(cover->Inside (pre_pos)==kInside ){
     154      theStep->GetTrack()->SetTrackStatus(fStopAndKill);
     155      fStatus=-4;
     156      double dist = cover->DistanceToOut (pre_pos, pre_dir);
     157      // intersection with virtual wall
    199158     
    200         if ( PostStepVolume.contains("FresnelLens") || PreStepVolume.contains("FresnelLens")||
    201             PreStepVolume.contains("mirror_phys") ||PostStepVolume.contains("mirror_phys") ){
    202             if(fStatus==0) fStatus++;
    203         }
    204 
    205         if(PostStepVolume=="absorber"){
    206            
    207             theStep->GetTrack()->SetTrackStatus(fStopAndKill);
    208             fStatus=(fStatus>0?-2:-1);
    209             if(fStatus==-1||fStatus==-2)PostStep->SetPosition(pre_pos);
    210         }
    211     }
     159      if(dist!=kInfinity ){
     160                                pre_pos+=pre_dir.unit()*dist;
     161                                PostStep->SetPosition(pre_pos);
     162                                PostStep->SetMomentumDirection(pre_dir);
     163                                post_pos=pre_pos;
     164                                post_dir=pre_dir;
     165                                double t = PostStep->GetLocalTime();
     166                                t+=dist/c_light;
     167                                PostStep->SetLocalTime(t);
     168                        }
     169    }
     170    else if(cover->Inside (pre_pos)==kSurface)
     171      PostStep->SetPosition(pre_pos);
     172    else PostStep->SetPosition(pre_pos);
     173   
     174  }
     175 
     176  if ( PostStep && fStatus!=-4 && PostStep->GetStepStatus()==fPostStepDoItProc &&
     177    theStep->GetTrack()->GetTrackStatus()==fStopAndKill ) fStatus=-3;
     178
     179  if(postvol){
     180     
     181    if ( PostStepVolume.contains("FresnelLens") || PreStepVolume.contains("FresnelLens")||
     182      PreStepVolume.contains("mirror_phys") ||PostStepVolume.contains("mirror_phys") ){
     183      if(fStatus==0) fStatus++;
     184    }
     185
     186    if(PostStepVolume=="absorber"){
     187 
     188      theStep->GetTrack()->SetTrackStatus(fStopAndKill);
     189      fStatus=(fStatus>0?-2:-1);
     190      if(fStatus==-1||fStatus==-2)PostStep->SetPosition(pre_pos);
     191    }
     192   
     193    if(PostStepVolume=="iris"){
     194      fStatus=-5;
     195    }
     196  }
    212197
    213198   
    214     if ( theStatus==Absorption || theStatus==Detection ){
    215         fStatus=-3;
    216        
    217     }
    218     G4TrackStatus TrackStatus = theStep->GetTrack()->GetTrackStatus();
    219 
    220     if(TrackStatus==fStopAndKill){
    221         fStepCounter=0;
    222         int history;
    223         if ( post_pos.z()<0 && post_dir.z()<0  ){
    224             history=kEarth;
    225             double dist = cover->DistanceToOut (pre_pos, pre_dir);
    226         // intersection with virtual wall
    227            
    228             if(dist!=kInfinity && dist!=0){
    229 
    230                 pre_pos+=pre_dir.unit()*dist;
    231                 PostStep->SetPosition(pre_pos);
    232                 PostStep->SetMomentumDirection(pre_dir);
    233                 post_pos=pre_pos;
    234                 post_dir=pre_dir;
    235                 double t = PostStep->GetLocalTime();
    236                 t+=dist/c_light;
    237                 PostStep->SetLocalTime(t);
    238             }
    239          
    240          }
    241          else if ( fStatus==-2 ) history=kG4FocalSurface;
    242          else if ( fStatus==-3) history=kOptics;
    243          else if ( fStatus==-4) history=kWalls;
    244          else history=keUndefined;
    245 
    246          fStacking->SetTrackHistory(history);
    247          fStatus=0;
    248     }
     199  if ( theStatus==Absorption || theStatus==Detection ){
     200    fStatus=-3;
     201
     202  }
     203  G4TrackStatus TrackStatus = theStep->GetTrack()->GetTrackStatus();
     204
     205  if(TrackStatus==fStopAndKill){
     206    fStepCounter=0;
     207    int history;
     208 
     209    if ( post_pos.z()<=0 && post_dir.z()<0 && !fTus){
     210      history=kEarth;
     211      double dist = cover->DistanceToOut (pre_pos, pre_dir);
     212  // intersection with virtual wall
     213      if(dist!=kInfinity && dist!=0){
     214        pre_pos+=pre_dir.unit()*dist;
     215        PostStep->SetPosition(pre_pos);
     216        PostStep->SetMomentumDirection(pre_dir);
     217        post_pos=pre_pos;
     218        post_dir=pre_dir;
     219        double t = PostStep->GetLocalTime();
     220        t+=dist/c_light;
     221        PostStep->SetLocalTime(t);
     222      }
     223
     224    }
     225    else if ( fStatus==-2 ) history=kG4FocalSurface;
     226    else if ( fStatus==-3) history=kOptics;
     227    else if ( fStatus==-4) history=kWalls;
     228    else if ( fStatus==-5) history=kIris;
     229    else history=keUndefined;
     230
     231    fStacking->SetTrackHistory(history);
     232    fStatus=0;
     233  }
    249234}

  • esaf_lal/branches/camille/packages/simulation/detector/G4Detector/optics/src/OptUserStackingAction.cc

    r16 r184  
    3737    kWalls                    = 5,
    3838    kEarth                    = 6,
    39     kG4FocalSurface           = 7
     39    kIris                     = 7,
     40    kG4FocalSurface           = 8
    4041};
    4142
    … …  
    7172    if (det.find("Tus")!=string::npos)ftus=true;
    7273    else ftus=false;
    73 
     74               
    7475    SetSender("OptStackingAction");
    7576}
    … …  
    7980    // fFiller is deleted by Geant:
    8081   // delete info;
     82   MsgForm(EsafMsg::Info,"fNTracked = %i.",fNTracked );
     83         MsgForm(EsafMsg::Info,"walls = %i.",fNwalls );
     84         MsgForm(EsafMsg::Info,"optics = %i.",fNoptics );
     85         MsgForm(EsafMsg::Info,"focal plane = %i.",fNfocalplane );
     86         MsgForm(EsafMsg::Info,"earth = %i.",fNearth );
     87   MsgForm(EsafMsg::Info,"iris = %i.",fNiris );
     88         MsgForm(EsafMsg::Info,"lost = %i.",fNlost );
    8189}
    8290
    … …  
    120128        fPhotonDef = particleTable->FindParticle("opticalphoton");
    121129    }
    122 
    123130    EVector& phdir=ph->dir;
    124131    if(ftus)phdir*=-1.;
    … …  
    127134                                                    G4ThreeVector(phdir[X],phdir[Y],phdir[Z]),
    128135                                                    energy);
    129 
    130136    EVector& phpos=ph->pos;
    131137    /// system of ESAF is copied from CLHEP, so no change is needed
    … …  
    146152
    147153    if(refl)return ph_tr;
    148     if ( fEv ) {
     154    if (fEv) {
    149155        EDetectorPhotonAdder a(ph,0,kTRUE);
    150156        fEv->Fill(a);
    … …  
    177183        double theConvConst=hbarc*twopi/eV;
    178184        ph->wl=theConvConst/energy*nm;
    179         ParentPhoton *p = new ParentPhoton(fNG4Tracked,x, y, dir.theta(),           dir.phi(),theConvConst/energy,time,(ParentPhotonType)2);
     185        ParentPhoton *p = new ParentPhoton(fNG4Tracked,x, y, dir.theta(),
     186                                           dir.phi(),theConvConst/energy,time,(ParentPhotonType)2);
    180187        ph->parent=p;
    181188        track->SetUserInformation(new G4TrackInformation(ph,0,1));
    … …  
    247254//______________________________________________________________________________
    248255void OptUserStackingAction::SetPupil(PhotonsOnPupil* p){
    249 
     256    printf("**********Set event\n");
     257    Reset();
     258    fEv = EEvent::GetCurrent();
    250259    if(!p) return;
    251260
    252     Reset();
     261   
    253262    fPupil=p;
    254263
    255264    fNTotal=fPupil->GetNphotons();
    256     fEv = EEvent::GetCurrent();
    257     fTimeFirstPhoton = fHuge;
    258     fTimeLastPhoton = -fHuge;
     265   
    259266}
    260267
    … …  
    262269void OptUserStackingAction::Reset(){
    263270
    264     fNG4focalplane=0;
    265     fNfocalplane=0;
    266     fNlost=0;
    267     fNTracked=0;
    268     fNG4Tracked=0;
    269     fProgress=0;
    270     forigin=0;
     271    fNG4focalplane = 0;
     272    fNfocalplane = 0;
     273    fNlost = 0;
     274                fNearth = 0;
     275                fNoptics = 0;
     276                fNwalls = 0;
     277    fNTracked = 0;
     278    fNiris = 0;
     279    fNG4Tracked = 0;
     280    fProgress = 0;
     281    forigin = 0;
     282   
     283    fNTotal = 0;
     284   
     285    fEv = 0;
     286   
    271287    fTimeFirstPhoton = fHuge;
    272288    fTimeLastPhoton = -fHuge;
    … …  
    290306        Photon* ph = (Photon*)this->UpdateCurrentPhoton((G4Track*)fTrackPost);
    291307
    292 
    293         if ( this->GoToElectronics(ph) ){
    294 
    295             return;
     308        if ( this->GoToElectronics(ph) ){         
     309          return;
    296310        }
    297 
     311       
    298312        switch ( ph->history ) {
    299             case kFocalPlane:
    300                 fNfocalplane++;
    301                 break;
    302             case kG4FocalSurface :
    303                 fNG4focalplane++;
    304                 break;
    305             case kOptics:
    306                 fNoptics++;
    307                 break;
    308             default:
    309                 fNlost++;
     313          case kFocalPlane:
     314            fNfocalplane++;
     315            break;
     316          case kG4FocalSurface :
     317            fNG4focalplane++;
     318            break;
     319          case kOptics:
     320            fNoptics++;
     321            break;
     322          case kWalls:
     323            fNwalls++;
     324            break;
     325          case kEarth:
     326            fNearth++;
     327            break;
     328          case kIris:
     329            fNiris++;
     330            break;
     331          default:
     332            fNlost++;
     333            printf("%i history *** x,y,z (%g %g %g)\n",ph->history, ph->pos[X],ph->pos[Y],ph->pos[Z]);
    310334        }
    311 
     335       
    312336        if( fEv  ) {
    313             EDetectorPhotonAdder a(ph,0,kFALSE);
    314             fEv->Fill(a);
    315         }
    316     }
     337//           bool is_new = kFALSE;
     338//           if(!fNTotal)is_new = kTRUE;
     339//           EDetectorPhotonAdder a(ph,0,is_new);
     340//           fEv->Fill(a);
     341        }
     342       
     343    }
     344   
     345   
    317346    if (fNTotal){
    318         fNTracked++;
    319         this->Progress(100.*fNTracked/fNTotal,"Photons");
     347      fNTracked++;
     348      this->Progress(100.*fNTracked/fNTotal,"Photons");
    320349    }
    321350}
    … …  
    342371    }
    343372    EVector intPoint;
    344     ph->history=kFocalPlane;
     373//     ph->history=kFocalPlane;
    345374
    346375    if(ftus){
    … …  
    350379    else fIdealFocalSurface->HitPosition(ph);
    351380
    352 
     381    if(ph->hitIfs)ph->history=kFocalPlane;
     382     
    353383    double DZ= fFocalPlane->Bottom()-ph->pos[Z];
    354384    // calculate impact point on the focal plane
    … …  
    356386    ph->time+=(ph->pos - intPoint).Mag()/Clight();
    357387    ph->pos=intPoint;
    358     Photon* newp=fFocalPlane->Transport(ph);
     388//     printf("1new!\n");
     389   
     390    Photon* newp = fFocalPlane->Transport(ph);
     391//     printf("0new!\n");
    359392    if ( !newp ) return false;
    360 
     393//     printf("new!\n");
    361394    stackManager->PushOneTrack( this->NewPhoton(newp,newp->id,1) );
    362395    return true;
    … …  
    370403    if (  (curn+1)>=n ) {
    371404        if ( fEv ) {
    372                 TArrayI photonHistory(8);
     405                TArrayI photonHistory(9);
    373406                photonHistory.AddAt(fNlost, 0);
    374407                photonHistory.AddAt(fNG4focalplane, kG4FocalSurface);

  • esaf_lal/branches/camille/packages/simulation/detector/GNUmakefile

    r16 r184  
    11# Dummy MakeFile
    2 # $Id: GNUmakefile 2856 2010-12-23 19:37:08Z biktem $
     2# $Id: GNUmakefile 3032 2013-07-03 11:52:38Z biktem $
    33# D. De Marco created Jan, 23 2002
    44#
    … …  
    99
    1010include $(ESAFINSTALL)/packages/config.gmk
    11 #include $(G4INSTALL)/config/binmake.gmk
    1211SUBDIRS = electronics optics tools
    1312
    14 ifdef G4INSTALL
     13ifdef USE_G4
    1514SUBDIRS += G4Detector
    1615endif
Note: See TracChangeset for help on using the changeset viewer.