//____________________________________________________________________________
/*!

\class   PropaMuon

\brief   A class for propagate muons 

\author  Carla Distefano <distefano_c@lns.infn.it>
         LNS-INFN, Catania

\created December 10, 2015

\cpright  Copyright (c) 2015-2019, The KM3NeT Collaboration 
          For the full text of the license see $GSEAGEN/LICENSE
*/
//____________________________________________________________________________


#include <iostream>
#include <cmath>
#include <algorithm>
#include <stdlib.h>


#include <TMath.h>
#include <TVector3.h>


#include "PropaUtils.h"
#include "PropaMuon.h"
#include "PropaIon.h"
#include "PropaBrem.h"
#include "PropaPair.h"
#include "PropaNucl.h"

using namespace std;
using namespace propamuon;

//____________________________________________________________________________
PropaMuon::PropaMuon(map<int,double> MediaComposition, double MediaDensity){
  double EMax=1.E11;
  this->Initialize(MediaComposition, MediaDensity, EMax);
}
//____________________________________________________________________________
PropaMuon::PropaMuon(map<int,double> MediaComposition, double MediaDensity, double EMax){
  this->Initialize(MediaComposition, MediaDensity, EMax);
}
//____________________________________________________________________________
void PropaMuon::Initialize(map<int,double> MediaComposition, double MediaDensity, double EMax){

  fVcut=3.E-3;
  fEneCut=1.;
  
  fEneMax=EMax;

  fMediaDensity = MediaDensity; // MediaDensity is g/cm^3 
  fMediaComposition = MediaComposition;
  this->ComputeMeanFreePaths();
  fRnd = new TRandom();
  cout<<"Initialization of PropaMuon... done."<<endl;
}
//___________________________________________________________________________
PropaMuon::~PropaMuon(){  

  delete fELossIon;
  delete fELossRad;

  delete fLambdaBrem;
  delete fLambdaPair;
  delete fLambdaNucl;
  delete fLambdaTot;

  delete fProbBremTot;
  delete fProbPairTot;
  delete fProbNuclTot;  

  delete fRnd;
}
//___________________________________________________________________________
void PropaMuon::ComputeMeanFreePaths(void){

  double EMin=fEneCut;
  double EMax=fEneMax;
  double LogEStep=0.1;
  
  int NEmuon=(log10(EMax)-log10(EMin))/LogEStep+1;
  double logE,Energy;
  
  double Vmin,Vmax; 
  double l_brem,l_pair,l_nucl,l_tot; 
  double b_brem,b_pair,b_nucl,b_tot;

  fLambdaBrem = new TGraph; 
  fLambdaPair = new TGraph;   
  fLambdaNucl = new TGraph;
  fLambdaTot = new TGraph;  
  
  fELossIon = new TGraph;
  fELossRad = new TGraph;
  
  PropaIon gPropaIon;
  PropaBrem gPropaBrem(2);   
  PropaPair gPropaPair;   
  PropaNucl gPropaNucl;   
  
  for(int iLogE=0; iLogE<NEmuon; iLogE++){   
    logE=TMath::Log10(EMin)+LogEStep*iLogE;
    Energy=pow(10.,logE);    
    Vmin=0.;
    Vmax=fVcut;
    
    b_brem = fMediaDensity*gPropaBrem.BCoeff(Energy,fMediaComposition,Vmin,Vmax);
    b_pair = fMediaDensity*gPropaPair.BCoeff(Energy,fMediaComposition,Vmin,Vmax);
    b_nucl = fMediaDensity*gPropaNucl.BCoeff(Energy,fMediaComposition,Vmin,Vmax);    
    b_tot  = b_brem+b_pair+b_nucl;
    fELossRad->SetPoint(iLogE,Energy,b_tot*100.);	// in Gev/m

    Vmin=fVcut;
    Vmax=1.;
    l_brem = gPropaBrem.Lambda(Energy,fMediaComposition,Vmin,Vmax)/fMediaDensity/100.; // in m
    l_pair = gPropaPair.Lambda(Energy,fMediaComposition,Vmin,Vmax)/fMediaDensity/100.; // in m
    l_nucl = gPropaNucl.Lambda(Energy,fMediaComposition,Vmin,Vmax)/fMediaDensity/100.; // in m
    l_tot=1./(1./l_brem+1./l_pair+1./l_nucl);
    fLambdaBrem->SetPoint(iLogE,Energy,l_brem); // in m
    fLambdaPair->SetPoint(iLogE,Energy,l_pair); // in m
    fLambdaNucl->SetPoint(iLogE,Energy,l_nucl); // in m
    fLambdaTot->SetPoint (iLogE,Energy,l_tot); // in m
  }  

  Vmin=fVcut;
  Vmax=1.;

  double sigma;

  double LogVStep=0.05;
  double V;
 
  double binEmin=log10(EMin)-LogEStep;
  double binEmax=log10(EMax)+LogEStep;

  int nBinE=(int)(binEmax-binEmin)/LogEStep;

  int nBinP=100;
  double binPmin=0.;
  double binPmax=1.;
  double PStep=(binPmax-binPmin)/(double)nBinP;
  
  int iPoint =0;

  fProbBremTot = new TH2D("","",nBinE,binEmin,binEmax,nBinP,binPmin,binPmax);
  fProbPairTot = new TH2D("","",nBinE,binEmin,binEmax,nBinP,binPmin,binPmax);
  fProbNuclTot = new TH2D("","",nBinE,binEmin,binEmax,nBinP,binPmin,binPmax);  

  for(int iLogE=0; iLogE<=nBinE; iLogE++){ 

    logE=binEmin+LogEStep*iLogE+LogEStep;
    Energy=pow(10.,logE);      

    TGraph * graph = new TGraph;
    iPoint=0;
    
    double SigmaTotBrem=gPropaBrem.SigmaTot(Energy,fMediaComposition,Vmin,Vmax);
    for(double LogV=log10(Vmin); LogV<=log10(Vmax); LogV+=LogVStep){
      V=pow(10.,LogV);
      sigma=gPropaBrem.SigmaTot(Energy,fMediaComposition,V,Vmax)/SigmaTotBrem;
      if(sigma<0)sigma=0.;	
      graph->SetPoint(iPoint,sigma,LogV);
      iPoint++; 
    }
    for(int i=1; i<=nBinP; i++  ){
      double P=binPmin+PStep/2.+(i-1)*PStep;
      double LogV=graph->Eval(P);
      fProbBremTot->SetBinContent(iLogE,i,LogV);
    }
    delete graph;
//    fProbBremTot->Print();

    graph = new TGraph;
    iPoint=0;
    
    double SigmaTotPair=gPropaPair.SigmaTot(Energy,fMediaComposition,Vmin,Vmax);
    for(double LogV=log10(Vmin); LogV<=log10(Vmax); LogV+=LogVStep){
      V=pow(10.,LogV);
     sigma=gPropaPair.SigmaTot(Energy,fMediaComposition,V,Vmax)/SigmaTotPair;
     if(sigma<0)sigma=0.;	
       graph->SetPoint(iPoint,sigma,LogV);
       iPoint++; 
     }
      for(int i=1; i<=nBinP; i++  ){
      double P=binPmin+PStep/2.+(i-1)*PStep;
      double LogV=graph->Eval(P);
      fProbPairTot->SetBinContent(iLogE,i,LogV);
    }
    delete graph;
      
    graph = new TGraph;
    iPoint=0;
    
    double SigmaTotNucl=gPropaNucl.SigmaTot(Energy,fMediaComposition,Vmin,Vmax);
    for(double LogV=log10(Vmin); LogV<=log10(Vmax); LogV+=LogVStep){
      V=pow(10.,LogV);
     sigma=gPropaNucl.SigmaTot(Energy,fMediaComposition,V,Vmax)/SigmaTotNucl;
     if(sigma<0)sigma=0.;	
       graph->SetPoint(iPoint,sigma,LogV);
       iPoint++; 
     }
     for(int i=1; i<=nBinP; i++  ){
       double P=binPmin+PStep/2.+(i-1)*PStep;
       double LogV=graph->Eval(P);
       fProbNuclTot->SetBinContent(iLogE,i,LogV);
     }
     delete graph;
  } 

  EMin=kMuonMass;
  NEmuon=(log10(EMax)-log10(EMin))/LogEStep+1;
  double a_ion;
  for(int iLogE=0; iLogE<NEmuon; iLogE++){   
    logE=TMath::Log10(EMin)+LogEStep*iLogE;
    Energy=pow(10.,logE); 
    a_ion=fMediaDensity*gPropaIon.ACoeff(Energy,fMediaComposition)*100.;  //  in GeV/m    
    fELossIon->SetPoint(iLogE,Energy,a_ion);
  }      
  
  // Radiation Lenght
  
  fRadLen=0.;  
  map<int,double>::iterator iter = fMediaComposition.begin();
  for ( ; iter != fMediaComposition.end(); ++iter) {    
    double Percentage=iter->second; 
    double Z=propamuon::PdgToZ(iter->first);
    double A=propamuon::PdgToA(iter->first);
    fRadLen += Percentage*fMediaDensity*(4.*kAem*Z*(Z+1)*kRe*kRe*log(183.*pow(Z,-1./3.))/A);    
  }  
  fRadLen=1./fRadLen/100.; // in m
  
  return;
}  

