source: JEM-EUSO/esaf_lal/branches/camille/packages/simulation/detector/G4Detector/optics/src/OptUserStackingAction.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: 12.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    kIris                     = 7,
40    kG4FocalSurface           = 8
41};
42
43void SetOptPhotonPolar(G4Track* track)
44{
45    G4double angle = G4UniformRand() * 360.0*deg;
46
47    G4ThreeVector normal (1., 0., 0.);
48    G4ThreeVector kphoton = track->GetMomentumDirection();
49    G4ThreeVector product = normal.cross(kphoton);
50    G4double modul2       = product*product;
51
52    G4ThreeVector e_perpend (0., 0., 1.);
53    if (modul2 > 0.) e_perpend = (1./std::sqrt(modul2))*product;
54    G4ThreeVector e_paralle    = e_perpend.cross(kphoton);
55
56    G4ThreeVector polar = std::cos(angle)*e_paralle + std::sin(angle)*e_perpend;
57    track->SetPolarization(polar);
58}
59
60//_______________________________________________________________________________
61OptUserStackingAction::OptUserStackingAction(IdealFocalSurface* fs) :
62fPhotonDef(0),
63fIdealFocalSurface(fs),
64fHuge(1e20),
65fNmax(-1),
66fTrackHistory(0)
67{
68    fFocalPlane=OpticsFactory::Get()->GetFocalPlane();
69    this->Reset();
70    ConfigFileParser* conf = Config::Get()->GetCF("G4","G4Detector");
71    string det=conf->GetStr("G4Detector.DetectorType");
72    if (det.find("Tus")!=string::npos)ftus=true;
73    else ftus=false;
74               
75    SetSender("OptStackingAction");
76}
77
78//_______________________________________________________________________________
79OptUserStackingAction::~OptUserStackingAction(){
80    // fFiller is deleted by Geant:
81   // 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 );
89}
90
91
92//_______________________________________________________________________________
93void OptUserStackingAction::CheckTrace(int id,
94                                           const G4ThreeVector& pos,
95                                           const G4ThreeVector& mom,
96                                           double lambda,
97                                           double time)
98{
99    NTraceLens::NPhoton ph_in;
100    ph_in.id=id;
101    ph_in.x[0]=pos[0]/mm;
102    ph_in.x[1]=pos[1]/mm;
103    ph_in.x[2]=pos[2]/mm;
104    ph_in.px[0]=mom[0];
105    ph_in.px[1]=mom[1];
106    ph_in.px[2]=mom[2];
107    ph_in.id=id;
108    ph_in.lambda=lambda;
109    ph_in.t=time;
110    ph_in.history=0;
111
112
113#ifdef TRACEDEBUG
114    cout<<"NOpticalSystem"<<endl;
115    cout<<"position = "<<p.pos<<endl;
116    cout<<"direction = "<<p.dir<<endl;
117    double theta=p.dir.Angle(EVector(0.,0.,1.));
118    cout<<"theta = "<<theta/deg<<endl;
119    cout<<"wl = "<<p.wl/nm<<endl;
120    cout<<"time = "<<p.time/ns<<endl;
121#endif /* DEBUG */
122}
123
124//_______________________________________________________________________________
125G4Track* OptUserStackingAction::NewPhoton(Photon* ph, int id, int refl){
126    if ( !fPhotonDef ) {
127        G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();
128        fPhotonDef = particleTable->FindParticle("opticalphoton");
129    }
130    EVector& phdir=ph->dir;
131    if(ftus)phdir*=-1.;
132    double energy=hbarc*twopi/ph->wl; // lambda in mm to energy in MeV
133    G4DynamicParticle* ph_dyn=new G4DynamicParticle(fPhotonDef,
134                                                    G4ThreeVector(phdir[X],phdir[Y],phdir[Z]),
135                                                    energy);
136    EVector& phpos=ph->pos;
137    /// system of ESAF is copied from CLHEP, so no change is needed
138    G4ThreeVector pos(phpos[X],phpos[Y],phpos[Z]);
139    G4ThreeVector mom(phdir[X],phdir[Y],phdir[Z]);
140    double time;
141    if(!ftus && !refl){
142        double shift=100.*mm;
143        pos-=mom.unit()*shift;
144        time=ph->time-shift/c_light;
145    }
146    else time=ph->time;
147    ph->time=time;
148    G4Track* ph_tr=new G4Track( ph_dyn, time, pos);
149    ph_tr->SetUserInformation(new G4TrackInformation(ph,refl,forigin));
150    ph_tr->SetTrackID(id);
151    SetOptPhotonPolar(ph_tr);
152
153    if(refl)return ph_tr;
154    if (fEv) {
155        EDetectorPhotonAdder a(ph,0,kTRUE);
156        fEv->Fill(a);
157    }
158    if ( ph->time > fTimeLastPhoton ) fTimeLastPhoton = ph->time;
159    if ( ph->time < fTimeFirstPhoton ) fTimeFirstPhoton = ph->time;
160    return ph_tr;
161}
162
163//_______________________________________________________________________________
164Photon* OptUserStackingAction::UpdateCurrentPhoton(G4Track* track){
165    G4Step* st = (G4Step*)track->GetStep();
166    G4StepPoint *PostStep = st->GetPostStepPoint();
167    G4ThreeVector pos =  PostStep->GetPosition();
168    G4ThreeVector dir =  PostStep->GetMomentumDirection();
169    double time = PostStep->GetLocalTime();
170    G4TrackInformation* iph = (G4TrackInformation*)track->GetUserInformation();
171    Photon* ph;
172    double x[3] = {0.,0.,0.};
173    double y[3] = {pos.x(),pos.y(),pos.z()};
174
175    if(!iph){
176        ph = new Photon();
177        forigin=1;
178
179        ph->time=0;
180        ph->parent=0;
181        ph->id=fNG4Tracked;
182        double energy=track->GetKineticEnergy();
183        double theConvConst=hbarc*twopi/eV;
184        ph->wl=theConvConst/energy*nm;
185        ParentPhoton *p = new ParentPhoton(fNG4Tracked,x, y, dir.theta(), 
186                                           dir.phi(),theConvConst/energy,time,(ParentPhotonType)2);
187        ph->parent=p;
188        track->SetUserInformation(new G4TrackInformation(ph,0,1));
189        fNG4Tracked++;
190
191    }
192    else ph = (Photon*)iph->fPhoton;
193    ph->pos.SetXYZ(pos.x(),pos.y(),pos.z());
194    ph->dir.SetXYZ(dir.x(),dir.y(),dir.z());
195    ph->time+=time;
196    ph->history=fTrackHistory;
197    return ph;
198}
199
200//_______________________________________________________________________________
201G4ClassificationOfNewTrack OptUserStackingAction::ClassifyNewTrack(const G4Track* aTrack){
202
203//     if ( !fFiller ) {
204//
205//         return fKill;
206//     }
207//     if ( aTrack==fFiller ) {
208//
209//         return fWaiting;
210//     }
211
212    return fUrgent;
213}
214
215//______________________________________________________________________________
216void OptUserStackingAction::NewStage(){
217
218    // get a new photon
219    if(!fNTotal)return;
220
221    Photon* photon = fPupil->GetPhoton();
222
223    // build an empty EDetPhoton in the event
224    if ( !photon || ( fNmax>0 && fNTracked>=fNmax ) ) {
225
226
227        if ( fEv ) {
228            TArrayI photonHistory(8);
229            photonHistory.AddAt(fNlost, 0);
230            photonHistory.AddAt(fNG4focalplane, kG4FocalSurface);
231            photonHistory.AddAt(fNfocalplane, kFocalPlane);
232            photonHistory.AddAt(fNoptics, kOptics);
233
234            EDetectorPhotonAdder a(photonHistory);
235            fEv->Fill(a);
236        }
237
238        MsgForm(EsafMsg::Info,"Number of events processed = %i.",fNTracked );
239        return;
240    }
241
242    forigin=0;
243    stackManager->PushOneTrack( NewPhoton(photon, fNTracked,0) );
244}
245
246
247//_______________________________________________________________________________
248void OptUserStackingAction::PrepareNewEvent(){
249
250    this->NewStage();
251}
252
253
254//______________________________________________________________________________
255void OptUserStackingAction::SetPupil(PhotonsOnPupil* p){
256    printf("**********Set event\n");
257    Reset();
258    fEv = EEvent::GetCurrent();
259    if(!p) return;
260
261   
262    fPupil=p;
263
264    fNTotal=fPupil->GetNphotons();
265   
266}
267
268//__________________________________________________________________
269void OptUserStackingAction::Reset(){
270
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   
287    fTimeFirstPhoton = fHuge;
288    fTimeLastPhoton = -fHuge;
289}
290
291//__________________________________________________________________
292EVector OptUserStackingAction::NewInteractionPoint(const Photon *p, double DZ) const {
293    EVector dir=p->dir.Unit();
294    EVector newPoint=p->pos+dir*(DZ/dir[Z]);
295#ifdef TRACEDEBUG
296    cout<<"new intPoint: "<<newPoint<<endl;
297#endif /* DEBUG */
298    return newPoint;
299}
300
301
302//-----------------------------------------------------------------
303void OptUserStackingAction::PostUserTrackingAction(const G4Track* fTrackPost){
304
305    if( fTrackPost->GetDefinition() == G4OpticalPhoton::OpticalPhotonDefinition()){
306        Photon* ph = (Photon*)this->UpdateCurrentPhoton((G4Track*)fTrackPost);
307
308        if ( this->GoToElectronics(ph) ){         
309          return;
310        }
311       
312        switch ( ph->history ) {
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]); 
334        }
335       
336        if( fEv  ) {
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   
346    if (fNTotal){
347      fNTracked++;
348      this->Progress(100.*fNTracked/fNTotal,"Photons");
349    }
350}
351//______________________________________________________________________________
352void OptUserStackingAction::Progress(int n,const char* p){
353    int prog =n;
354
355    if ( prog%10==0 && prog>=fProgress ) {
356        Msg(EsafMsg::Info).SetProgress(fProgress);             //fProgress
357        Msg(EsafMsg::Info) << Form("Tracking %s:  ",p) << MsgCount;
358        fProgress+=10;
359
360    }
361
362
363}
364//______________________________________________________________________________
365bool OptUserStackingAction::GoToElectronics(Photon* ph){
366
367    if ( ph->history!= kG4FocalSurface){
368//         Msg(EsafMsg::Debug) << Form("photon isn't on FocalPlane %d",ph->history) << MsgDispatch;
369//         Msg(EsafMsg::Debug) << "@ new photon " << MsgDispatch;
370        return false;
371    }
372    EVector intPoint;
373//     ph->history=kFocalPlane;
374
375    if(ftus){
376        ph->posOnIfs= ph->pos;
377        ph->hitIfs=kTRUE;
378    }
379    else fIdealFocalSurface->HitPosition(ph);
380
381    if(ph->hitIfs)ph->history=kFocalPlane;
382     
383    double DZ= fFocalPlane->Bottom()-ph->pos[Z];
384    // calculate impact point on the focal plane
385    intPoint=NewInteractionPoint(ph, DZ);
386    ph->time+=(ph->pos - intPoint).Mag()/Clight();
387    ph->pos=intPoint;
388//     printf("1new!\n");
389   
390    Photon* newp = fFocalPlane->Transport(ph);
391//     printf("0new!\n");
392    if ( !newp ) return false;
393//     printf("new!\n");
394    stackManager->PushOneTrack( this->NewPhoton(newp,newp->id,1) );
395    return true;
396}
397
398void OptUserStackingAction::EndOfEventAction(const G4Event* anEvent){
399    int n=G4RunManager::GetRunManager()->GetCurrentRun()->GetNumberOfEventToBeProcessed();
400    int curn=G4RunManager::GetRunManager()->GetCurrentRun()->GetNumberOfEvent();
401    this->Progress(100.*curn/n,"events");
402
403    if (  (curn+1)>=n ) {
404        if ( fEv ) {
405                TArrayI photonHistory(9);
406                photonHistory.AddAt(fNlost, 0);
407                photonHistory.AddAt(fNG4focalplane, kG4FocalSurface);
408                photonHistory.AddAt(fNfocalplane, kFocalPlane);
409                photonHistory.AddAt(fNoptics, kOptics);
410                EDetectorPhotonAdder a(photonHistory);
411                fEv->Fill(a);
412        }
413    }
414
415}
416
417void OptUserStackingAction::PreUserTrackingAction(const G4Track* fTrackPre){
418
419}
420void OptUserStackingAction::BeginOfEventAction(const G4Event* anEvent){
421}
Note: See TracBrowser for help on using the repository browser.