//$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
//             LOW ENERGY MUON SPIN RELAXATION, ROTATION, RADIATION 
//
//  ID    : LEMuSRMuonPhysics.cc , v 1.3
//  AUTHOR: Taofiq PARAISO
//  DATE  : 2006-01-19 16:15
//$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
//
//  &           &&&&&&&&&&                            &&&&&&&         &&&&&&&&
//  &           &                                   &&      &&       &       &&
//  &           &                   &      &       &                &       &&
//  &           &&&&&&&            &      &         &&&&&&         &&&&&&&&  
//  &           &                 &      &&                &      &      &&
//  &           &                &&     &  &     &&      &&      &        &
//  &&&&&&&&&&  &&&&&&&&&&      &  &&&&&    &&    &&&&&&&       &        &&   
//                             &
//                            &
//                           &
//                         &  
//                               MuonPhysics   
//$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$

#include "LEMuSRMuonPhysics.hh"
 
#include "globals.hh"
#include "G4ios.hh"
#include <iomanip>   


LEMuSRMuonPhysics::LEMuSRMuonPhysics(const G4String& name)
                   :  G4VPhysicsConstructor(name)
{
}

LEMuSRMuonPhysics::~LEMuSRMuonPhysics()
{
} 

#include "G4ParticleDefinition.hh"
#include "G4ParticleTable.hh"

#include "G4MuonPlus.hh"
#include "G4Muonium.hh"
#include "G4MuonMinus.hh"
#include "G4TauMinus.hh"
#include "G4TauPlus.hh"
#include "G4NeutrinoTau.hh"
#include "G4AntiNeutrinoTau.hh"
#include "G4NeutrinoMu.hh"
#include "G4AntiNeutrinoMu.hh"
#include "LEMuSRMuonDecayChannel.hh"
#include "G4UserSpecialCuts.hh"
#include "G4DecayTable.hh"
#include "G4MuonDecayChannelWithSpin.hh"




void LEMuSRMuonPhysics::ConstructParticle()
{
  //! Muon
  G4MuonPlus::MuonPlusDefinition();
  G4Muonium::MuoniumDefinition();
  G4MuonMinus::MuonMinusDefinition();
  G4NeutrinoMu::NeutrinoMuDefinition();
  G4AntiNeutrinoMu::AntiNeutrinoMuDefinition();
 
  //! Tau
  G4TauMinus::TauMinusDefinition();
  G4TauPlus::TauPlusDefinition();
  G4NeutrinoTau::NeutrinoTauDefinition();
  G4AntiNeutrinoTau::AntiNeutrinoTauDefinition();

  
  //! The decay process
  /*!
   * Although the asymmetric decay process is a feature of new Geant4 versions, 
   * it is recommended to use the original \lemu implementation. 
*/
#ifdef ASYM_USE_LEMU
  G4DecayTable* MuonPlusDecayTable = new G4DecayTable();
  MuonPlusDecayTable -> Insert(new LEMuSRMuonDecayChannel("mu+",1.00));
  G4MuonPlus::MuonPlusDefinition() -> SetDecayTable(MuonPlusDecayTable);

  /*!
    Half of the muonium have a completely isotropic decay. 
   */
  G4DecayTable* MuoniumDecayTable = new G4DecayTable();
  MuoniumDecayTable -> Insert(new LEMuSRMuonDecayChannel("Mu",0.5));
  MuoniumDecayTable -> Insert(new G4MuonDecayChannel("Mu",0.5)); 
  G4Muonium::MuoniumDefinition() -> SetDecayTable(MuoniumDecayTable);

  G4DecayTable* MuonMinusDecayTable = new G4DecayTable();
  MuonMinusDecayTable -> Insert(new LEMuSRMuonDecayChannel("mu-",1.00));
  G4MuonMinus::MuonMinusDefinition() -> SetDecayTable(MuonMinusDecayTable);
  G4cout<<"ASYM_USE_LEMU\n";
#else
  G4DecayTable* MuonPlusDecayTable = new G4DecayTable();
  MuonPlusDecayTable -> Insert(new G4MuonDecayChannel("mu+",1.00));
  G4MuonPlus::MuonPlusDefinition() -> SetDecayTable(MuonPlusDecayTable);
  
  G4DecayTable* MuoniumDecayTable = new G4DecayTable();
  MuoniumDecayTable -> Insert(new G4MuonDecayChannel("Mu",1.));
  G4Muonium::MuoniumDefinition() -> SetDecayTable(MuoniumDecayTable);
  
  G4DecayTable* MuonMinusDecayTable = new G4DecayTable();
  MuonMinusDecayTable -> Insert(new G4MuonDecayChannelWithSpin("mu-",1.00));
  G4MuonMinus::MuonMinusDefinition() -> SetDecayTable(MuonMinusDecayTable);
  G4cout<<"ASYM_USE_G4\n";
#endif
   

}


