source: JEM-EUSO/esaf_cc_at_lal/packages/simulation/detector/G4Detector/optics/src/OptUserStackingAction.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: 11.3 KB
Line 
1#include "OptUserStackingAction.hh"
2
3#include "OpticsFactory.hh"
4#include "IdealFocalSurface.hh"
5#include "EDetectorPhotonAdder.hh"
6#include "EVector.hh"
7#include "EEvent.hh"
8#include "EConst.hh"
9#include "PhotonsOnPupil.hh"
10#include "utils.hh"
11#include "ConfigFileParser.hh"
12#include "EsafConfigurable.hh"
13
14#include <G4ParticleTable.hh>
15#include <G4DynamicParticle.hh>
16#include <G4Track.hh>
17#include <Randomize.hh>
18#include <G4StackManager.hh>
19#include <G4RunManager.hh>
20#include <G4Run.hh>
21#include <G4TrackInformation.hh>
22#include <G4VUserTrackInformation.hh>
23#include <TArrayI.h>
24#include <globals.hh>
25
26#include "ParentPhoton.hh"
27
28using namespace EConst;
29// Detector subsystem identifier
30// redefine to avoid kUndefined being defined twice
31enum EDetectorSystem {
32    keUndefined               = 0,
33    kNotInside                = 1,
34    kBaffle                   = 2,
35    kOptics                   = 3,
36    kFocalPlane               = 4,
37    kWalls                    = 5,
38    kEarth                    = 6,
39    kG4FocalSurface           = 7
40};
41
42void SetOptPhotonPolar(G4Track* track)
43{
44    G4double angle = G4UniformRand() * 360.0*deg;
45
46    G4ThreeVector normal (1., 0., 0.);
47    G4ThreeVector kphoton = track->GetMomentumDirection();
48    G4ThreeVector product = normal.cross(kphoton);
49    G4double modul2       = product*product;
50
51    G4ThreeVector e_perpend (0., 0., 1.);
52    if (modul2 > 0.) e_perpend = (1./std::sqrt(modul2))*product;
53    G4ThreeVector e_paralle    = e_perpend.cross(kphoton);
54
55    G4ThreeVector polar = std::cos(angle)*e_paralle + std::sin(angle)*e_perpend;
56    track->SetPolarization(polar);
57}
58
59//_______________________________________________________________________________
60OptUserStackingAction::OptUserStackingAction(IdealFocalSurface* fs) :
61fPhotonDef(0),
62fIdealFocalSurface(fs),
63fHuge(1e20),
64fNmax(-1),
65fTrackHistory(0)
66{
67    fFocalPlane=OpticsFactory::Get()->GetFocalPlane();
68    this->Reset();
69    ConfigFileParser* conf = Config::Get()->GetCF("G4","G4Detector");
70    string det=conf->GetStr("G4Detector.DetectorType");
71    if (det.find("Tus")!=string::npos)ftus=true;
72    else ftus=false;
73
74    SetSender("OptStackingAction");
75}
76
77//_______________________________________________________________________________
78OptUserStackingAction::~OptUserStackingAction(){
79    // fFiller is deleted by Geant:
80   // delete info;
81}
82
83
84//_______________________________________________________________________________
85void OptUserStackingAction::CheckTrace(int id,
86                                           const G4ThreeVector& pos,
87                                           const G4ThreeVector& mom,
88                                           double lambda,
89                                           double time)
90{
91    NTraceLens::NPhoton ph_in;
92    ph_in.id=id;
93    ph_in.x[0]=pos[0]/mm;
94    ph_in.x[1]=pos[1]/mm;
95    ph_in.x[2]=pos[2]/mm;
96    ph_in.px[0]=mom[0];
97    ph_in.px[1]=mom[1];
98    ph_in.px[2]=mom[2];
99    ph_in.id=id;
100    ph_in.lambda=lambda;
101    ph_in.t=time;
102    ph_in.history=0;
103
104
105#ifdef TRACEDEBUG
106    cout<<"NOpticalSystem"<<endl;
107    cout<<"position = "<<p.pos<<endl;
108    cout<<"direction = "<<p.dir<<endl;
109    double theta=p.dir.Angle(EVector(0.,0.,1.));
110    cout<<"theta = "<<theta/deg<<endl;
111    cout<<"wl = "<<p.wl/nm<<endl;
112    cout<<"time = "<<p.time/ns<<endl;
113#endif /* DEBUG */
114}
115
116//_______________________________________________________________________________
117G4Track* OptUserStackingAction::NewPhoton(Photon* ph, int id, int refl){
118    if ( !fPhotonDef ) {
119        G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();
120        fPhotonDef = particleTable->FindParticle("opticalphoton");
121    }
122
123    EVector& phdir=ph->dir;
124    if(ftus)phdir*=-1.;
125    double energy=hbarc*twopi/ph->wl; // lambda in mm to energy in MeV
126    G4DynamicParticle* ph_dyn=new G4DynamicParticle(fPhotonDef,
127                                                    G4ThreeVector(phdir[X],phdir[Y],phdir[Z]),
128                                                    energy);
129
130    EVector& phpos=ph->pos;
131    /// system of ESAF is copied from CLHEP, so no change is needed
132    G4ThreeVector pos(phpos[X],phpos[Y],phpos[Z]);
133    G4ThreeVector mom(phdir[X],phdir[Y],phdir[Z]);
134    double time;
135    if(!ftus && !refl){
136        double shift=100.*mm;
137        pos-=mom.unit()*shift;
138        time=ph->time-shift/c_light;
139    }
140    else time=ph->time;
141    ph->time=time;
142    G4Track* ph_tr=new G4Track( ph_dyn, time, pos);
143    ph_tr->SetUserInformation(new G4TrackInformation(ph,refl,forigin));
144    ph_tr->SetTrackID(id);
145    SetOptPhotonPolar(ph_tr);
146
147    if(refl)return ph_tr;
148    if ( fEv ) {
149        EDetectorPhotonAdder a(ph,0,kTRUE);
150        fEv->Fill(a);
151    }
152    if ( ph->time > fTimeLastPhoton ) fTimeLastPhoton = ph->time;
153    if ( ph->time < fTimeFirstPhoton ) fTimeFirstPhoton = ph->time;
154    return ph_tr;
155}
156
157//_______________________________________________________________________________
158Photon* OptUserStackingAction::UpdateCurrentPhoton(G4Track* track){
159    G4Step* st = (G4Step*)track->GetStep();
160    G4StepPoint *PostStep = st->GetPostStepPoint();
161    G4ThreeVector pos =  PostStep->GetPosition();
162    G4ThreeVector dir =  PostStep->GetMomentumDirection();
163    double time = PostStep->GetLocalTime();
164    G4TrackInformation* iph = (G4TrackInformation*)track->GetUserInformation();
165    Photon* ph;
166    double x[3] = {0.,0.,0.};
167    double y[3] = {pos.x(),pos.y(),pos.z()};
168
169    if(!iph){
170        ph = new Photon();
171        forigin=1;
172
173        ph->time=0;
174        ph->parent=0;
175        ph->id=fNG4Tracked;
176        double energy=track->GetKineticEnergy();
177        double theConvConst=hbarc*twopi/eV;
178        ph->wl=theConvConst/energy*nm;
179        ParentPhoton *p = new ParentPhoton(fNG4Tracked,x, y, dir.theta(),           dir.phi(),theConvConst/energy,time,(ParentPhotonType)2);
180        ph->parent=p;
181        track->SetUserInformation(new G4TrackInformation(ph,0,1));
182        fNG4Tracked++;
183
184    }
185    else ph = (Photon*)iph->fPhoton;
186    ph->pos.SetXYZ(pos.x(),pos.y(),pos.z());
187    ph->dir.SetXYZ(dir.x(),dir.y(),dir.z());
188    ph->time+=time;
189    ph->history=fTrackHistory;
190    return ph;
191}
192
193//_______________________________________________________________________________
194G4ClassificationOfNewTrack OptUserStackingAction::ClassifyNewTrack(const G4Track* aTrack){
195
196//     if ( !fFiller ) {
197//
198//         return fKill;
199//     }
200//     if ( aTrack==fFiller ) {
201//
202//         return fWaiting;
203//     }
204
205    return fUrgent;
206}
207
208//______________________________________________________________________________
209void OptUserStackingAction::NewStage(){
210
211    // get a new photon
212    if(!fNTotal)return;
213
214    Photon* photon = fPupil->GetPhoton();
215
216    // build an empty EDetPhoton in the event
217    if ( !photon || ( fNmax>0 && fNTracked>=fNmax ) ) {
218
219
220        if ( fEv ) {
221            TArrayI photonHistory(8);
222            photonHistory.AddAt(fNlost, 0);
223            photonHistory.AddAt(fNG4focalplane, kG4FocalSurface);
224            photonHistory.AddAt(fNfocalplane, kFocalPlane);
225            photonHistory.AddAt(fNoptics, kOptics);
226
227            EDetectorPhotonAdder a(photonHistory);
228            fEv->Fill(a);
229        }
230
231        MsgForm(EsafMsg::Info,"Number of events processed = %i.",fNTracked );
232        return;
233    }
234
235    forigin=0;
236    stackManager->PushOneTrack( NewPhoton(photon, fNTracked,0) );
237}
238
239
240//_______________________________________________________________________________
241void OptUserStackingAction::PrepareNewEvent(){
242
243    this->NewStage();
244}
245
246
247//______________________________________________________________________________
248void OptUserStackingAction::SetPupil(PhotonsOnPupil* p){
249
250    if(!p) return;
251
252    Reset();
253    fPupil=p;
254
255    fNTotal=fPupil->GetNphotons();
256    fEv = EEvent::GetCurrent();
257    fTimeFirstPhoton = fHuge;
258    fTimeLastPhoton = -fHuge;
259}
260
261//__________________________________________________________________
262void OptUserStackingAction::Reset(){
263
264    fNG4focalplane=0;
265    fNfocalplane=0;
266    fNlost=0;
267    fNTracked=0;
268    fNG4Tracked=0;
269    fProgress=0;
270    forigin=0;
271    fTimeFirstPhoton = fHuge;
272    fTimeLastPhoton = -fHuge;
273}
274
275//__________________________________________________________________
276EVector OptUserStackingAction::NewInteractionPoint(const Photon *p, double DZ) const {
277    EVector dir=p->dir.Unit();
278    EVector newPoint=p->pos+dir*(DZ/dir[Z]);
279#ifdef TRACEDEBUG
280    cout<<"new intPoint: "<<newPoint<<endl;
281#endif /* DEBUG */
282    return newPoint;
283}
284
285
286//-----------------------------------------------------------------
287void OptUserStackingAction::PostUserTrackingAction(const G4Track* fTrackPost){
288
289    if( fTrackPost->GetDefinition() == G4OpticalPhoton::OpticalPhotonDefinition()){
290        Photon* ph = (Photon*)this->UpdateCurrentPhoton((G4Track*)fTrackPost);
291
292
293        if ( this->GoToElectronics(ph) ){
294
295            return;
296        }
297
298        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++;
310        }
311
312        if( fEv  ) {
313            EDetectorPhotonAdder a(ph,0,kFALSE);
314            fEv->Fill(a);
315        }
316    }
317    if (fNTotal){
318        fNTracked++;
319        this->Progress(100.*fNTracked/fNTotal,"Photons");
320    }
321}
322//______________________________________________________________________________
323void OptUserStackingAction::Progress(int n,const char* p){
324    int prog =n;
325
326    if ( prog%10==0 && prog>=fProgress ) {
327        Msg(EsafMsg::Info).SetProgress(fProgress);             //fProgress
328        Msg(EsafMsg::Info) << Form("Tracking %s:  ",p) << MsgCount;
329        fProgress+=10;
330
331    }
332
333
334}
335//______________________________________________________________________________
336bool OptUserStackingAction::GoToElectronics(Photon* ph){
337
338    if ( ph->history!= kG4FocalSurface){
339//         Msg(EsafMsg::Debug) << Form("photon isn't on FocalPlane %d",ph->history) << MsgDispatch;
340//         Msg(EsafMsg::Debug) << "@ new photon " << MsgDispatch;
341        return false;
342    }
343    EVector intPoint;
344    ph->history=kFocalPlane;
345
346    if(ftus){
347        ph->posOnIfs= ph->pos;
348        ph->hitIfs=kTRUE;
349    }
350    else fIdealFocalSurface->HitPosition(ph);
351
352
353    double DZ= fFocalPlane->Bottom()-ph->pos[Z];
354    // calculate impact point on the focal plane
355    intPoint=NewInteractionPoint(ph, DZ);
356    ph->time+=(ph->pos - intPoint).Mag()/Clight();
357    ph->pos=intPoint;
358    Photon* newp=fFocalPlane->Transport(ph);
359    if ( !newp ) return false;
360
361    stackManager->PushOneTrack( this->NewPhoton(newp,newp->id,1) );
362    return true;
363}
364
365void OptUserStackingAction::EndOfEventAction(const G4Event* anEvent){
366    int n=G4RunManager::GetRunManager()->GetCurrentRun()->GetNumberOfEventToBeProcessed();
367    int curn=G4RunManager::GetRunManager()->GetCurrentRun()->GetNumberOfEvent();
368    this->Progress(100.*curn/n,"events");
369
370    if (  (curn+1)>=n ) {
371        if ( fEv ) {
372                TArrayI photonHistory(8);
373                photonHistory.AddAt(fNlost, 0);
374                photonHistory.AddAt(fNG4focalplane, kG4FocalSurface);
375                photonHistory.AddAt(fNfocalplane, kFocalPlane);
376                photonHistory.AddAt(fNoptics, kOptics);
377                EDetectorPhotonAdder a(photonHistory);
378                fEv->Fill(a);
379        }
380    }
381
382}
383
384void OptUserStackingAction::PreUserTrackingAction(const G4Track* fTrackPre){
385
386}
387void OptUserStackingAction::BeginOfEventAction(const G4Event* anEvent){
388}
Note: See TracBrowser for help on using the repository browser.