---

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 "ReactorConstants.hh"
#include "ReactorFortran.hh"
#include "MuonPropagator.hh"
#include "MyHeSource.hh"
#include "MyLiSource.hh"

MuonPropagator::MuonPropagator(int newRmaxgen, double newDepth){

  Rmaxgen = newRmaxgen;
  Depth = newDepth;

  //cout << "MuonPropagator::MuonPropagator(" << Rmaxgen << "," 
  //     << Depth << ")" << endl;

  static double re  = k.Relectron;
  static double pi  = k.pi;
  static double me  = k.Melectron;
  static double mmu = k.Mmuon;
  static double Na  = k.Navogadro;
  static double c   = k.c;
  static double density = k.density;
  //
  // For electron density
  //
  static double nel[2];
  static double A[2] = {k.A0,k.A1};
  static double Z[2] = {k.Z0,k.Z1};
  static double f[2] = {k.f0,k.f1};
  static double I[2] = {k.I0,k.I1};
  //
  static double ncar;
  //
  // create the energy and cosine of zenith angle (theta) arrays
  bool static first = true;
  double static Energy[Mdim],CosTheta[Mdim];
  double flux;
  int static muonstop;
  float static fspace[200];
  if(first){

    if(Depth == 500){
      // read in the muon energy in GeV and angle
      fstream fin("data/Totb500a.dat", ios::in);
      //cout << "setup file: data/Totb500a.dat" << endl;

      double dummy;
      for(int i=0;i<Mdim;i++){
        fin >> dummy >> Energy[i] >> dummy >> dummy >> dummy >> CosTheta[i];
        //cout << i << "  " << Energy[i] << "  " << CosTheta[i] << endl;
      }

      flux = 0.145;  // per m**2 per sec
    }
    else{
      cout << "  ERROR: no muon data 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.;
    float q = 0.;
    funlxp_(IsotopeDistance,fspace,xlo,xhi);

    Weight = 0.;

    // to count stopped muons
    muonstop = 0;

    first = false;
  }

  /* get the muon energy and cosine of the zenith angle
     the spectra in Totb500a.dat are random so start 
     with the first entry and cycle through the file
  */
  static int index = 0;
  if(index < Mdim){
    //
    // muon energy in MeV
    energy = Energy[index]*1000.;
    cosg = CosTheta[index];
    //cout << " muon energy " << energy << " MeV, cos(theta) " << cosg << endl;
    index++;
  }
  else{
    cout << "  WARNING: end of muon input file reached at index " << index 
         << ".  Starting over." << endl;
    index = 0;
  }

  /* 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;
  double r[6],t[3];
  double x[2],y[2],z[2];
  double delta = Rmaxgen*20.;
  bool passed = false;
  //
  while(!passed){

    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*k.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{
      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
  //
  int len=2;
  float rw[len];
  ranlux_(rw,len);

  double newphig = 2*k.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;
  len = 1;
  funlux_(fspace,lateraldistance,len);
  //
  // 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];

  /* 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]);
  //
  // second rotation
  //
  angle = pi - phig;
  double pt[3];
  double* ptptr = pt;
  ptptr = Rotate('Z',angle,qt[0],qt[1],qt[2]);
  //
  // 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];
  //
  // 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 x , y , z;
  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];

  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(p);
}
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;
}

---