#include "G4ProcessManager.hh"

void LEMuSRMuonPhysics::ConstructProcess()
{
  G4ProcessManager * pManager = 0;

  // Muon Plus Physics
  pManager = G4MuonPlus::MuonPlus()->GetProcessManager();
     
//! Syntax:   pManager->AddProcess(&fprocess, atrest, prestep, poststep). -1 means that the process is not executed under correspondant condition.
  pManager->AddProcess(&fMuAtRestSpinRotation, 1, -1, -1); 
  pManager->AddProcess(&fMuFormation,-1,  -1, 1);
  pManager->AddProcess(&fMuPlusMultipleScattering,-1,  1, -1);// may crash when enabled as post step process
  pManager->AddProcess(&fMuPlusIonisation,        -1,  2, 2);
  pManager->AddProcess(&fMuPlusBremsstrahlung,    -1,  3, 3);
  pManager->AddProcess(&fMuPlusPairProduction,    -1,  4, 4);
  
  pManager->AddProcess(new G4UserSpecialCuts(),-1,-1,6); 
  //  pManager ->SetProcessOrderingToLast(&fDepolarization, idxAtRest);
  //  pManager ->SetProcessOrdering(&fDepolarization, idxPostStep);

  //! the last process: decay
  pManager->AddProcess(&fDecayProcess);
  pManager ->SetProcessOrderingToLast(&fDecayProcess, idxAtRest);
  pManager ->SetProcessOrdering(&fDecayProcess, idxPostStep);


  //! Muonium Physics is the same as muon plus physics
  pManager = G4Muonium::Muonium()->GetProcessManager();

  pManager->AddProcess(&fMuAtRestSpinRotation,     1, -1, -1);
  pManager->AddProcess(&fMuPlusIonisation,        -1,  2, 2);
  pManager->AddProcess(&fMuoniumScatt,            -1, -1, 3);
  pManager->AddProcess(&fMuPlusBremsstrahlung,    -1,  4, 4);
  pManager->AddProcess(&fMuPlusPairProduction,    -1,  5, 5);

  pManager->AddProcess(new G4UserSpecialCuts(),-1,-1,6); 

  // the last process: decay
  pManager->AddProcess(&fDecayProcess);
  pManager->SetProcessOrderingToLast(&fDecayProcess, idxAtRest);
  pManager->SetProcessOrdering(&fDecayProcess, idxPostStep);
 
 //! Muon Minus Physics
  pManager = G4MuonMinus::MuonMinus()->GetProcessManager();
   
  pManager->AddProcess(&fMuMinusMultipleScattering,-1,  1, 1);
  pManager->AddProcess(&fMuMinusIonisation,        -1,  2, 2);
  pManager->AddProcess(&fMuMinusBremsstrahlung,    -1,  3, 3);
  pManager->AddProcess(&fMuMinusPairProduction,    -1,  4, 4);

  pManager->AddRestProcess(&fMuMinusCaptureAtRest);
  pManager->AddProcess(&fDecayProcess);
  pManager->SetProcessOrderingToLast(&fDecayProcess, idxAtRest);
  pManager->SetProcessOrdering(&fDecayProcess, idxPostStep);

  //! Tau Plus Physics
  pManager = G4TauPlus::TauPlus()->GetProcessManager();
   
  pManager->AddProcess(&fTauPlusMultipleScattering, -1, 1, 1);
  pManager->AddProcess(&fTauPlusIonisation,         -1, 2, 2);

  //! Tau Minus Physics
  pManager = G4TauMinus::TauMinus()->GetProcessManager();

  pManager->AddProcess(&fTauMinusMultipleScattering, -1, 1, 1);
  pManager->AddProcess(&fTauMinusIonisation,         -1, 2, 2);

}