source: JEM-EUSO/esaf_lal/tags/v1_r0/esaf/packages/simulation/lighttoeuso/src/TestLightToEuso.cc @ 117

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

ESAF version compilable on mac OS
File size: 18.1 KB
Line 
1// ESAF : Euso Simulation and Analysis Framework
2// $Id: TestLightToEuso.cc 2828 2009-04-01 14:29:16Z biktem $
3// M. Pallavicini created Jul,  2 2002
4
5#include "TestLightToEuso.hh"
6#include "EsafRandom.hh"
7#include "ListPhotonsOnPupil.hh"
8#include "utils.hh"
9#include "MCTruth.hh"
10#include "EEventTruthAdder.hh"
11#include "EEvent.hh"
12#include "EConst.hh"
13#include "DetectorGeometry.hh"
14#include "EsafSpectrum.hh"
15
16#include "TMath.h"
17
18using namespace sou;
19using namespace TMath;
20using namespace EConst;
21
22ClassImp(TestLightToEuso)
23
24//______________________________________________________________________________
25//
26// This class brings test fPhotons to the Euso Pupil
27// The type of events supported are the following:
28//     SPOT         light from a single point at the same time
29//     EXTSPOT      light from a single point at the same time
30//     RANDSPOT     single point randomly chosen
31//     CIRCLE       light making a circle in focal surface (theta constant,
32//                  phi 0-360 degrees)
33//     EXTCIRCLE    light making a thick circle in focal surface (theta
34//                  constant, phi 0-360 degrees)
35//     RADIUS       light making a radius (phi constant, theta 0-30 degrees)
36//     TRACK        light making a track in focal surface with peaked pattern
37//     SHOWERTRACK  track in the atmosphere
38//
39// The parameters needed are the following:
40//     TestLightToEuso.Type       see above
41//     TestLightToEuso.Duration   event duration in micro-seconds
42//     TestLightToEuso.Photons    number of fPhotons generated per event
43//     TestLightToEuso.Theta1     constant if CIRCLE, SPOT or EXTSPOT,
44//                                starting point if EXTCIRCLE, RADIUS or TRACK
45//     TestLightToEuso.Theta2     ignored if CIRCLE, SPOT or EXTSPOT,
46//                                end point if RADIU, EXTCIRCLE or TRACK
47//     TestLightToEuso.Phi1       starting point if CIRCLE, EXTCIRCLE or TRACK,
48//                                constant if RADIUS, SPOT or EXTSPOT
49//     TestLightToEuso.Phi2       end point if CIRCLE, EXTCIRCLE or TRACK,
50//                                ignored if RADIUS
51//                    or SPOT, radius of the beam in EXTSPOT
52//
53
54//______________________________________________________________________________
55TestLightToEuso::TestLightToEuso() : LightToEuso("TEST"), fTruth(NULL) {
56    //
57    // Constructor
58    //
59
60    fPhOnPupil = new ListPhotonsOnPupil();//(&fPhotons);
61}
62
63//______________________________________________________________________________
64TestLightToEuso::~TestLightToEuso() {
65    //
66    // Destructor
67    //
68
69    Reset();
70    delete fPhOnPupil;
71    if ( fTruth ) delete fTruth;
72}
73
74//______________________________________________________________________________
75Bool_t TestLightToEuso::ClearMemory() {
76    //
77    // Release all the memory hold in the arrays
78    //
79
80    return fPhOnPupil ? fPhOnPupil->ClearMemory() : kTRUE;
81}
82
83//______________________________________________________________________________
84void TestLightToEuso::Reset() {
85    //
86    // Reset internal list of fPhotons
87    //
88
89    fPhOnPupil->Clear();
90}
91
92//______________________________________________________________________________
93PhotonsOnPupil *TestLightToEuso::Get(const DetectorGeometry* dg) {
94    //
95    // Returns the list of fPhotons for the pupil.
96    // The type of photon distribution is determined
97    // by TestLightToEuso.Type
98    //
99
100    Reset();
101    string type = Conf()->GetStr("TestLightToEuso.Type");
102    Double_t th1 = TMath::DegToRad()*Conf()->GetNum("TestLightToEuso.Theta1");
103    Double_t th2 = TMath::DegToRad()*Conf()->GetNum("TestLightToEuso.Theta2");
104    Double_t ph1 = TMath::DegToRad()*Conf()->GetNum("TestLightToEuso.Phi1");
105    Double_t ph2 = TMath::DegToRad()*Conf()->GetNum("TestLightToEuso.Phi2");
106    Int_t nPhotons = (Int_t)Conf()->GetNum("TestLightToEuso.Photons");
107    Double_t tm = Conf()->GetNum("TestLightToEuso.Duration");
108    TRandom* rndm = EsafRandom::Get();
109
110    if ( type == "SPOT" ) {              // spot in fixed position
111        return MakeSpot(th1,ph1,nPhotons,tm);
112    }
113    else if ( type == "RANDSPOT" ) {     // spot in random position within Phi1,2 Theta1,2
114        Double_t cmax = Cos(th2);
115        Double_t cmin = Cos(th1);
116        if ( th1 > th2 ) {
117            Msg(EsafMsg::Info) << "In TestLightToEuso Theta1 and Theta2 will be switched\n";
118            cmax = cmin;
119            cmin = Cos(th2);
120        }
121        Double_t diff = cmax - cmin;
122        Double_t xc = cmin + rndm->Rndm()*diff;
123        Double_t th = acos(xc);
124        if ( ph2 > ph1 )
125            diff = ph2 - ph1;
126        else {
127            Msg(EsafMsg::Info) << "In TestLightToEuso Phi1 and Phi2 will be switched\n";
128            diff = ph1 - ph2;
129            Double_t temp = ph2;
130            ph2 = ph1;
131            ph1 = temp;
132        }
133        Double_t ph = rndm->Rndm()*diff + ph1;
134        return MakeSpot(th,ph,nPhotons,tm);
135    }
136    else if ( type == "EXTSPOT" ) {
137        return MakeExtSpot(th1,ph1,RadToDeg()*ph2,nPhotons,tm);
138    }
139    else if ( type == "CIRCLE" ) {
140        return MakeCircle(ph1,ph2,th1,nPhotons,tm);
141    }
142    else if ( type == "EXTCIRCLE" ) {
143        return MakeExtCircle(ph1,ph2,th1,th2,dg->GetRadius(),nPhotons,tm);
144    }
145    else if ( type == "DIFFUSE" ) {
146        return MakeDiffuse(th1,th2,TMath::RadToDeg()*ph1,nPhotons,tm);
147    }
148    else if ( type == "CFLUXDIFFUSE" ) {
149        return MakeCFluxDiffuse(th1,th2,TMath::RadToDeg()*ph1,nPhotons,tm);
150    }
151    else if ( type == "RADIUS" ) {
152        return MakeRadius(ph1,th1,th2,nPhotons,tm);
153    }
154    else if ( type == "EXTRADIUS" ) {
155        return MakeExtRadius(ph1,th1,th2,dg->GetRadius(),nPhotons,tm);
156    }
157    else if ( type == "TRACK" ) {
158        return MakeTrack(ph1,ph2,th1,th2,nPhotons,tm);
159    }
160    else if ( type == "SHOWERTRACK" ) {
161
162        return MakeShowerTrack( dg );
163    }
164    throw runtime_error("Invalid parameter TestLightToEuso.Type = " + type );
165    return (PhotonsOnPupil*)0;
166}
167
168//______________________________________________________________________________
169PhotonsOnPupil *TestLightToEuso::MakeSpot( Double_t th,  Double_t  ph, Int_t nPhotons,
170        Double_t  tm) {
171    //
172    // Produce light in a single spot
173    //
174
175    Msg(EsafMsg::Info) << "SPOT CALLED th="<<th<<" ph="<<ph<<" Photons="<<nPhotons<< MsgDispatch;
176    Double_t wl = 357.*nm;
177    Double_t x[3] = {0.,0.,0.};
178    Double_t y[3] = {0.,0.,0.};
179    Double_t t;
180    TRandom* rndm = EsafRandom::Get();
181    for(Int_t i=0; i<nPhotons; i++) {
182       t = rndm->Rndm()*tm*microsecond;
183       ParentPhoton *p = new ParentPhoton(i,x,y,th,ph,wl,t,Fluorescence);
184       fPhOnPupil->Add( p );
185    }
186    return fPhOnPupil;
187}
188
189//______________________________________________________________________________
190PhotonsOnPupil *TestLightToEuso::MakeExtSpot( Double_t theta, Double_t phi,
191        Double_t radius, Int_t nPhotons, Double_t tm) {
192    //
193    // Produce light in an extensive spot
194    //
195
196    Msg(EsafMsg::Info) << "EXTSPOT CALLED theta="<<theta<<" phi="<<phi<<" radius="<<radius<<" Photons="<<nPhotons<< MsgDispatch;
197    Double_t wl = 357.*nm;
198    wl = (Conf()->GetNum("TestLightToEuso.wl"))*nm;
199    Msg(EsafMsg::Debug)<<Form("Generating %i photons with wl=%f and theta=%.2f (extspot)",nPhotons,wl,theta*TMath::RadToDeg())<<MsgDispatch;
200
201    Double_t x[3] = {0.,0.,0.};
202    Double_t y[3] = {0.,0.,0.};
203    Double_t t;
204    TRandom* rndm = EsafRandom::Get();
205    for(Int_t i=0; i<nPhotons; i++) {
206        t = rndm->Rndm()*tm*microsecond;
207        Double_t radius1 = Sqrt(rndm->Rndm())*radius;
208        Double_t th2 = rndm->Rndm()*TMath::TwoPi();
209        y[0] = radius1*Cos(th2)*mm;
210        y[1] = radius1*Sin(th2)*mm;
211        ParentPhoton *p = new ParentPhoton(i,x,y,theta,phi,wl,t,Fluorescence);
212        fPhOnPupil->Add( p );
213    }
214    return fPhOnPupil;
215}
216
217//______________________________________________________________________________
218PhotonsOnPupil *TestLightToEuso::MakeCircle( Double_t ph1,  Double_t  ph2,  Double_t  th, Int_t nPhotons,  Double_t tm) {
219    //
220    // Produce light in a circle
221    //
222    TRandom* rndm = EsafRandom::Get();
223    Msg(EsafMsg::Info) << "CIRCLE CALLED with ph1=" << ph1 << " ph2=" << ph2 << " th=" << th << MsgDispatch;
224    Double_t wl = 357.*nm;
225    Double_t x[3] = {0.,0.,0.};
226    Double_t y[3] = {0.,0.,0.};
227    Double_t t;
228    for(Int_t i=0; i<nPhotons; i++) {
229        t = rndm->Rndm()*tm*microsecond;
230        Double_t ph = ph1 + rndm->Rndm()*(ph2-ph1);
231        ParentPhoton *p = new ParentPhoton(i,x,y,th,ph,wl,t,Fluorescence);
232        fPhOnPupil->Add( p );
233    }
234    return fPhOnPupil;
235}
236
237//______________________________________________________________________________
238PhotonsOnPupil *TestLightToEuso::MakeExtCircle( Double_t ph1,  Double_t ph2,
239        Double_t  th1,  Double_t th2, Double_t radius, Int_t nPhotons,
240        Double_t  tm) {
241    //
242    // Produce light in a extensive circle (fPhotons with th1 < th < th2)
243    //
244
245    TRandom *rndm = EsafRandom::Get();
246    Msg(EsafMsg::Info) << "EXTCIRCLE CALLED with ph1=" << ph1 << " ph2=" << ph2 <<
247        " th1=" << th1 << " th2=" << th2 << MsgDispatch;
248    Double_t wl = 357.*nm;
249    Double_t x[3] = {0.,0.,0.};
250    Double_t y[3] = {0.,0.,0.};
251    Double_t t,theta1,theta2;
252    if ( th1 > th2 ) {
253        theta1=th2;
254        theta2=th1;
255    } else {
256        theta1=th1;
257        theta2=th2;
258    }
259    for(Int_t i=0; i<nPhotons; i++) {
260        t = rndm->Rndm()*tm*microsecond;
261        Double_t th = theta1+(theta2-theta1)*Sqrt((theta2-theta1)/(Pi()/6)*rndm->Rndm());
262        Double_t ph = ph1 + rndm->Rndm()*(ph2-ph1);
263        ParentPhoton *p = new ParentPhoton(i,x,y,th,ph,wl,t,Fluorescence);
264        fPhOnPupil->Add( p );
265    }
266    return fPhOnPupil;
267
268}
269
270//______________________________________________________________________________
271PhotonsOnPupil *TestLightToEuso::MakeDiffuse( Double_t th1,  Double_t  th2,
272         Double_t radius, Int_t nPhotons,  Double_t  tm) {
273    //
274    //
275    //
276
277    TRandom *rndm = EsafRandom::Get();
278    Msg(EsafMsg::Info) << "DIFFUSE CALLED with th1=" << th1 << " th2=" << th2 <<
279        " radius=" << radius <<  MsgDispatch;
280    Double_t wl = 357.*nm;
281    Double_t x[3] = {0.,0.,0.};
282    Double_t y[3] = {0.,0.,0.};
283    Double_t t(0),theta1,theta2;
284    if ( th1 > th2 ) {
285        theta1=th2;
286        theta2=th1;
287    } else {
288        theta1=th1;
289        theta2=th2;
290    }
291
292    Double_t c1 = Cos(theta1);
293    Double_t c2 = Cos(theta2);
294    for(Int_t i=0; i<nPhotons; i++) {
295//        Double_t th = theta1+(theta2-theta1)*Sqrt((theta2-theta1)/(Pi()/6)*rndm->Rndm());
296        Double_t th = ACos(rndm->Uniform(c1,c2));
297        Double_t ph = TwoPi()*rndm->Rndm();
298
299        Double_t rho = Sqrt(rndm->Rndm())*radius;
300        Double_t alpha = rndm->Rndm()*TwoPi();
301        y[0] = rho*Cos(alpha)*mm;
302        y[1] = rho*Sin(alpha)*mm;
303
304        ParentPhoton *p = new ParentPhoton(i,x,y,th,ph,wl,t,Fluorescence);
305        fPhOnPupil->Add( p );
306    }
307    return fPhOnPupil;
308}
309
310//______________________________________________________________________________
311PhotonsOnPupil *TestLightToEuso::MakeCFluxDiffuse( Double_t th1,  Double_t
312        th2,  Double_t radius, Int_t nPhotons,  Double_t  tm) {
313    TRandom *rndm = EsafRandom::Get();
314    Msg(EsafMsg::Info) << "C FLUX-DIFFUSE CALLED with th1=" << th1 << " th2=" << th2 <<
315        " radius=" << radius << MsgDispatch;
316    Double_t wl = 357.*nm;
317    Double_t x[3] = {0.,0.,0.};
318    Double_t y[3] = {0.,0.,0.};
319    Double_t t(0),theta1(0),theta2(0);
320    if ( th1 > th2 ) {
321        theta1=th2;
322        theta2=th1;
323    } else {
324        theta1=th1;
325        theta2=th2;
326    }
327
328    Double_t c1 = Cos(2*theta1);
329    Double_t c2 = Cos(2*theta2);
330    for(Int_t i(0); i<nPhotons; i++) {
331        Double_t th = 0.5*ACos(rndm->Uniform(c1,c2));
332        Double_t ph = TMath::TwoPi()*rndm->Rndm();
333
334        Double_t rho = Sqrt(rndm->Rndm())*radius;
335        Double_t alpha = rndm->Rndm()*TMath::TwoPi();
336        y[0] = rho*Cos(alpha)*mm;
337        y[1] = rho*Sin(alpha)*mm;
338
339        ParentPhoton *p = new ParentPhoton(i,x,y,th,ph,wl,t,Fluorescence);
340        fPhOnPupil->Add( p );
341    }
342    return fPhOnPupil;
343}
344
345//______________________________________________________________________________
346PhotonsOnPupil *TestLightToEuso::MakeRadius( Double_t ph,  Double_t th1,
347        Double_t  th2, Int_t nPhotons,  Double_t  tm) {
348    //
349    // Produce light along a focal surface radius.
350    //
351
352    TRandom* rndm = EsafRandom::Get();
353    Msg(EsafMsg::Info) << "RADIUS CALLED" << MsgDispatch;
354    Double_t wl = 357.*nm;
355    Double_t x[3] = {0.,0.,0.};
356    Double_t y[3] = {0.,0.,0.};
357    Double_t t;
358    for(Int_t i=0; i<nPhotons; i++) {
359       t = rndm->Rndm()*tm*microsecond;
360       Double_t th = th1 + rndm->Rndm()*(th2-th1); // uniform in theta, not costheta
361       ParentPhoton *p = new ParentPhoton(i,x,y,th,ph,wl,t,Fluorescence);
362       fPhOnPupil->Add( p );
363    }
364    return fPhOnPupil;
365}
366
367//______________________________________________________________________________
368PhotonsOnPupil *TestLightToEuso::MakeExtRadius( Double_t ph,  Double_t th1,
369        Double_t  th2, Double_t radius, Int_t nPhotons, Double_t tm) {
370    //
371    // Produce light along a focal surface radius.
372    //
373
374    TRandom* rndm = EsafRandom::Get();
375    Msg(EsafMsg::Info) << "EXTRADIUS CALLED" << MsgDispatch;
376    Double_t wl = 357.*nm;
377    Double_t x[3] = {0.,0.,0.};
378    Double_t y[3] = {0.,0.,0.};
379    Double_t t;
380    for(Int_t i=0; i<nPhotons; i++) {
381        t = rndm->Rndm()*tm*microsecond;
382        Double_t th = th1 + rndm->Rndm()*(th2-th1); // uniform in theta, not costheta
383        Double_t rho = Sqrt(rndm->Rndm())*radius;
384        Double_t alpha = rndm->Rndm()*TwoPi();
385        y[0] = rho*Cos(alpha)*mm;
386        y[1] = rho*Sin(alpha)*mm;
387
388        ParentPhoton *p = new ParentPhoton(i,x,y,th,ph,wl,t,Fluorescence);
389        fPhOnPupil->Add( p );
390    }
391    return fPhOnPupil;
392}
393
394
395//______________________________________________________________________________
396PhotonsOnPupil *TestLightToEuso::MakeTrack( Double_t ph1,  Double_t  ph2,
397        Double_t th1,  Double_t  th2, Int_t nPhotons,  Double_t  tm) {
398    //
399    // Produce light along a track
400    //
401
402    //TRandom* rndm = EsafRandom::Get();
403    return 0;
404}
405
406//______________________________________________________________________________
407PhotonsOnPupil *TestLightToEuso::MakeShowerTrack(const DetectorGeometry* geo) {
408    // produce light along a track
409
410    Msg(EsafMsg::Info) << "SHOWERTRACK CALLED" << MsgDispatch;
411
412    TRandom* rndm = EsafRandom::Get();
413
414    Double_t iX(0),iY(0),iZ(0);
415    Double_t tLength = Conf()->GetNum("TestLightToEuso.ShowerTrack.fTrackLength")*km;
416    Double_t tTheta = Conf()->GetNum("TestLightToEuso.ShowerTrack.fTheta")*TMath::DegToRad();
417    Double_t tPhi = Conf()->GetNum("TestLightToEuso.ShowerTrack.fPhi")*TMath::DegToRad();
418    string sInitPos = Conf()->GetStr("TestLightToEuso.ShowerTrack.fPosInit");
419    Int_t tProfileCode = Nint(Conf()->GetNum("TestLightToEuso.ShowerTrack.fProfileCode"));
420
421    if ( sInitPos == "fixed") {
422        // fixed int point
423        iX = Conf()->GetNum("TestLightToEuso.ShowerTrack.fXInit")*km;
424        iY = Conf()->GetNum("TestLightToEuso.ShowerTrack.fYInit")*km;
425        iZ = Conf()->GetNum("TestLightToEuso.ShowerTrack.fZInit")*km;
426
427    } else if (sInitPos == "random") {
428        // random interacion point
429        Double_t rInitMax = Conf()->GetNum("TestLightToEuso.ShowerTrack.fRInitMax")*km;
430        Double_t radInit = TMath::Sqrt(rndm->Rndm())*rInitMax;
431        Double_t phiInit = rndm->Rndm()*TMath::TwoPi();
432        iX = radInit*TMath::Cos(phiInit);
433        iY = radInit*TMath::Sin(phiInit);
434        // Z init between 0 and 50 km
435        iZ = rndm->Rndm()*Conf()->GetNum("TestLightToEuso.ShowerTrack.fZInitMax")*km;
436
437    } else if (sInitPos=="zfixed") {
438        // Z coordinate of the interacion point is kept fixed
439        Double_t rInitMax = Conf()->GetNum("TestLightToEuso.ShowerTrack.fRInitMax")*km;
440        Double_t radInit = TMath::Sqrt(rndm->Rndm())*rInitMax;
441        Double_t phiInit = rndm->Rndm()*TMath::TwoPi();
442        Msg(EsafMsg::Info) << rInitMax << " " << phiInit << " " << radInit << MsgDispatch;
443        iX = radInit*TMath::Cos(phiInit);
444        iY = radInit*TMath::Sin(phiInit);
445        iZ = Conf()->GetNum("TestLightToEuso.ShowerTrack.fZInit")*km;
446
447    } else {
448        Msg(EsafMsg::Info) << "TestLightToEuso.cfg:"
449             << " wrong value in TestLightToEuso.ShowerTrack.fPosInit:"
450             << sInitPos << MsgDispatch;
451        Msg(EsafMsg::Panic) << "ESAF Error: Wrong value in TestLightToEuso.cfg" << MsgDispatch;
452    }
453
454    Int_t nPhotons = (Int_t)Conf()->GetNum("TestLightToEuso.ShowerTrack.fPhotons");
455    const TVector3& detPos = geo->GetPos();
456    Double_t detRad = geo->GetRadius();
457
458    TVector3 initPos(iX, iY, iZ);
459    TVector3 diffPos;
460    diffPos.SetMagThetaPhi(tLength,tTheta,tPhi);
461    diffPos *= -1;
462    TVector3 endPos = diffPos + initPos;
463    TVector3 dir = (endPos-initPos).Unit();
464    TVector3 dL = (endPos-initPos)*(1./nPhotons);
465
466    fTruth = new MCTruth;
467    fTruth->SetFirstInt(initPos, 0);
468    fTruth->SetThetaPhi(tTheta, tPhi);
469    fTruth->SetEarthImpact(endPos);
470    fTruth->SetShowerMax(0.5*(endPos-initPos)+initPos, 0);
471
472    if ( EEvent::GetCurrent() ) {
473        EEventTruthAdder a(fTruth);
474
475        EEvent::GetCurrent()->Fill( a );
476    }
477
478    EsafSpectrum *shSpec = new EsafSpectrum("nitro10km");
479
480    Double_t wl;
481    Double_t t;
482    Double_t x[3] = {0.,0.,0.};
483    Double_t y[3] = {0.,0.,0.};
484    Double_t l;
485
486    for(Int_t i(0); i<nPhotons; i++) {
487        // get a random number, gaussian distributed between 0 and 1 sigma
488        switch ( tProfileCode) {
489            case 1:
490                for(l = rndm->Gaus(0.5,1/4.); l < 0 || l > 1;l = rndm->Gaus(0.5,1/4.) );
491                break;
492            case 0:
493            default:
494                l = rndm->Rndm();
495        }
496
497        TVector3 shPos = initPos+diffPos*l;
498        TVector3 inDir = -((shPos-detPos).Unit());
499        Double_t radius = Sqrt(rndm->Rndm())*detRad;
500        Double_t phi  = rndm->Rndm()*TMath::TwoPi();
501
502        x[0] = (shPos-detPos)[0];
503        x[1] = (shPos-detPos)[1];
504        x[2] = (shPos-detPos)[2];
505
506        y[0] = radius*Cos(phi)*mm;
507        y[1] = radius*Sin(phi)*mm;
508        y[2] = 0*mm;
509
510        wl = shSpec->GetLambda();
511
512        t = (l*tLength+(shPos-detPos).Mag())/Clight();
513        ParentPhoton *p = new ParentPhoton(i,x,y,inDir.Theta(),inDir.Phi(),wl,t,Fluorescence);
514        fPhOnPupil->Add( p );
515    }
516
517
518    return fPhOnPupil;
519}
Note: See TracBrowser for help on using the repository browser.