//$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$//*
//             LOW ENERGY MUON SPIN RELAXATION, ROTATION, RADIATION   Geant4 SIMULATION
//  ID    : LEMuSRMuonDecayChannel.cc , v 1.2
//  AUTHOR: Taofiq PARAISO based on G4MuonDecayChannel $Id$
//  DATE  : 2004-07-13 11:15
// 
//  add muonium decay, PARAISO 07/04/2005                  
//
//$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$//
//
//  &           &&&&&&&&&&                            &&&&&&&         &&&&&&&&
//  &           &                                   &&      &&       &       &&
//  &           &                   &      &       &                &       &&
//  &           &&&&&&&            &      &         &&&&&&         &&&&&&&&  
//  &           &                 &      &&                &      &      &&
//  &           &                &&     &  &     &&      &&      &        &
//  &&&&&&&&&&  &&&&&&&&&&      &  &&&&&    &&    &&&&&&&       &        &&   
//                             &
//                            &
//                           &
//                         &  
//                          MUON DECAY CHANNEL.CC
//$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$//



#include "Randomize.hh"
#include "G4ios.hh"
#include "G4ParticleDefinition.hh"
#include "G4DecayProducts.hh"
#include "G4VDecayChannel.hh"
#include "LEMuSRMuonDecayChannel.hh"
#include "Randomize.hh"
#include "G4LorentzVector.hh"
#include "G4LorentzRotation.hh"
#include "G4Transform3D.hh"
#include "G4UnitsTable.hh"

#include "LEMuSRAtRestSpinRotation.hh"


LEMuSRMuonDecayChannel::LEMuSRMuonDecayChannel(const G4String& theParentName, 
				       G4double        theBR)
                   :G4VDecayChannel("LEMuSR Muon Decay",1)
{

  // set names for daughter particles
  if (theParentName == "mu+") {
    SetBR(theBR);
    SetParent("mu+");
    SetNumberOfDaughters(3);
    SetDaughter(0, "e+");
    SetDaughter(1, "nu_e");
    SetDaughter(2, "anti_nu_mu");

  } else if (theParentName == "mu-") {
    SetBR(theBR);
    SetParent("mu-");
    SetNumberOfDaughters(3);
    SetDaughter(0, "e-");
    SetDaughter(1, "anti_nu_e");
    SetDaughter(2, "nu_mu");
  }
  else if (theParentName == "Mu") {
    SetBR(theBR);
    SetParent("Mu");
    SetNumberOfDaughters(3);
    SetDaughter(0, "e+");
    SetDaughter(1, "nu_e");
    SetDaughter(2, "anti_nu_mu");
  }  
  else {

#ifdef G4VERBOSE 
    if (GetVerboseLevel()>0) {
      G4cout << "LEMuSRMuonDecayChannel:: constructor :";
      G4cout << " parent particle is not muon but ";
      G4cout << theParentName << G4endl;
    }
#endif
    
  }
  
  // get the random number engine
  theEngine = HepRandom::getTheEngine();

  
  theParentPolarization=G4ThreeVector(0,0,0);
  pointer = this; 



}

LEMuSRMuonDecayChannel*  LEMuSRMuonDecayChannel::pointer=0;
LEMuSRMuonDecayChannel*  LEMuSRMuonDecayChannel::GetInstance()
{
  return pointer;
} 



LEMuSRMuonDecayChannel::~LEMuSRMuonDecayChannel()
{

}


G4DecayProducts *LEMuSRMuonDecayChannel::DecayItPolarized(G4double mass,G4ThreeVector polar)
{
  SetParentPolarization(polar);

#ifdef G4VERBOSE 
  if (GetVerboseLevel()>2) {
    G4cout << "LEMuSRMuonDecayChannel:: theParentPolarization is" <<theParentPolarization  <<G4endl;
  }
#endif
  return DecayIt(mass);
}