//___________________________________________________________________________
bool PropaMuon::Propagate(PropaTrack TrackIn, double Depth){
  
  PropaTrack TrackOut=TrackIn;
  PropaTrack TrackInt=TrackOut; //To store intermediate states
  double rnd;

  fPath=0.;
  fRange=0.;
  fPathTot=0.;
  
  fDepth=Depth;
  fDepthToDo=Depth;
  fDepthDone=0.;

  if(fDepth<=0)return false;  

  if(TrackIn.E>fEneMax){
    cout<<"WARNING PropaMuon: Muon energy greater than maximum value, muon not propagated..."<<endl;
    return false;
  }
//cout<<"parto propamuon"<<fDepth<<" "<<TrackIn.E<<" "<<TrackIn.V3<<endl;
  while(TrackOut.E>kMuonMass){
    TrackInt=TrackOut; //Last final state is now TrackInt
    if(TrackOut.E<=fEneCut){
    // only Ion up to stopping
//    cout<<"case 1 "<<fPath<<" "<<fDepthToDo<<endl;
    fPath=fDepthToDo;
    
    TrackOut = this->ContinuousLosses(TrackOut);
    break;
    }    
    else{
      // simulation of the interaction point --> we extract random the path 
      double TotalLambda=fLambdaTot->Eval(TrackOut.E);
      rnd = fRnd->Rndm();
      fPath = -TotalLambda*TMath::Log(rnd);  
      
      if(fPath>fDepthToDo){
        // interaction point deeper than depth, continuous losses up to Depth and break
//	    cout<<"case 2 "<<fPath<<" "<<fDepthToDo<<endl;

	fPath=fDepthToDo;
	TrackOut = this->ContinuousLosses(TrackOut);
	break;
      }
      else {
        // continuous energy losses up to the next interaction point	
//	    cout<<"case 3 "<<fPath<<" "<<fDepthToDo<<endl;

        TrackOut = this->ContinuousLosses(TrackOut);
        // hard rad process
      	if(TrackOut.E>fEneCut)TrackOut=StochasticLosses(TrackOut);

	this->ComputeRange(TrackInt,TrackOut); //Now calculates distance travelled between last state and this one
	fDepthToDo-=fRange;			// we have to improve calculation of fDepthToDo
	if(fDepthToDo<0.01)break;		
      
      }    
    
    }
  }



/*
  while(TrackOut.E>kMuonMass){
  
    if(TrackOut.E<=fEneCut){
      // only Ion up to stopping 
      double CosTheta=TrackIn.D1*TrackOut.D1+TrackIn.D2*TrackOut.D2+TrackIn.D3*TrackOut.D3;
      if(CosTheta==0)fPath=0;
      else fPath = (fDepth-fDepthDone)/CosTheta;	/// dobbiamo dirgli per quanto propagare....
      TrackOut = this->ContinuousLosses(TrackOut);
      this->ComputeRange(TrackIn,TrackOut);
//      this->ComputeDepthDone(TrackIn,TrackOut);
      break;
    }
    else {
      // simulation of the interaction point --> we extract random the path 
      double TotalLambda=fLambdaTot->Eval(TrackOut.E);
      rnd = fRnd->Rndm();
      fPath = -TotalLambda*TMath::Log(rnd);   
//      double CosTheta=TrackIn.D1*TrackOut.D1+TrackIn.D2*TrackOut.D2+TrackIn.D3*TrackOut.D3;
//      double DepthStep=fPath*CosTheta;
      
      
      if(fPath>fDepthToDo){
        // interaction point deeper than depth, continuous losses up to Depth and break
	fPath=fDepthToDo;
        TrackOut = this->ContinuousLosses(TrackOut);
        this->ComputeRange(TrackIn,TrackOut);
//	this->ComputeDepthDone(TrackIn,TrackOut);	
        break;
      }
      else{
        // continuous energy losses up to the next interaction point	
        TrackOut = this->ContinuousLosses(TrackOut);
        this->ComputeRange(TrackIn,TrackOut);
//	this->ComputeDepthDone(TrackIn,TrackOut);
	
	if(TrackOut.E>fEneCut){
          // hard rad process
	  TrackOut=StochasticLosses(TrackOut);
	  this->ComputeRange(TrackIn,TrackOut);
//	  this->ComputeDepthDone(TrackIn,TrackOut);
	}
	
	this->ComputeDepthToDo();
	
      }
    }
    
  }
*/
//cout<<"propa "<<inter<<endl;
  this->ComputeRange(TrackIn,TrackOut);
//  this->ComputeDepthDone(TrackIn,TrackOut);

/*
  if(TrackOut.E>kMuonMass){
  cout<<endl<<"Summary "<<endl;

  cout<<"Depth"<<endl;
  cout<<fDepth<<endl;
  
  cout<<"DepthDone"<<endl;
  cout<<fDepthDone<<endl;
  
  cout<<"Range"<<endl;
  cout<<fRange<<" "<<fPathTot<<endl;

  cout<<"Energy"<<endl;
  cout<<TrackIn.E<<" "<<TrackOut.E<<endl;

  cout<<"Position"<<endl;
  cout<<TrackIn.V1<<" "<<TrackOut.V1<<endl;
  cout<<TrackIn.V2<<" "<<TrackOut.V2<<endl;
  cout<<TrackIn.V3<<" "<<TrackOut.V3<<endl;

  cout<<"Direction"<<endl;
  cout<<TrackIn.D1<<" "<<TrackOut.D1<<endl;
  cout<<TrackIn.D2<<" "<<TrackOut.D2<<endl;
  cout<<TrackIn.D3<<" "<<TrackOut.D3<<endl;

  cout<<"Time"<<endl;
  cout<<TrackIn.T<<" "<<TrackOut.T<<endl;
  }
*/

  TrackOut.Pdg=TrackIn.Pdg;
  TrackOut.Id=TrackIn.Id;
  
  fTrackOut=TrackOut;
  if(TrackOut.E>kMuonMass)return true;
  else return false;

}

