---

MuonPropagator Class Reference

MuonPropogator handles muon transit for ReactorEvent. More...

#include <MuonPropagator.hh>

List of all members.

Public Methods
	MuonPropagator (double newRmaxgen=RMAXGEN, double newDepth=450., double newEnergy=0., double newCosTheta=1.)
	MuonPropagator constructor. More...
	MuonPropagator (int test=0, double newRmaxgen=RMAXGEN, double newEnergy=0., double newX=0., double newY=0., double newZ=0., double newVx=0., double newVy=0., double newVz=0.)
	~MuonPropagator ()
	MuonPropagator destructor.
	MuonPropagator (const MuonPropagator &Muon)
	MuonPropagator copy constructor.
MuonPropagator &	operator= (const MuonPropagator &rhs)
	MuonPropagator overloaded = operator.
double *	Intersection (double, double, double, double, double, double, double, double, double, double)
double *	Rotate (char, double, double, double, double)
double	RadCorr (double, double)
double	GetMuonE () const
double	GetMuonCosTheta () const
double	GetMuonPhi () const
int	GetNumNeutron () const
	Returns the neutron multiplicity.
double	GetNeutronKE (int &n) const
double	GetNeutronTime (int &n) const
double *const	GetNeutronVertexXYZ ()
int	GetNumElectron () const
	Returns the electron multiplicity.
double	GetElectronKE (int &n) const
double	GetElectronTime (int &n) const
double *const	GetElectronVertexXYZ ()
int	GetNumGamma () const
	Returns the number of photons.
double	GetGammaKE (int &n) const
double	GetGammaTime (int &n) const
double *const	GetGammaVertexXYZ ()
double	GetWeight () const
Static Public Methods
float	IsotopeDistance (const float &r)
	The lateral activation profile of muons passing through the detector to produce radioactive isotopes, which describes the probability that an isotope is produced at a given distance r (in cm) perpendicular to the muon track.

---

Detailed Description

MuonPropogator handles muon transit for ReactorEvent.

It returns particle information from the muon itself, as well as any byproducts of the muon shower, including neutrons and electrons from 9Li/8He beta decay.

Definition at line 17 of file MuonPropagator.hh.

---

Constructor & Destructor Documentation

MuonPropagator::MuonPropagator (  double    newRmaxgen = RMAXGEN, double    newDepth = 450., double    newEnergy = 0., double    newCosTheta = 1. )

