fiedl / hole-ice-study

This project aims to incorporate the effects of hole ice into the clsim photon propagation simulation of the icecube neutrino observatory.
Other
4 stars 2 forks source link

Make sure that plane wave photon positions are distributed properly #41

Open fiedl opened 6 years ago

fiedl commented 6 years ago

This is how it is currently done: https://github.com/fiedl/hole-ice-study/blob/master/scripts/AngularAcceptance/lib/GeneratePhotonFlashModule.py#L77

bildschirmfoto 2018-03-01 um 15 55 29

Looks good, but better double check.

Here is a working example: http://code.icecube.wisc.edu/projects/icecube/browser/IceCube/sandbox/gluesenkamp/slow-mp-gen/trunk/private/slow-mp-gen/I3SlowMpGen.cxx

#include "icetray/I3Frame.h"
#include "dataclasses/I3Double.h"
#include "dataclasses/physics/I3Particle.h"
#include "dataclasses/physics/I3MCTree.h"
#include "dataclasses/physics/I3MCTreeUtils.h"
#include "sim-services/I3GeoShifter.h"
#include "dataclasses/I3Position.h"
#include "dataclasses/I3Direction.h"
#include "dataclasses/I3Constants.h"
#include "phys-services/I3RandomService.h"
 
#include <cmath>
using namespace std;

I3Position I3SlowMpGen::RandomizePosition(I3RandomServicePtr random,
--
I3Direction dir, double diskDist,
double diskRad)
{
double z =  diskDist/I3Units::m;
double r(NAN);
double theta(NAN);
 
r = sqrt(random->Uniform(0,pow(diskRad/I3Units::m,2)));
theta = random->Uniform(0,2*I3Constants::pi);
 
log_debug("Chosen radius is %f (m)",r);
log_debug("Chosen theta is %f (rad)",theta);
 
if(isnan(r)\|\|isnan(theta))
{
log_fatal("Radius (%f) or theta (%f) not set in Randomize Position",r,theta);
}
 
double x = r*cos(theta);
double y = r*sin(theta);
log_debug("(x,y) position of particle before disk rotation is (%f,%f) in m",x,y);
 
/**
* Rotates Center of disk to:
* (diskDist,zenith,0) then
* (diskDist,zenith,azimuth)
* in spherical coordinates
*/
I3Position p(x*I3Units::m,y*I3Units::m,z*I3Units::m);
p.RotateY(dir.GetZenith());
p.RotateZ(dir.GetAzimuth());
log_debug("Final (x,y,z) positon is (%f,%f,%f) in m",p.GetX(),p.GetY(),p.GetZ());
 
return p;
}