//___________________________________________________________________________
void PropaMuon::ComputeRange(PropaTrack TrackIn, PropaTrack TrackOut){
  fRange=TMath::Power((TrackIn.V1-TrackOut.V1),2)+TMath::Power((TrackIn.V2-TrackOut.V2),2)+TMath::Power((TrackIn.V3-TrackOut.V3),2);
  fRange=TMath::Sqrt(fRange);
  return;
}
//___________________________________________________________________________
void PropaMuon::ComputeDepthDone(PropaTrack TrackIn, PropaTrack TrackOut){  
  fDepthDone=TrackIn.D1*(TrackOut.V1-TrackIn.V1)+TrackIn.D2*(TrackOut.V2-TrackIn.V2)+TrackIn.D3*(TrackOut.V3-TrackIn.V3);
  fDepthDone=TMath::Abs(fDepthDone);  
  return;
}
//___________________________________________________________________________
double PropaMuon::PathOnlyContinuous(double Energy){

  double a=fELossIon->Eval(Energy);
  double path;
  
  if(Energy<=fEneCut){
    path=(Energy-kMuonMass)/a;
  }
  else {
    double b=fELossRad->Eval(Energy);
    path=1./b*TMath::Log((a+b*Energy)/(a+b*kMuonMass));  
  }

  return path;
}
//___________________________________________________________________________
PropaTrack PropaMuon::ContinuousLosses(PropaTrack TrackIn){

  PropaTrack TrackOut;

  double EnergyF = this->ContinuousEnergyLosses(TrackIn.E);
  
  TrackOut=this->MultipleScattering(TrackIn);

  // check if propagated for a depth greater than the requested one because of multiple scattering
  this->ComputeDepthDone(TrackIn, TrackOut);
  if(fDepthDone>fDepthToDo){
    TrackOut.V1 -= TrackOut.D1*(fDepthDone-fDepthToDo);
    TrackOut.V2 -= TrackOut.D2*(fDepthDone-fDepthToDo);
    TrackOut.V3 -= TrackOut.D3*(fDepthDone-fDepthToDo);
    fPath-=fDepthDone-fDepthToDo;
    EnergyF = this->ContinuousEnergyLosses(TrackIn.E);
  }

  TrackOut.E=EnergyF;
  
  
  fPathTot+=fPath;

  // computing propagation time
  double GammaMean=(TrackIn.E+TrackOut.E)/2./kMuonMass;
  double Beta=TMath::Sqrt(GammaMean*GammaMean-1.)/GammaMean;
  double LightSpeed=0.299792458;// in m/ns 
  TrackOut.T=fPath/(Beta*LightSpeed)+TrackIn.T; // in ns

  return TrackOut; 
}
//___________________________________________________________________________
double PropaMuon::ContinuousEnergyLosses(double Energy){
  // We assume a and b constant and eveluated at the initial energy
  
  double EnergyF;
  double PathStop;
  
  double a=fELossIon->Eval(Energy);
  double b=0.;
  
  PathStop=this->PathOnlyContinuous(Energy);
  
  if(fPath>=PathStop){
    fPath=PathStop;
    EnergyF=kMuonMass;
  }
  else {
    if(Energy<=fEneCut){
      EnergyF=Energy-a*fPath;    
    }
    else{
      b=fELossRad->Eval(Energy);
      EnergyF=1./b*( (a+b*Energy)*TMath::Exp(-b*fPath)-a );
    }  
  }  


  return EnergyF;
}
//__________________________________________________________________________  
PropaTrack PropaMuon::MultipleScattering(PropaTrack track){

// Nuclear Instruments and methods in Physics Research B58 (7991) 6-10
// 
 
  PropaTrack trackf;

  double Rnd1, Rnd2;
  double sx, tx, sy, ty, sz, tz, ax, ay, az;

  double Theta0 = fPath/fRadLen;  
  double Momentum2= track.E*track.E-kMuonMass2;

  double Sqrt2=TMath::Sqrt(2);
  double Sqrt3=TMath::Sqrt(3);

  double Beta = 1./sqrt(1+kMuonMass2/Momentum2);
  Theta0 = 13.6/(sqrt(Momentum2)*Beta)*sqrt(Theta0)*(1.+0.088*log10(Theta0));

  Rnd1 = Sqrt2*Theta0*TMath::ErfInverse(2.*(fRnd->Rndm()-0.5));
  Rnd2 = Sqrt2*Theta0*TMath::ErfInverse(2.*(fRnd->Rndm()-0.5));

  sx = (Rnd1/Sqrt3+Rnd2)/2;
  tx = Rnd2;
  
  Rnd1 = Sqrt2*Theta0*TMath::ErfInverse(2.*(fRnd->Rndm()-0.5));
  Rnd2 = Sqrt2*Theta0*TMath::ErfInverse(2.*(fRnd->Rndm()-0.5));

  sy = (Rnd1/Sqrt3+Rnd2)/2;
  ty = Rnd2;

  sz = sqrt(max(1.-(sx*sx+sy*sy), 0.));
  tz = sqrt(max(1.-(tx*tx+ty*ty), 0.));

  double theta=acos(track.D3);
  double phi=atan2(track.D2,track.D1);
  
  double sinth, costh,sinph,cosph;

  sinth = sin(theta);  
  costh = cos(theta);
  
  sinph = sin(phi);
  cosph = cos(phi);

  ax = sinth*cosph*sz+costh*cosph*sx-sinph*sy;
  ay = sinth*sinph*sz+costh*sinph*sx+cosph*sy;
  az = costh*sz-sinth*sx;

  trackf.V1 = track.V1 + ax*fPath;		
  trackf.V2 = track.V2 + ay*fPath;
  trackf.V3 = track.V3 + az*fPath;

  ax = sinth*cosph*tz+costh*cosph*tx-sinph*ty;
  ay = sinth*sinph*tz+costh*sinph*tx+cosph*ty;
  az = costh*tz-sinth*tx;

  costh = az;
  sinth = sqrt(max(1-costh*costh, 0.));

  if(sinth!=0){
    sinph   =   ay/sinth;
    cosph   =   ax/sinth;
  }

  trackf.D1=sinth*cosph;
  trackf.D2=sinth*sinph;
  trackf.D3=costh;

  return trackf;
}