G4DecayProducts *LEMuSRMuonDecayChannel::DecayIt(G4double) 
{
  
  if(theParentPolarization == G4ThreeVector(0,0,0))
  {
    // Generate random theParentPolarization direction

    G4double cost = 1. - 2.*G4UniformRand();
    G4double sint = std::sqrt((1.-cost)*(1.+cost));

    G4double phi = twopi*G4UniformRand();
    G4double sinp = std::sin(phi);
    G4double cosp = std::cos(phi);

    G4double px = sint*cosp;
    G4double py = sint*sinp;
    G4double pz = cost;

    theParentPolarization.setX(px);
    theParentPolarization.setY(py);
    theParentPolarization.setZ(pz);

  }
  //          assumes the pure V-A coupling
  //              gives incorrect energy spectrum for neutrinos
  
#ifdef G4VERBOSE 
  if (GetVerboseLevel()>2) {
    G4cout << "LEMuSRMuonDecayChannel::DecayIt ";
  }
#endif
  
  if (parent == 0) FillParent();  
  if (daughters == 0) FillDaughters();
  
  // parent mass
  G4double parentmass = parent->GetPDGMass();
  
  //daughters'mass
  G4double daughtermass[3]; 
  G4double sumofdaughtermass = 0.0;
  for (G4int index=0; index<3; index++){
    daughtermass[index] = daughters[index]->GetPDGMass();
    sumofdaughtermass += daughtermass[index];
  }
  
  //create parent G4DynamicParticle at rest
  G4ThreeVector dummy;
  G4DynamicParticle * parentparticle = new G4DynamicParticle( parent, dummy, 0.0);
  //create G4Decayproducts
  G4DecayProducts *products = new G4DecayProducts(*parentparticle);
  delete parentparticle;
  
  
  //                             Calculate Positron energy

  G4double daughtermomentum[3];
  G4double energy;
  G4double x; // maximal positron energy is 52.831MeV
  G4double r;
  
  do {
    do {
      r = G4UniformRand();
      x = G4UniformRand();
    } while (r > (3.0 - 2.0*x)*x*x);    
    
    
    energy = x*52.831*MeV;
  } while (energy <0.0|| energy > 52.831);
  
  
  
  //                                Anglular Distribution
  
  G4double energymax=52.831*MeV;
  G4double E=energy/energymax;
   
  G4double D = (2*E-1)/(3-2*E);
  
  // theta
  G4double K = G4UniformRand();
  costheta = 1/D*(-1.+ sqrt(1.-2*D*(2*K-1)+D*D));
  theta = acos(costheta);
  sintheta=sin(theta);

  // phi
  phi  = 2.0*M_PI*G4UniformRand()*rad;
  sinphi = sin(phi);
  cosphi = cos(phi);

  // rotation angles
 G4double px = sintheta*sinphi;
 G4double py = sintheta*cosphi;
 G4double pz = costheta;
                                                                                        
 G4ThreeVector direction0(px,py,pz);
                                                                                          
 direction0.rotateUz(theParentPolarization); 

#ifdef G4VERBOSE 
  if (GetVerboseLevel()>2) {
    G4cout << "LEMuSRMuonDecayChanel::DecayIt \n";
    G4cout <<"e+ momentum direction: " <<  direction0 <<G4endl;
    G4cout <<"energy: " <<energy <<G4endl; 
  }
#endif
  

  daughtermomentum[0] = sqrt(energy*energy + 2.0*energy* daughtermass[0]);
  
  G4DynamicParticle * daughterparticle 
         = new G4DynamicParticle( daughters[0], direction0*daughtermomentum[0]);
  daughterparticle->SetPolarization(theParentPolarization.x(),theParentPolarization.y(),theParentPolarization.z());
  products->PushProducts(daughterparticle);



  // daughter 1 ,2 (neutrinos)
  // create neutrinos in the C.M frame of two neutrinos
  G4double energy2 = parentmass*(1.0 - x/2.0); 
  G4double vmass   = sqrt((energy2-daughtermomentum[0])*(energy2+daughtermomentum[0]));
  G4double beta = -1.0*daughtermomentum[0]/energy2;
  G4double costhetan = 2.*G4UniformRand()-1.0;
  G4double sinthetan = sqrt((1.0-costhetan)*(1.0+costhetan));
  G4double phin  = 2.0*M_PI*G4UniformRand()*rad;
  G4double sinphin = sin(phin);
  G4double cosphin = cos(phin);

  G4ThreeVector direction1(sinthetan*cosphin,sinthetan*sinphin,costhetan);
  G4DynamicParticle * daughterparticle1 
         = new G4DynamicParticle( daughters[1], direction1*(vmass/2.));
  G4DynamicParticle * daughterparticle2
         = new G4DynamicParticle( daughters[2], direction1*(-1.0*vmass/2.));

  // boost to the muon rest frame
  G4LorentzVector p4;
  p4 = daughterparticle1->Get4Momentum();
  p4.boost( direction0.x()*beta, direction0.y()*beta, direction0.z()*beta);
  daughterparticle1->Set4Momentum(p4);
  p4 = daughterparticle2->Get4Momentum();
  p4.boost( direction0.x()*beta, direction0.y()*beta, direction0.z()*beta);
  daughterparticle2->Set4Momentum(p4);
  products->PushProducts(daughterparticle1);
  products->PushProducts(daughterparticle2);
  daughtermomentum[1] = daughterparticle1->GetTotalMomentum();
  daughtermomentum[2] = daughterparticle2->GetTotalMomentum();
    
 // output message
#ifdef G4VERBOSE 
    if (GetVerboseLevel()>2) {
      G4cout << "LEMuSRMuonDecayChanel::DecayIt ";
      G4cout << "  THETA aNGLE ::" <<theta <<G4endl;
      G4cout << "  PHI aNGLE ::" <<phi <<G4endl;
    
      G4cout << "  create decay products in rest frame: " <<G4endl;
      products->DumpInfo();
    }
#endif


  return products;
}