source: JEM-EUSO/esaf_cc_at_lal/packages/simulation/detector/G4Detector/optics/src/OpEUSODiffraction.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: 7.1 KB
Line 
1#include "OpEUSODiffraction.hh"
2
3#include <G4GeometryTolerance.hh>
4#include <G4Step.hh>
5#include <G4DynamicParticle.hh>
6#include <G4Material.hh>
7#include <G4OpticalPhoton.hh>
8#include <G4TransportationManager.hh>
9#include <G4Navigator.hh>
10
11#include <TGraph.h>
12#include "Ntrace_optics.hh"
13using namespace NTraceLens;
14
15//______________________________________________________________________________
16OpEUSODiffraction::OpEUSODiffraction(const G4String& processName, G4ProcessType type) :
17G4VDiscreteProcess(processName, type) {
18        this->SetSender("OpEUSODiffraction");
19        Msg(EsafMsg::Debug) << "OpEUSODiffraction is created." << MsgDispatch;
20
21        kCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
22        const char* lensdir = Conf()->GetStr("G4EusoGeometryShape.lens_dir").data();
23        char* name=fOpticalSystemParam->name;
24        const char* dindex_fname = Form("%sDOE_%s.dat",lensdir,name);
25        theDindex=new TGraph(dindex_fname);
26        if ( !theDindex || theDindex->IsZombie() ){
27            Msg(EsafMsg::Panic) << "File with DOE not found(" << dindex_fname << ")" << MsgDispatch;
28        }
29        double y;
30        theDindex->GetPoint(theDindex->GetN()-1,theMaximumPerp,y);
31
32        theLambda0=357.0; ///nm
33
34}
35
36//______________________________________________________________________________
37OpEUSODiffraction::~OpEUSODiffraction(){
38    delete theDindex;
39}
40
41//______________________________________________________________________________
42G4bool OpEUSODiffraction::IsApplicable(const G4ParticleDefinition& aParticleType){
43   return ( &aParticleType == G4OpticalPhoton::OpticalPhoton() );
44}
45
46//______________________________________________________________________________
47G4VParticleChange*
48OpEUSODiffraction::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep)
49{
50    ///
51    /// partly taken from the G4OpBoundaryProcess::PostStepDoIt() implementation
52    ///
53    aParticleChange.Initialize(aTrack);
54
55    if (aTrack.GetStepLength()<=kCarTolerance*0.5){
56            theStatus = edStepTooSmall;
57            return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
58    }
59
60    G4StepPoint* pPreStepPoint  = aStep.GetPreStepPoint();
61    G4StepPoint* pPostStepPoint = aStep.GetPostStepPoint();
62
63    if ( pPostStepPoint->GetStepStatus()!=fGeomBoundary ){
64        theStatus = edNotAtBoundary;
65        return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
66    }
67
68    ///dynamic particle parameters
69    const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
70    thePhotonMomentum = aParticle->GetTotalMomentum();
71    OldMomentum       = aParticle->GetMomentumDirection();
72    OldPolarization   = aParticle->GetPolarization();
73
74    ///previous and next material parameters
75    Material1 = pPreStepPoint->GetMaterial();
76    Material2 = pPostStepPoint->GetMaterial();
77    if ( Material1==Material2 ){
78        theStatus = edSameMaterial;
79        return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
80    }
81
82    G4MaterialPropertyVector* Rindex;
83    G4MaterialPropertiesTable* aMaterialPropertiesTable = Material1->GetMaterialPropertiesTable();
84    if ( aMaterialPropertiesTable ) {
85        Rindex = aMaterialPropertiesTable->GetProperty("EUSO_RINDEX");
86    } else {
87        theStatus = edNoINDEX;
88        return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
89    }
90    if ( !aMaterialPropertiesTable->GetProperty("EUSO_DINDEX") ){
91        theStatus = edNoINDEX;
92        return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
93    }
94    if ( Rindex ) {
95        Rindex1 = Rindex->GetProperty(thePhotonMomentum);
96    } else {
97        theStatus = edNoINDEX;
98        return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
99    }
100
101    aMaterialPropertiesTable = Material2->GetMaterialPropertiesTable();
102    if (aMaterialPropertiesTable) {
103        Rindex = aMaterialPropertiesTable->GetProperty("RINDEX");
104    } else {
105        theStatus = edNoINDEX;
106        return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
107    }
108
109    if ( Rindex ) {
110        Rindex2 = Rindex->GetProperty(thePhotonMomentum);
111    } else {
112        theStatus = edNoINDEX;
113        return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
114    }
115
116    G4ThreeVector theGlobalPoint = pPostStepPoint->GetPosition();
117    double perp=theGlobalPoint.perp();
118    if ( perp > theMaximumPerp ){
119        return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
120    }
121    Dindex1 = theDindex->Eval( perp );
122
123///adopted from Ltrace_optF1v4
124/*
125 * ------------------------------------------------------------------------
126 *  calculate the diffracted ray
127 *  Input:
128 *        theGlobalPoint   array x,y,z of photon
129 *        nnn[3]        normal vector to the diffractive surface at the intersection
130 *        Rindex1            refractive index of xxx
131 *        Rindex2            refractive index of xxx
132 *        Dindex1
133 *        wl            wavelength in mm
134 *        m             index
135 *  Output:
136 *        in_ph[3]
137 * ------------------------------------------------------------------------
138 */
139     G4ThreeVector nnnR(theGlobalPoint.x()/perp,theGlobalPoint.y()/perp,0.0);
140//     G4ThreeVector nnn(.0,.0,1.0);
141
142    G4Navigator* theNavigator = G4TransportationManager::GetTransportationManager()->GetNavigatorForTracking();
143    G4ThreeVector nnn;   // Normal points back into volume
144    G4bool valid;
145    nnn = theNavigator->GetLocalExitNormal(&valid);
146    if ( !valid ) {
147        printf("Error: invalid normal\n");
148        return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
149    }
150
151    Rindex1=Rindex2=1.0003;                             // as in Ltrace, I don't know why
152    double dir=OldMomentum.z()>0.?1.0:-1.0;
153    double CosIN=OldMomentum*nnn;
154    double CosD=nnn*nnnR;
155
156    double m=1.;                                        // as in Ltrace
157    double wl=hbarc*twopi/aTrack.GetTotalEnergy()*1.e6; // energy in MeV to lambda in nm
158    double c1=Rindex1/Rindex2;
159    double c2=wl/theLambda0*Dindex1*dir*m/Rindex2;
160    G4ThreeVector ppp=(OldMomentum-nnn*CosIN)*c1+(nnnR-nnn*CosD)*c2;
161    double norm=ppp.mag2();
162
163    //#define DEBUG 1
164    if(norm>1.0) {
165        #ifdef DEBUG
166              printf("norm>1, return\n");
167        #endif
168        return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
169    }
170    norm=sqrt(1.0-norm);
171
172    NewMomentum=ppp+nnn*norm*dir;
173#ifdef DEBUG
174    cout<<"our"<<endl;
175    printf("(%.3f , %.3f , %.3f)\n", theGlobalPoint.x(),theGlobalPoint.y(),theGlobalPoint.z());
176    printf("(%.3f , %.3f , %.3f)\n", nnnR.x(),nnnR.y(),nnnR.z());
177    printf("(%.3f , %.3f , %.3f)\n", Rindex1,Dindex1,CosIN);
178    printf("ppp: (%.3f , %.3f , %.3f)\n", ppp.x(),ppp.y(),ppp.z());
179    printf("Diffractive: (%.3f , %.3f , %.3f)=>", OldMomentum.x(),OldMomentum.y(),OldMomentum.z());
180    printf("(%.3f , %.3f , %.3f)\n", NewMomentum.x(),NewMomentum.y(),NewMomentum.z());
181    printf("%g -> %g \n", OldMomentum.mag(), NewMomentum.mag());
182#endif
183
184    aParticleChange.ProposeMomentumDirection(NewMomentum);
185
186    return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
187
188}
189
190
191//______________________________________________________________________________
192G4double OpEUSODiffraction::GetMeanFreePath(const G4Track& ,
193                                              G4double ,
194                                              G4ForceCondition* condition)
195{
196        *condition = Forced;
197        return DBL_MAX;
198}
Note: See TracBrowser for help on using the repository browser.