MuonPropagator constructor. Generates a muon from a given energy and zenith angle. The muon starts at a random point above the detector and at a random azimuthal angle to decide the path of the muon track. Energy from the muon inside the detector is created incrementally along the muon tragectory as optical photons. One heavy isotope 9Li or 8He is produced per muon from the muon spallation and forced to beta decay to create and electron and neutron. The events are weighted by the 9Li/8He production cross sections and the decay branching ratio of the 9Li/8He beta decay. The neutron, electron, and photon energies, positions and times are saved. Definition at line 402 of file MuonPropagator.cpp. References MyHeSource::GetHeBetaEnergy, MyHeSource::GetHeDecayTime, MyHeSource::GetHeNeutronEnergy, MyLiSource::GetLiBetaEnergy, MyLiSource::GetLiDecayTime, MyLiSource::GetLiNeutronEnergy, and IsotopeDistance. { Rmaxgen = newRmaxgen; Depth = newDepth; energy = newEnergy; cosg = newCosTheta; //cout << "MuonPropagator::MuonPropagator(" << Rmaxgen << "," // << Depth << "," << energy << "," << cosg << ")" << endl; static double re = RELECTRON; static double me = MELECTRON; static double mmu = MMUON; static double Na = NAVOGADRO; // // For electron density // static double nel[2]; static double A[2] = {A0,A1}; static double Z[2] = {Z0,Z1}; static double f[2] = {f0,f1}; static double I[2] = {I0,I1}; // static double ncar; // // calculate the fluxes and densities // bool static first = true; double flux; int static muonstop; float static fspace[200]; // if(first){ // flux information if(Depth == 300.) flux = 0.374; // per m**2 per sec else if(Depth == 450.) flux = 0.160; // per m**2 per sec else if(Depth == 500.) flux = 0.145; // per m**2 per sec else{ cout << " ERROR: no muon flux for depth of " << Depth << " mwe" << endl; return; } // electron density for(int i=0;i<2;i++) nel[i] = Na*density*f[i]*Z[i]/A[i]; //cout << " nel " << nel[0]+nel[1] << endl; // 12C atom density ncar = Na*density*f[1]/A[1]; //cout << " ncar " << ncar << endl; // prepare the lateral activation profile float xlo = 0.; float xhi = 50.; funlxp_(IsotopeDistance,fspace,xlo,xhi); Weight = 0.; muonstop = 0; first = false; } /* pick a starting point: define a square above the detector through which the muons must pass. The square is centered on the z-axis at Rmaxgen and the area is A = (Rmaxgen+delta)**2. Randomly pick a point on the square to start the muon and use cos(theta) and a randomly selected azimuthal angle to decide how the muon hits the detector. Only keep events where the muon hits the scint. */ int inter = 0; double r[6],t[3]; double x[2],y[2],z[2]; double delta = Rmaxgen*20.; bool passed = false; // while(!passed){ const int len=3; float rv[len]; ranlux_(rv,len); r[0] = (2*rv[0]-1)*delta; r[1] = (2*rv[1]-1)*delta; r[2] = Rmaxgen; phig = 2*PI*rv[2]; // pick the direction t[0] = cos(phig)*sqrt(1-cosg*cosg); t[1] = sin(phig)*sqrt(1-cosg*cosg); t[2] = cosg; //for(int i=0;i<3;i++) cout << "t["<<i<<"] " << t[i] << " "; //cout << endl; // find the far end of the line a unit distance // from (x,y,z) along the t direction for(int i=0;i<3;i++) r[i+3] = r[i]-t[i]; // //for(int i=0;i<6;i++) cout << "r["<<i<<"] " << r[i] << " "; //cout << endl; // solve for intersections with the scint region double p[7]; double* ptr = p; int num; ptr = Intersection(r[0], r[1], r[2], // first end of line r[3], r[4], r[5], // second end of line 0, 0, 0, Rmaxgen); // detector centered at 0 num = (int)*ptr++; if(num>1){ inter = 2; for(int i=0; i<inter; i++){ x[i] = *ptr++; y[i] = *ptr++; z[i] = *ptr++; } passed = true; } } /* for(int i=0; i<inter; i++){ cout << " X-coordinate of point " << i << " : " << x[i] << endl; cout << " Y-coordinate of point " << i << " : " << y[i] << endl; cout << " Z-coordinate of point " << i << " : " << z[i] << endl; } */ // path length double L = sqrt((x[1]-x[0])*(x[1]-x[0])+ (y[1]-y[0])*(y[1]-y[0])+ (z[1]-z[0])*(z[1]-z[0])); //cout << " path length " << L << endl; /* iterate through the muon energy loss dE/dx in dx cm steps, making a gamma at the center of the dx cm segment of the line along the muon path equal to the energy lost in that segment; always check that the muon is still in the detector and has not stopped */ double step=5.0; double Emin=0.01; // double range = 0.; double time = 0.; // // define energy and position arrays // int maxgammafromisotope = 5; int maxgamma = int(Rmaxgen*2/step)+maxgammafromisotope; double* egamma = new double[maxgamma]; double* rgamma = new double[maxgamma*3]; double* tgamma = new double[maxgamma]; // double* egptr=egamma; double* rgptr=rgamma; double* tgptr=tgamma; // // Begin energy loss loop // double dE=0; double eloss; int ngamma = 0; // bool done = false; while(!done){ // zero the energy loss for this segment eloss = 0; // // check that the muon hasn't stopped // if(energy-dE > Emin){ // // check that we aren't leaving the detector // if(range+step < L){ // // iterate range and time // double gamma = 1+(energy-dE)/mmu; double beta = sqrt(1-1/gamma/gamma); range+=step; time+=step/beta/c; // // increment energy loss // double te = gamma-1; double tc = Emin/mmu; double tmax = te; double tup = tc; if(tc>tmax)tup=tmax; // for (int k=0;k<2;k++){ eloss+=2*PI*re*re*me*nel[k]*step/beta/beta* (log(2*(gamma+1)/(I[k]*I[k]/me/me)) +RadCorr(te,tup)); } // // add the energy loss to dE and set the "photon" energy // if(dE+eloss > energy){ *egptr++=energy-dE; // the muon stops, limit dE and egamma dE=energy; } else{ *egptr++=eloss; // the muon is still travelling dE+=eloss; } // // put the photon at the center of the segment // *rgptr++ = x[1] - (range-(step/2))*t[0]; *rgptr++ = y[1] - (range-(step/2))*t[1]; *rgptr++ = z[1] - (range-(step/2))*t[2]; // *tgptr++ = time; // // count the photon ngamma++; } else{ /* we have reached the end of the detector; put one final gamma to cover the last remaining segment of energy loss */ double laststep = L-range; // double gamma = 1+(energy-dE)/mmu; double beta = sqrt(1-1/gamma/gamma); range+=laststep; time+=laststep/beta/c; // // increment energy loss // double te = gamma-1; double tc = Emin/mmu; double tmax = te; double tup = tc; if(tc>tmax)tup=tmax; // for (int k=0;k<2;k++){ eloss+=2*PI*re*re*me*nel[k]*laststep/beta/beta* (log(2*(gamma+1)/(I[k]*I[k]/me/me)) +RadCorr(te,tup)); } // // add the energy loss to dE and set the "photon" energy // if(dE+eloss > energy){ *egptr++=energy-dE; // the muon stops, limit dE and egamma dE=energy; } else{ *egptr++=eloss; // the muon is still travelling dE+=eloss; } // // put the photon at the center of the segment // *rgptr++ = x[1] - (range-(laststep/2))*t[0]; *rgptr++ = y[1] - (range-(laststep/2))*t[1]; *rgptr++ = z[1] - (range-(laststep/2))*t[2]; // *tgptr++ = time; // // count the last photon ngamma++; // // finish the loop done = true; } } else{ /* the muon has stopped in the detector, make one last photon with energy equal the mass of the muon */ muonstop++; // set the "photon" energy // *egptr++=mmu; // // put the photon at point the muon stopped // *rgptr++ = x[1] - range*t[0]; *rgptr++ = y[1] - range*t[1]; *rgptr++ = z[1] - range*t[2]; // *tgptr++ = time; // // count the last photon ngamma++; // // finish the loop done = true; /* cout << "Muon " << muonstop << " has stopped at " << x[1] - range*t[0] << ", " << y[1] - range*t[1] << ", " << z[1] - range*t[2] << " in " << time << " ns" << endl; */ } } //cout << " muon range " << range << " cm, time " << time // << " ns and energy loss " << dE << " MeV" << endl; /* now deal with isotope production; start with the cross-section for 9Li+8He prodution (as a function of incident muon energy in MeV) and generate an isotope; this isn't perfect because the muon changes energy throughout its path, altering the XC, but it's okay for now */ double XC; double mfp; double alpha = 0.73; // from the Hagner paper double XCmeas = 2.12; double Emeas = 190000.; // 190 GeV double ubarntocmsq = 1e-30; //cout << " sf " << (XCmeas/pow(Emeas,alpha))*ubarntocmsq << endl; // XC = (XCmeas/pow(Emeas,alpha))*pow(energy,alpha)*ubarntocmsq; //cout << " XC " << XC << " cm^2" << endl; // // compute mean free path to produce 9Li+8He // int ntrial = 0; passed = false; while(!passed){ mfp = 1/ncar/XC; // // skip if ranlux doesn't have the necessary precision // double test = 1/pow(10.0,(range/mfp)); if(test > 0.9999998){ mfp = range/2.; ntrial = 0; } else{ const int len=1; float rv[len]; ranlux_(rv,len); mfp*=log(1/rv[0]); } // check that the isotope is inside the detector if(mfp <= range) passed = true; // count the attempt ntrial++; } //cout << " mfp " << mfp << " cm in " << ntrial << " trials" << endl; /* start at the mean free path distance from where the muon entered the detector along the muon track in the direction of the muon flight to produce the isotope; define a new coordinate system with the muon track being the +z'-axis, find the perpendicular distance from the muon track via Hagner, and pick a random azimuthal angle; thus the perpendicular distance and azimuthal angle give the new x' and y' values for the isotope and the point of isotope production along the muon track is z' = 0; then rotate the new coordinate system into the detector x, y, and z */ double ri[3]; ri[0] = x[1] - mfp*t[0]; ri[1] = y[1] - mfp*t[1]; ri[2] = z[1] - mfp*t[2]; // // pick the direction, relative to point ri // const int len=2; float rw[len]; ranlux_(rw,len); double newphig = 2*PI*rw[0]; double newcosg = 0.; double newt[3]; newt[0] = cos(newphig)*sqrt(1-newcosg*newcosg); newt[1] = sin(newphig)*sqrt(1-newcosg*newcosg); newt[2] = newcosg; // // select the lateral distance from the muon track, relative to ri // float lateraldistance; int len2 = 1; funlux_(fspace,lateraldistance,len2); // // define point rt by x',y', and z' in the new coordinate system // double rt[3]; rt[0] = lateraldistance*newt[0]; rt[1] = lateraldistance*newt[1]; rt[2] = lateraldistance*newt[2]; //for(int i=0; i<3; i++)cout << " rt[" << i << "] " << rt[i]; //cout << endl; /* always assume that the x'-axis of the new coordinate system is in the direction chosen by angle phig of the muon track above; thus we rotate point rt around the y'-axis (which is now chosen always perpendicular to the detector coordinate system's x-y plane) by the angle given by cosg above subtracted from pi; we then rotate the new point qt around the z'-axis by angle pi minus phig to get the new point pt, which is rt in the detector coordinate system */ double angle = PI - acos(cosg); double qt[3]; double* qtptr = qt; qtptr = Rotate('Y',angle,rt[0],rt[1],rt[2]); for(int i=0; i<3; i++) qt[i] = *qtptr++; // // second rotation // angle = PI - phig; double pt[3]; double* ptptr = pt; ptptr = Rotate('Z',angle,qt[0],qt[1],qt[2]); for(int i=0; i<3; i++) pt[i] = *ptptr++; // // add the new lateral displacement to the point on the muon // track at which the isotope is produced // for(int i=0; i<3; i++)ri[i] += pt[i]; //for(int i=0; i<3; i++)cout << " ri[" << i << "] " << ri[i]; //cout << endl; // // decide if we have 9Li or 8He // char* isotope; if(rw[1] <= 0.5) isotope = "9Li"; else isotope = "8He"; // // make energy and position arrays // int maxneutron = 1; double* eneutron = new double[maxneutron]; double* rneutron = new double[maxneutron*3]; double* tneutron = new double[maxneutron]; // double* enptr = eneutron; double* rnptr = rneutron; double* tnptr = tneutron; // int maxelectron = 1; double* eelectron = new double[maxelectron]; double* relectron = new double[maxelectron*3]; double* telectron = new double[maxelectron]; // double* eeptr = eelectron; double* reptr = relectron; double* teptr = telectron; // // let the isotope decay // int nneutron = 0; int nelectron = 0; // if(isotope == "9Li"){ MyLiSource LiSource(9,Rmaxgen); nneutron = 1; *enptr++ = LiSource.GetLiNeutronEnergy(); for(int i=0; i<3; i++) *rnptr++ = ri[i]; *tnptr++ = LiSource.GetLiDecayTime(); nelectron = 1; *eeptr++ = LiSource.GetLiBetaEnergy(); for(int i=0; i<3; i++) *reptr++ = ri[i]; *teptr++ = LiSource.GetLiDecayTime(); // Weight = 1./float(ntrial)*LiSource.GetLiWeight(); Weight = 1./float(ntrial); } else if(isotope == "8He"){ MyHeSource HeSource(8,Rmaxgen); nneutron = 1; *enptr++ = HeSource.GetHeNeutronEnergy(); for(int i=0; i<3; i++) *rnptr++ = ri[i]; *tnptr++ = HeSource.GetHeDecayTime(); nelectron = 1; *eeptr++ = HeSource.GetHeBetaEnergy(); for(int i=0; i<3; i++) *reptr++ = ri[i]; *teptr++ = HeSource.GetHeDecayTime(); // Weight = 1./float(ntrial)*HeSource.GetHeWeight(); Weight = 1./float(ntrial); } else{ cout << " ERROR: unknown isotope " << isotope << ". Options are '9Li' and '8He'." << endl; } // // transfer temp data to storage // Nneutron = nneutron; Nelectron = nelectron; Ngamma = ngamma; Eneutron = new double[nneutron]; Rneutron = new double[nneutron*3]; Tneutron = new double[nneutron]; Eelectron = new double[nelectron]; Relectron = new double[nelectron*3]; Telectron = new double[nelectron]; Egamma = new double[ngamma]; Rgamma = new double[ngamma*3]; Tgamma = new double[ngamma]; double* ptr; ptr = Eneutron+nneutron; while(ptr>Eneutron)*--ptr=*--enptr; ptr = Rneutron+3*nneutron; while(ptr>Rneutron)*--ptr=*--rnptr; ptr = Tneutron+nneutron; while(ptr>Tneutron)*--ptr=*--tnptr; ptr = Eelectron+nelectron; while(ptr>Eelectron)*--ptr=*--eeptr; ptr = Relectron+3*nelectron; while(ptr>Relectron)*--ptr=*--reptr; ptr = Telectron+nelectron; while(ptr>Telectron)*--ptr=*--teptr; ptr = Egamma+ngamma; while(ptr>Egamma)*--ptr=*--egptr; ptr = Rgamma+3*ngamma; while(ptr>Rgamma)*--ptr=*--rgptr; ptr = Tgamma+ngamma; while(ptr>Tgamma)*--ptr=*--tgptr; /* cout << Nneutron << " neutrons" << endl; for(int n=0; n<Nneutron; n++){ cout << " neutron " << n << " energy " << Eneutron[n] << " MeV, and x, y, z:"; for(int i=0; i<3; i++) cout << " " << Rneutron[3*n+i]; cout << endl << " and time " << Tneutron[n] << " nsec" << endl; } cout << Nelectron << " electrons" << endl; for(int n=0; n<Nelectron; n++){ cout << " electron " << n << " energy " << Eelectron[n] << " MeV, and x, y, z:"; for(int i=0; i<3; i++) cout << " " << Relectron[3*n+i]; cout << endl << " and time " << Telectron[n] << " nsec" << endl; } cout << Ngamma << " photons" << endl; for(int n=0; n<Ngamma; n++){ cout << " photon " << n << " energy " << Egamma[n] << " MeV, and x, y, z:"; for(int i=0; i<3; i++) cout << " " << Rgamma[3*n+i]; cout << endl << " and time " << Tgamma[n] << " nsec" << endl; } */ // clean up delete[] eneutron; delete[] rneutron; delete[] tneutron; delete[] eelectron; delete[] relectron; delete[] telectron; delete[] egamma; delete[] rgamma; delete[] tgamma; }

---

The documentation for this class was generated from the following files:

• MuonPropagator.hh
• MuonPropagator.cpp

---