//___________________________________________________________________________
PropaTrack PropaMuon::StochasticLosses(PropaTrack TrackIn){

  PropaTrack TrackOut=TrackIn;
  ScatterOut scatt;
  
  // select the stochastic radiative process
  string ProcessName = this->SelectProcess(TrackOut.E);
  // compute V

  if(ProcessName.compare("Ion")!=0){

    StochasOut stoch=this->ComputeV(ProcessName,TrackOut.E);
    TrackOut.E-=stoch.V*TrackOut.E;

    if(ProcessName.compare("Nucl")==0){
      scatt=this->NuclScattering(TrackIn.E,stoch);
    }
    else{
      scatt=this->BremPairScattering(TrackIn.E,stoch,ProcessName);  
    }


    TrackOut.V1=0.;
    TrackOut.V2=0.;
    TrackOut.V3=0.;

    TrackOut.D1=TMath::Sin(scatt.Theta)*TMath::Cos(scatt.Phi);
    TrackOut.D2=TMath::Sin(scatt.Theta)*TMath::Sin(scatt.Phi);
    TrackOut.D3=TMath::Cos(scatt.Theta);

    TrackOut=this->RotoTranslation(TrackIn, TrackOut);

  }

  return TrackOut; 
}
//___________________________________________________________________________
string PropaMuon::SelectProcess(double Energy){
  string ProcessName;

  double BremLambda=fLambdaBrem->Eval(Energy);
  double PairLambda=fLambdaPair->Eval(Energy) ;
  double NuclLambda=fLambdaNucl->Eval(Energy);
  
  double TotalLambda=fLambdaTot->Eval(Energy);

  double BremProb=TotalLambda/BremLambda;
  double PairProb=TotalLambda/PairLambda;
  double NuclProb=TotalLambda/NuclLambda;

  double TotalProb=BremProb+PairProb+NuclProb;   


  BremProb/=TotalProb;
  PairProb/=TotalProb;
  NuclProb/=TotalProb;

  double rnd = fRnd->Rndm();
  

  if(rnd<BremProb)ProcessName="Brem";
  else if(rnd<BremProb+PairProb)ProcessName="Pair";
  else ProcessName="Nucl";

//  cout<<ProcessName<<endl; 
 
  return ProcessName;
}
//___________________________________________________________________________
StochasOut PropaMuon::ComputeV(string ProcessName, double Energy){
  
  StochasOut stoch;
  stoch.V=0.;
  stoch.PdgCode=0;
  
  if(ProcessName.compare("Ion")==0)return stoch;
    
  // tagert
  double rnd=fRnd->Rndm();  
  double PercTot=0.;
  map<int,double>::iterator iter = fMediaComposition.begin();
  for ( ; iter != fMediaComposition.end(); ++iter) {    
    PercTot += iter->second;
    if(rnd<PercTot){
      stoch.PdgCode=iter->first;
      break;
    }      
  }  
  
  rnd=fRnd->Rndm();
  if(ProcessName.compare("Brem")==0)
    stoch.V=fProbBremTot->Interpolate(TMath::Log10(Energy),rnd); 
  else if(ProcessName.compare("Pair")==0)
    stoch.V=fProbPairTot->Interpolate(TMath::Log10(Energy),rnd); 
  else if(ProcessName.compare("Nucl")==0)
    stoch.V=fProbNuclTot->Interpolate(TMath::Log10(Energy),rnd); 

  stoch.V=TMath::Power(10.,stoch.V);  

  return stoch;
}
//___________________________________________________________________________
ScatterOut PropaMuon::BremPairScattering(double Energy, StochasOut stoch, string ProcessName){

// From A. Van Ginneken, NIM A251 (1986) 21

  ScatterOut  scatt;
  
  scatt.Theta=0.;
  scatt.Phi=0.;
  
  double Theta2Mean=0.;

  if(ProcessName.compare("Brem")==0){
    double k1=0.092*TMath::Power(Energy,-1./3.);
    double k2=0.052/Energy*TMath::Power(propamuon::PdgToZ(stoch.PdgCode),-1./4.);
    double k3=0.22*TMath::Power(Energy,-0.92);
    double k4=0.26*TMath::Power(Energy,-0.91);
    double n=0.81*TMath::Power(Energy,1./2.)/(TMath::Power(Energy,1./2.)+1.8);
    double Theta1=TMath::Max(TMath::Min(k1*TMath::Sqrt(stoch.V),k2),k3*stoch.V);
    double Theta2=k4*TMath::Power(stoch.V,1.+n)*TMath::Power(1-stoch.V,-n); 
    double k5=k4*TMath::Power(0.5,1.+n)*TMath::Power(0.5,-n)/TMath::Power(0.5,-0.5);
    double Theta3=k5*TMath::Power(1-stoch.V,-1./2.);
    
    if(stoch.V<=0.5)Theta2Mean=Theta1;
    else {
      if(Theta2<0.2)Theta2Mean=Theta2;
      else Theta2Mean=Theta3;
    }
  
  }
  else if(ProcessName.compare("Pair")==0){
    Theta2Mean=2.3+TMath::Log(Energy)/(Energy*(1-stoch.V))*TMath::Power((stoch.V-2.*kElectronMass)/(Energy*stoch.V),2);
    Theta2Mean*= TMath::Min(8.9E-4* TMath::Power(stoch.V,1./4.)*(1.+1.5E-5*Energy)+0.032*stoch.V/(stoch.V+1),0.1);    
  }
  
  Theta2Mean=TMath::Power(Theta2Mean,2);
   
  double rnd;
  while(1){
    rnd=fRnd->Rndm();
    scatt.Theta=TMath::Sqrt(-Theta2Mean*TMath::Log(rnd));
    if(scatt.Theta<=kPi)break;
  }
  
  rnd=fRnd->Rndm();
  scatt.Phi=2.*kPi*rnd;

  return scatt;
}
//___________________________________________________________________________
ScatterOut PropaMuon::NuclScattering(double Energy, StochasOut stoch){
  ScatterOut  scatt;

  scatt.Theta=0.;
  scatt.Phi=0.;
  
  double m02=0.4;   // Gev^2

  double DEnergy=(1.-stoch.V)*Energy;
  double DEnergy2=TMath::Power(DEnergy,2);

  double tmin=TMath::Power(kMuonMass*(1-stoch.V),2)/stoch.V;
  double tmax=2.*kProtonMass*DEnergy;

  double t1=TMath::Min(DEnergy2,m02);

  double rnd=fRnd->Rndm();
  double trnd=tmax*t1/((tmax+t1)*TMath::Power((tmax*(tmin+t1))/(tmin*(tmax+t1)),rnd)-tmax);
  
  scatt.Theta=(trnd-tmin)/(4.*(stoch.V*Energy*Energy-kMuonMass*kMuonMass)-2.*tmin);  

  if(scatt.Theta<0)scatt.Theta*=-1.;

  scatt.Theta=TMath::Sqrt(scatt.Theta);
  scatt.Theta=2.*TMath::ASin(scatt.Theta);

  rnd=fRnd->Rndm();
  scatt.Phi=2.*kPi*rnd;

  return scatt;
}   
//___________________________________________________________________________
PropaTrack PropaMuon::RotoTranslation(PropaTrack TrackIn, PropaTrack TrackOut){

  double Theta= TMath::ACos(TrackIn.D3);
  double Phi=   TMath::ATan2(TrackIn.D2,TrackIn.D1);

  double X0=TrackIn.V1;
  double Y0=TrackIn.V2;
  double Z0=TrackIn.V3;

  TrackIn=TrackOut;

  double l1=TMath::Cos(Theta)*TMath::Cos(Phi)*TMath::Cos(Phi)+TMath::Sin(Phi)*TMath::Sin(Phi);
  double m1=TMath::Cos(Theta)*TMath::Cos(Phi)*TMath::Sin(Phi)-TMath::Cos(Phi)*TMath::Sin(Phi);
  double n1=-TMath::Sin(Theta)*TMath::Cos(Phi);
  
  double l2=TMath::Cos(Theta)*TMath::Cos(Phi)*TMath::Sin(Phi)-TMath::Cos(Phi)*TMath::Sin(Phi);
  double m2=TMath::Cos(Theta)*TMath::Sin(Phi)*TMath::Sin(Phi)+TMath::Cos(Phi)*TMath::Cos(Phi);
  double n2=-TMath::Sin(Theta)*TMath::Sin(Phi);
  
  double l3=TMath::Sin(Theta)*TMath::Cos(Phi);
  double m3=TMath::Sin(Theta)*TMath::Sin(Phi);
  double n3=TMath::Cos(Theta);


  // we have to rotate vertex and direction

  TrackOut.V1=l1*TrackIn.V1+l2*TrackIn.V2+l3*TrackIn.V3+X0;
  TrackOut.V2=m1*TrackIn.V1+m2*TrackIn.V2+m3*TrackIn.V3+Y0;
  TrackOut.V3=n1*TrackIn.V1+n2*TrackIn.V2+n3*TrackIn.V3+Z0;
/*
  TrackOut.V1=TrackIn.V1;
  TrackOut.V2=TrackIn.V2;
  TrackOut.V3=TrackIn.V3;
*/

  TrackOut.D1=l1*TrackIn.D1+l2*TrackIn.D2+l3*TrackIn.D3;
  TrackOut.D2=m1*TrackIn.D1+m2*TrackIn.D2+m3*TrackIn.D3;
  TrackOut.D3=n1*TrackIn.D1+n2*TrackIn.D2+n3*TrackIn.D3;
  
  return TrackOut;
} 
//___________________________________________________________________________
PropaTrack PropaMuon::GetPropaTrack(void){return fTrackOut;};
//___________________________________________________________________________
double PropaMuon::GetRange(void){return fRange;};
//___________________________________________________________________________
double PropaMuon::GetTotPath(void){return fPathTot;};