---

MuonPropagator.cpp

Go to the documentation of this file.

#include <math.h>
#include <iostream.h>
#include <fstream.h>   // file I/O
#include <iomanip.h>   // format manipulation
#include <string>
#include <stdio.h>
#include "ReactorConstants.hh"
#include "ReactorFortran.hh"
#include "MuonPropagator.hh"
#include "MyHeSource.hh"
#include "MyLiSource.hh"

MuonPropagator::MuonPropagator(int test, double newRmaxgen, double newEnergy,
                               double newX, double newY, double newZ,
                               double newVx, double newVy, double newVz){

  Rmaxgen = newRmaxgen;
  energy = newEnergy;

  double ri[3];
  ri[0] = newX;
  ri[1] = newY;
  ri[2] = newZ;

  double pi[3];
  pi[0] = newVx;
  pi[1] = newVy;
  pi[2] = newVz;

  /*
  cout << "MuonPropagator::MuonPropagator(" 
       << Rmaxgen << ","  << energy << "," 
       << ri[0] << "," << ri[1] << "," << ri[2] << ","
       << pi[0] << "," << pi[1] << "," << pi[2] << ")" << 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
  //
  int static muonstop;
  bool static first = true;
  if(first){

    // 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;
      
    Weight = 0.;
    muonstop = 0;
  
    first = false;
  }
  //
  // find where the muon enters and leaves the detector
  //
  int inter = 0;
  double x[2],y[2],z[2];

  // create the unit vector along the direction of momentum
  //
  double pdotp = sqrt(pi[0]*pi[0] + pi[1]*pi[1] + pi[2]*pi[2]);
  double t[3];
  for(int i=0;i<3;i++) t[i] = pi[i]/pdotp;
  //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
  //
  double rf[3];
  for(int i=0;i<3;i++) rf[i] = ri[i]-t[i];
  //for(int i=0;i<3;i++) cout << " rf[" << i << "] " << rf[i] << " ";
  //cout << endl;

  // solve for intersections with the scint region
  double p[7];
  double* pptr = p;
  int num;

  pptr = Intersection(ri[0], ri[1], ri[2],  // first end of line
                      rf[0], rf[1], rf[2],  // second end of line
                      0, 0, 0, Rmaxgen);    // detector centered at 0

  num = (int)*pptr++;
  if(num>1){

    inter = 2;
    for(int i=0; i<inter; i++){
      x[i] = *pptr++;
      y[i] = *pptr++;
      z[i] = *pptr++;
    }
  }
  /*
  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;

  // don't yet do isotope or other secondary production

  //
  // 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;
}


/*
  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.
*/
MuonPropagator::MuonPropagator(double newRmaxgen,
                               double newDepth,
                               double newEnergy,
                               double newCosTheta){

  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;
}

MuonPropagator::~MuonPropagator(){
  delete[] Eneutron;
  delete[] Rneutron;
  delete[] Tneutron;
  delete[] Eelectron;
  delete[] Relectron;
  delete[] Telectron;
  delete[] Egamma;
  delete[] Rgamma;
  delete[] Tgamma;
}

MuonPropagator::MuonPropagator(const MuonPropagator& Muon){

  Rmaxgen = Muon.Rmaxgen;

  Nmuon = Muon.Nmuon;
  energy = Muon.energy;
  cosg = Muon.cosg;
  phig = Muon.phig;

  Nneutron = Muon.Nneutron;
  Nelectron = Muon.Nelectron;
  Ngamma = Muon.Ngamma;

  double* pold;
  double* pnew;

  Eneutron = new double[Nneutron];
  pold = Muon.Eneutron+Nneutron;
  pnew = Eneutron+Nneutron;
  while(pnew>Eneutron)*--pnew=*--pold;

  Rneutron = new double[3*Nneutron];
  pold = Muon.Rneutron+3*Nneutron;
  pnew = Rneutron+3*Nneutron;
  while(pnew>Rneutron)*--pnew=*--pold;

  Tneutron = new double[Nneutron];
  pold = Muon.Tneutron+Nneutron;
  pnew = Tneutron+Nneutron;
  while(pnew>Tneutron)*--pnew=*--pold;

  Eelectron = new double[Nelectron];
  pold = Muon.Eelectron+Nelectron;
  pnew = Eelectron+Nelectron;
  while(pnew>Eelectron)*--pnew=*--pold;

  Relectron = new double[3*Nelectron];
  pold = Muon.Relectron+3*Nelectron;
  pnew = Relectron+3*Nelectron;
  while(pnew>Relectron)*--pnew=*--pold;

  Telectron = new double[Nelectron];
  pold = Muon.Telectron+Nelectron;
  pnew = Telectron+Nelectron;
  while(pnew>Telectron)*--pnew=*--pold;

  Egamma = new double[Ngamma];
  pold = Muon.Egamma+Ngamma;
  pnew = Egamma+Ngamma;
  while(pnew>Egamma)*--pnew=*--pold;

  Rgamma = new double[3*Ngamma];
  pold = Muon.Rgamma+3*Ngamma;
  pnew = Rgamma+3*Ngamma;
  while(pnew>Rgamma)*--pnew=*--pold;

  Tgamma = new double[Ngamma];
  pold = Muon.Tgamma+Ngamma;
  pnew = Tgamma+Ngamma;
  while(pnew>Tgamma)*--pnew=*--pold;

  Weight = Muon.Weight;
}    

MuonPropagator& MuonPropagator::operator=(const MuonPropagator& rhs){

  if(this == &rhs)return *this;

  Rmaxgen = rhs.Rmaxgen;

  Nmuon = rhs.Nmuon;
  energy = rhs.energy;
  cosg = rhs.cosg;
  phig = rhs.phig;

  Nneutron = rhs.Nneutron;
  Nelectron = rhs.Nelectron;
  Ngamma = rhs.Ngamma;

  double* pold;
  double* pnew;

  delete[] Eneutron;
  Eneutron = new double[Nneutron];
  pold = rhs.Eneutron+Nneutron;
  pnew = Eneutron+Nneutron;
  while(pnew>Eneutron)*--pnew=*--pold;

  delete[] Rneutron;
  Rneutron = new double[3*Nneutron];
  pold = rhs.Rneutron+3*Nneutron;
  pnew = Rneutron+3*Nneutron;
  while(pnew>Rneutron)*--pnew=*--pold;

  delete[] Tneutron;
  Tneutron = new double[Nneutron];
  pold = rhs.Tneutron+Nneutron;
  pnew = Tneutron+Nneutron;
  while(pnew>Tneutron)*--pnew=*--pold;

  delete[] Eelectron;
  Eelectron = new double[Nelectron];
  pold = rhs.Eelectron+Nelectron;
  pnew = Eelectron+Nelectron;
  while(pnew>Eelectron)*--pnew=*--pold;

  delete[] Relectron;
  Relectron = new double[3*Nelectron];
  pold = rhs.Relectron+3*Nelectron;
  pnew = Relectron+3*Nelectron;
  while(pnew>Relectron)*--pnew=*--pold;

  delete[] Telectron;
  Telectron = new double[Nelectron];
  pold = rhs.Telectron+Nelectron;
  pnew = Telectron+Nelectron;
  while(pnew>Telectron)*--pnew=*--pold;

  delete[] Egamma;
  Egamma = new double[Ngamma];
  pold = rhs.Egamma+Ngamma;
  pnew = Egamma+Ngamma;
  while(pnew>Egamma)*--pnew=*--pold;

  delete[] Rgamma;
  Rgamma = new double[3*Ngamma];
  pold = rhs.Rgamma+3*Ngamma;
  pnew = Rgamma+3*Ngamma;
  while(pnew>Rgamma)*--pnew=*--pold;

  delete[] Tgamma;
  Tgamma = new double[Ngamma];
  pold = rhs.Tgamma+Ngamma;
  pnew = Tgamma+Ngamma;
  while(pnew>Tgamma)*--pnew=*--pold;

  Weight = rhs.Weight;

  return *this;
}

double* MuonPropagator::Intersection(double x1, double y1, double z1,
                                     double x2, double y2, double z2,
                                     double x3, double y3, double z3, double r){

  // x1,y1,z1  P1 coordinates (point of line)
  // x2,y2,z2  P2 coordinates (point of line)
  // x3,y3,z3, r  P3 coordinates and radius (sphere)
  // x,y,z   intersection coordinates
  //
  // This function returns a pointer array which first index indicates
  // the number of intersection point, followed by coordinate pairs.

  double a, b, c, mu, i ;
  // array containing number of points and its coordinates
  double point[7];
  double* p = point;

  a = (x2-x1)*(x2-x1) + (y2-y1)*(y2-y1) + (z2-z1)*(z2-z1);
  b = 2*((x2-x1)*(x1-x3) + (y2-y1)*(y1-y3) + (z2-z1)*(z1-z3));
  c = (x3*x3 + y3*y3 + z3*z3 + x1*x1 + y1*y1 + z1*z1 - 
       2*(x3*x1 + y3*y1 + z3*z1) - r*r);
  i = b*b - 4*a*c;
  //cout << "intersection " << i << endl;

  if(i < 0.0){
    // no intersection
    p[0] = 0.0;
  }
  if(i ==0.0){
    // one intersection
    p[0] = 1.0;

    mu = -b/(2*a) ;
    p[1] = x1 + mu*(x2-x1);
    p[2] = y1 + mu*(y2-y1);
    p[3] = z1 + mu*(z2-z1);
  }
  if(i>0.0){
    // two intersections
    p[0] = 2.0;

    // first intersection
    mu = (-b + sqrt(b*b - 4*a*c))/(2*a);
    p[1] = x1 + mu*(x2-x1);
    p[2] = y1 + mu*(y2-y1);
    p[3] = z1 + mu*(z2-z1);

    // second intersection
    mu = (-b - sqrt(b*b - 4*a*c))/(2*a);
    p[4] = x1 + mu*(x2-x1);
    p[5] = y1 + mu*(y2-y1);
    p[6] = z1 + mu*(z2-z1);
  }
  return(p);
}

double* MuonPropagator::Rotate(char axis,double angle,
                               double x1, double y1, double z1){

  // axis      axis about which to rotate
  // angle     angle of rotation in radians
  // x1,y1,z1  cartesian coordinates of point to rotate

  double x2,y2,z2;
  double p[3];
  double* ptr = p;

  if(axis == 'X'){
    x2 = x1;
    y2 = cos(angle)*y1-sin(angle)*z1;
    z2 = sin(angle)*y1+cos(angle)*z1;
  }
  else if(axis == 'Y'){
    x2 = cos(angle)*x1+sin(angle)*z1;
    y2 = y1;
    z2 = -sin(angle)*x1+cos(angle)*z1;
  }
  else if(axis == 'Z'){
    x2 = cos(angle)*x1-sin(angle)*y1;
    y2 = sin(angle)*x1+cos(angle)*y1;
    z2 = z1;
  }
  else{
    cout << "ERROR: Unknown axis " << axis 
         << ".  Axis of rotation must be 'X', 'Y', or 'Z'." << endl;
    x2 = 0;
    y2 = 0;
    z2 = 0;
  }

  p[0] = x2;
  p[1] = y2;
  p[2] = z2;

  return(ptr);
}

double MuonPropagator::RadCorr(double t,double tup){

  if(t<=0 || tup<=0 || t<=tup) return 0;
  //
  double y = 1/(2+t);
  //
  return log(t*tup)
    -tup*tup*(t+2*tup-3*tup*tup*y/2
              -(tup-exp(3*log(tup))/3)*y*y
              -(tup*tup/2
                -t*exp(3*log(tup))/2
                +exp(4*log(tup))/4)*exp(3*log(y)));
}

float MuonPropagator::IsotopeDistance(const float& r){

  // from Hagner et al, Astropart Phys 14, 37 (2000)
  float phi = 1/166.0*exp(-0.0403*r*r)+1/1223.0*exp(-0.098*r);
  return phi;
}

---