// PhysicsList.cc
// modified by d.satoh on Dec.24.2009.
// G4MultipleScattering class is removed.

#include "PhysicsList.hh"

#include "G4ParticleDefinition.hh"
#include "G4ParticleWithCuts.hh"
#include "G4ProcessManager.hh"
#include "G4ProcessVector.hh"
#include "G4ParticleTypes.hh"
//#include "G4ParticleTable.hh"
#include "G4IonTable.hh"

#include "G4UnitsTable.hh"
#include "G4SystemOfUnits.hh"
#include "G4PhysicalConstants.hh"

PhysicsList::PhysicsList() : G4VUserPhysicsList()
{
  defaultCutValue = 1.0 *mm;
  SetVerboseLevel(1);
}

//PhysicsList::~PhysicsList()
//{;}

#include "G4BosonConstructor.hh"
#include "G4LeptonConstructor.hh"
#include "G4MesonConstructor.hh"
#include "G4BaryonConstructor.hh"
#include "G4ShortLivedConstructor.hh"
#include "G4IonConstructor.hh"

// -------------------------------------------------------
void PhysicsList::ConstructParticle()
{
  // Bosons
  G4BosonConstructor::ConstructParticle();

  // Leptons
  G4LeptonConstructor::ConstructParticle();

  // Mesons
  G4MesonConstructor::ConstructParticle();

  // Baryons
  G4BaryonConstructor::ConstructParticle();

  // Shortlived
  G4ShortLivedConstructor::ConstructParticle();

  // Ions
  G4IonConstructor::ConstructParticle();
}
// -------------------------------------------------------

void PhysicsList::ConstructProcess()
{
  AddTransportation();

  ConstructEM();
  ConstructGeneral();
  ConstructHadronProcess();
  ConstructIonProcess();
}

#include "G4ComptonScattering.hh"
#include "G4GammaConversion.hh"
#include "G4PhotoElectricEffect.hh"

#include "G4eMultipleScattering.hh"
#include "G4MuMultipleScattering.hh"
#include "G4hMultipleScattering.hh"

#include "G4UAtomicDeexcitation.hh"
#include "G4NuclearStopping.hh"
#include "G4LossTableManager.hh"

#include "G4eplusAnnihilation.hh"
#include "G4eIonisation.hh"
#include "G4eBremsstrahlung.hh"

#include "G4MuIonisation.hh"
#include "G4MuBremsstrahlung.hh"
#include "G4MuPairProduction.hh"

#include "G4hIonisation.hh"

// --------------------------
void PhysicsList::ConstructEM()
{
  auto theParticleIterator = GetParticleIterator();
  theParticleIterator->reset();
  while( (*theParticleIterator)() ){
    G4ParticleDefinition* particle = theParticleIterator->value();
    G4ProcessManager* pmanager = particle->GetProcessManager();
    G4String particleName = particle->GetParticleName();
    
    if (particleName == "gamma") {
      // gamma
      pmanager->AddDiscreteProcess(new G4GammaConversion());
      pmanager->AddDiscreteProcess(new G4ComptonScattering());      
      pmanager->AddDiscreteProcess(new G4PhotoElectricEffect());

    } else if (particleName == "e-") {
      //electron
      G4VProcess* theeminusMultipleScattering = new G4eMultipleScattering();
      G4VProcess* theeminusIonisation         = new G4eIonisation();
      G4VProcess* theeminusBremsstrahlung     = new G4eBremsstrahlung();
      // add processes
      pmanager->AddProcess(theeminusMultipleScattering);
      pmanager->AddProcess(theeminusIonisation);
      pmanager->AddProcess(theeminusBremsstrahlung);   
      // set ordering for AlongStepDoIt
      pmanager->SetProcessOrdering(theeminusMultipleScattering, idxAlongStep,1);
      pmanager->SetProcessOrdering(theeminusIonisation,         idxAlongStep,2);
      // set ordering for PostStepDoIt
      pmanager->SetProcessOrdering(theeminusMultipleScattering, idxPostStep,1);
      pmanager->SetProcessOrdering(theeminusIonisation,         idxPostStep,2);
      pmanager->SetProcessOrdering(theeminusBremsstrahlung,     idxPostStep,3);
      
    } else if (particleName == "e+") {
      //positron
      G4VProcess* theeplusMultipleScattering = new G4eMultipleScattering();
      G4VProcess* theeplusIonisation         = new G4eIonisation();
      G4VProcess* theeplusBremsstrahlung     = new G4eBremsstrahlung();
      G4VProcess* theeplusAnnihilation       = new G4eplusAnnihilation();
      // add processes
      pmanager->AddProcess(theeplusMultipleScattering);
      pmanager->AddProcess(theeplusIonisation);
      pmanager->AddProcess(theeplusBremsstrahlung);
      pmanager->AddProcess(theeplusAnnihilation);
      // set ordering for AtRestDoIt
      pmanager->SetProcessOrderingToFirst(theeplusAnnihilation, idxAtRest);
      // set ordering for AlongStepDoIt
      pmanager->SetProcessOrdering(theeplusMultipleScattering, idxAlongStep,1);
      pmanager->SetProcessOrdering(theeplusIonisation,         idxAlongStep,2);
      // set ordering for PostStepDoIt
      pmanager->SetProcessOrdering(theeplusMultipleScattering, idxPostStep,1);
      pmanager->SetProcessOrdering(theeplusIonisation,         idxPostStep,2);
      pmanager->SetProcessOrdering(theeplusBremsstrahlung,     idxPostStep,3);
      pmanager->SetProcessOrdering(theeplusAnnihilation,       idxPostStep,4);
  
    } else if( particleName == "mu+" || particleName == "mu-" ) {
      //muon  
      G4VProcess* aMultipleScattering = new G4MuMultipleScattering();
      G4VProcess* aBremsstrahlung     = new G4MuBremsstrahlung();
      G4VProcess* aPairProduction     = new G4MuPairProduction();
      G4VProcess* anIonisation        = new G4MuIonisation();
      // add processes
      pmanager->AddProcess(anIonisation);
      pmanager->AddProcess(aMultipleScattering);
      pmanager->AddProcess(aBremsstrahlung);
      pmanager->AddProcess(aPairProduction);
      // set ordering for AlongStepDoIt
      pmanager->SetProcessOrdering(aMultipleScattering, idxAlongStep,1);
      pmanager->SetProcessOrdering(anIonisation,        idxAlongStep,2);
      // set ordering for PostStepDoIt
      pmanager->SetProcessOrdering(aMultipleScattering, idxPostStep,1);
      pmanager->SetProcessOrdering(anIonisation,        idxPostStep,2);
      pmanager->SetProcessOrdering(aBremsstrahlung,     idxPostStep,3);
      pmanager->SetProcessOrdering(aPairProduction,     idxPostStep,4);

    } else if ((!particle->IsShortLived()) &&
	       (particle->GetParticleType() != "nucleus") &&
	       (particle->GetPDGCharge() != 0.0) &&
	       (particle->GetParticleName() != "chargedgeantino")) {
     // all others charged particles except geantino     
     G4VProcess* aMultipleScattering = new G4hMultipleScattering();
     G4VProcess* anIonisation        = new G4hIonisation();
     // add processes
     pmanager->AddProcess(anIonisation);
     pmanager->AddProcess(aMultipleScattering);
     // set ordering for AlongStepDoIt
     pmanager->SetProcessOrdering(aMultipleScattering, idxAlongStep,1);
     pmanager->SetProcessOrdering(anIonisation,        idxAlongStep,2);
     // set ordering for PostStepDoIt
     pmanager->SetProcessOrdering(aMultipleScattering, idxPostStep,1);
     pmanager->SetProcessOrdering(anIonisation,        idxPostStep,2);
     }
  }
  // Deexcitation
  G4VAtomDeexcitation *de = new G4UAtomicDeexcitation();
  G4LossTableManager::Instance()->SetAtomDeexcitation(de);
  de->SetFluo(true);
}

#include "G4Decay.hh"
#include "G4RadioactiveDecay.hh"

void PhysicsList::ConstructGeneral()
{
  // Add Decay Process
  G4Decay* theDecayProcess = new G4Decay();
  auto theParticleIterator = GetParticleIterator();
  theParticleIterator->reset();
  while( (*theParticleIterator)() ){
    G4ParticleDefinition* particle = theParticleIterator->value();
    G4ProcessManager* pmanager = particle->GetProcessManager();
    if (theDecayProcess->IsApplicable(*particle)) { 
      pmanager->AddProcess(theDecayProcess);
      // set ordering for PostStepDoIt and AtRestDoIt
      pmanager->SetProcessOrdering(theDecayProcess, idxPostStep);
      pmanager->SetProcessOrdering(theDecayProcess, idxAtRest);
    }
  }
  const G4IonTable *theIonTable = G4IonTable::GetIonTable()->GetIonTable();
  G4RadioactiveDecay *theRadioactiveDecay = new G4RadioactiveDecay();

  for (G4int i=0; i<theIonTable->Entries(); i++) 
    {
      G4ParticleDefinition *particle = theIonTable->GetParticle(i);
      G4String particleName = particle->GetParticleName();
      
      if (particleName == "GenericIon") {
          G4ProcessManager* pmanager = particle->GetProcessManager();
          pmanager->AddProcess(theRadioactiveDecay);
          pmanager->SetProcessOrdering(theRadioactiveDecay, idxPostStep);
          pmanager->SetProcessOrdering(theRadioactiveDecay, idxAtRest);
        } 
    }
}

#include "G4GenericIon.hh"
#include "G4HadronElasticProcess.hh"
//#include "G4LElastic.hh"
#include "G4HadronElastic.hh"

#include "G4HadronCaptureProcess.hh"
#include "G4NeutronRadCapture.hh"
#include "G4NeutronCaptureXS.hh"
#include "G4NeutronHPCaptureData.hh"

#include "G4HadronInelasticProcess.hh"
#include "G4HadronInelasticQBBC.hh"
#include "G4NeutronInelasticXS.hh"
#include "G4BGGNucleonInelasticXS.hh"

G4ThreadLocal PhysicsList::ThreadPrivate* 
PhysicsList::tpdata=0;

void PhysicsList::ConstructHadron()
{
  G4bool quasiElasticFTF= false;   // Use built-in quasi-elastic (not add-on)
  G4bool quasiElasticQGS= true;    // For QGS, it must use it.

  const G4double maxFTFP = 25.0*GeV;
  const G4double minFTFP =  9.5*GeV;
  const G4double maxBIC  =  9.9*GeV;
  const G4double maxBERT =  5.0*GeV;
  const G4double maxHP   = 19.9*MeV;

  tpdata->theNeutrons=new G4NeutronBuilder( false ); // Fission off
  tpdata->theNeutrons->RegisterMe(tpdata->theQGSPNeutron=new G4QGSPNeutronBuilder(quasiElasticQGS));
  tpdata->theNeutrons->RegisterMe(tpdata->theFTFPNeutron=new G4FTFPNeutronBuilder(quasiElasticFTF));
  tpdata->theFTFPNeutron->SetMinEnergy(minFTFP);
  tpdata->theFTFPNeutron->SetMaxEnergy(maxFTFP);

  tpdata->theNeutrons->RegisterMe(tpdata->theBinaryNeutron=new G4BinaryNeutronBuilder);
  tpdata->theBinaryNeutron->SetMinEnergy(maxHP);
  tpdata->theBinaryNeutron->SetMaxEnergy(maxBIC);

  tpdata->theNeutrons->RegisterMe(tpdata->theHPNeutron=new G4NeutronHPBuilder);

  tpdata->thePro=new G4ProtonBuilder;
  tpdata->thePro->RegisterMe(tpdata->theQGSPPro=new G4QGSPProtonBuilder(quasiElasticQGS));
  tpdata->thePro->RegisterMe(tpdata->theFTFPPro=new G4FTFPProtonBuilder(quasiElasticFTF));
  tpdata->theFTFPPro->SetMinEnergy(minFTFP);
  tpdata->theFTFPPro->SetMaxEnergy(maxFTFP);
  tpdata->thePro->RegisterMe(tpdata->theBinaryPro=new G4BinaryProtonBuilder);
  tpdata->theBinaryPro->SetMaxEnergy(maxBIC);
  
  tpdata->thePiK=new G4PiKBuilder;
  tpdata->thePiK->RegisterMe(tpdata->theQGSPPiK=new G4QGSPPiKBuilder(quasiElasticQGS));
  tpdata->thePiK->RegisterMe(tpdata->theFTFPPiK=new G4FTFPPiKBuilder(quasiElasticFTF));
  tpdata->theFTFPPiK->SetMaxEnergy(maxFTFP);
  tpdata->thePiK->RegisterMe(tpdata->theBertiniPiK=new G4BertiniPiKBuilder);
  tpdata->theBertiniPiK->SetMaxEnergy(maxBERT);
}            

PhysicsList::~PhysicsList() 
{
   delete tpdata->theHPNeutron;
   delete tpdata->theBinaryNeutron;
   delete tpdata->theFTFPNeutron;
   delete tpdata->theQGSPNeutron;
   delete tpdata->theNeutrons;
   delete tpdata->theBertiniPiK;
   delete tpdata->theFTFPPiK;
   delete tpdata->theQGSPPiK;
   delete tpdata->thePiK;
   delete tpdata->theBinaryPro;
   delete tpdata->theFTFPPro;
   delete tpdata->theQGSPPro;
   delete tpdata->thePro;
   delete tpdata->xsNeutronCaptureXS;

  //Note that here we need to set to 0 the pointer
  //since tpdata is static and if thread are "reused"
  //it can be problematic
   delete tpdata; tpdata = 0;
}

#include "G4ChipsElasticModel.hh"
#include "G4CrossSectionDataSetRegistry.hh"
#include "G4ChipsNeutronElasticXS.hh"
#include "G4ChipsProtonElasticXS.hh"

void PhysicsList::ConstructHadronProcess()
{
  G4ProcessManager* pManager = 0;
  if ( tpdata == 0 ) tpdata = new ThreadPrivate;
  ConstructHadron();
  tpdata->theNeutrons->Build();
  tpdata->thePro->Build();
  tpdata->thePiK->Build();
  
  // Proton Elastic Process
  G4ChipsElasticModel* ProtonElasticModel = new G4ChipsElasticModel();
  G4HadronElasticProcess* ProtonElasticProcess = new G4HadronElasticProcess();
  ProtonElasticProcess->AddDataSet(G4CrossSectionDataSetRegistry::Instance()->GetCrossSectionDataSet(G4ChipsProtonElasticXS::Default_Name()));
  ProtonElasticProcess->RegisterMe(ProtonElasticModel);

  pManager = G4Proton::Proton()->GetProcessManager();
  pManager->AddDiscreteProcess(ProtonElasticProcess);
  
  // Neutron Elastic Process
  G4ChipsElasticModel* NeutronElasticModel = new G4ChipsElasticModel();
  G4HadronElasticProcess* NeutronElasticProcess = new G4HadronElasticProcess();
  NeutronElasticProcess->AddDataSet(G4CrossSectionDataSetRegistry::Instance()->GetCrossSectionDataSet(G4ChipsNeutronElasticXS::Default_Name()));
  NeutronElasticProcess->RegisterMe(NeutronElasticModel);
    
  // Neutron Capture
  G4HadronicProcess* capture = 0;
  pManager = G4Neutron::Neutron()->GetProcessManager();
  pManager->AddDiscreteProcess(NeutronElasticProcess);
  G4ProcessVector*  pv = pManager->GetProcessList();
  for ( size_t i=0; i < static_cast<size_t>(pv->size()); ++i ) {
    if ( fCapture == ((*pv)[i])->GetProcessSubType() ) {
      capture = static_cast<G4HadronicProcess*>((*pv)[i]);
    }
  }  
  if ( ! capture ) {
    capture = new G4HadronCaptureProcess("nCapture");
    pManager->AddDiscreteProcess(capture);
  }
  tpdata->xsNeutronCaptureXS = new G4NeutronCaptureXS();
  capture->AddDataSet(tpdata->xsNeutronCaptureXS);
  capture->AddDataSet( new G4NeutronHPCaptureData );
  G4NeutronRadCapture* theNeutronRadCapture = new G4NeutronRadCapture(); 
  theNeutronRadCapture->SetMinEnergy(19.9*MeV); 
  capture->RegisterMe( theNeutronRadCapture );
}

#include "G4DeuteronInelasticProcess.hh"
#include "G4TritonInelasticProcess.hh"
#include "G4AlphaInelasticProcess.hh"

#include "G4hIonisation.hh"
#include "G4ionIonisation.hh"
#include "G4hMultipleScattering.hh"

#include "G4BinaryLightIonReaction.hh"
#include "G4QMDReaction.hh"
#include "G4FTFBuilder.hh"
#include "G4HadronicInteraction.hh"
#include "G4HadronicInteractionRegistry.hh"

#include "G4TripathiCrossSection.hh"
#include "G4IonsShenCrossSection.hh"

#include "G4ComponentGGNuclNuclXsc.hh"
#include "G4CrossSectionInelastic.hh"
#include "G4CrossSectionElastic.hh"

void PhysicsList::ConstructIonProcess()
{
  G4ProcessManager* pManager = 0;

  const G4double emaxFTFP = 1.*TeV;
  const G4double emaxQMD =  10.*GeV;
  const G4double eminQMD  =  10.*MeV;
  const G4double eminBIC =  0.*MeV;
  const G4double overlap   = 10*MeV;

  G4HadronicInteraction* p = G4HadronicInteractionRegistry::Instance()->FindModel("PRECO");
  G4PreCompoundModel* thePreCompound = static_cast<G4PreCompoundModel*>(p); 
  if(!thePreCompound) { thePreCompound = new G4PreCompoundModel; }

  // MCS and ionization
  G4hMultipleScattering* fIonMultipleScattering=new G4hMultipleScattering();
  G4ionIonisation* fIonIonisation=new G4ionIonisation();

  // Binary Light Ion Cascade
  G4BinaryLightIonReaction* theBLIModel= new G4BinaryLightIonReaction(thePreCompound);
  theBLIModel->SetMinEnergy(eminBIC);
  theBLIModel->SetMaxEnergy(eminQMD+overlap);
  
  // QMD
  G4QMDReaction* theQMDModel= new G4QMDReaction; 
  theQMDModel->SetMinEnergy(eminQMD);
  theQMDModel->SetMaxEnergy(emaxQMD);
  //theQMDModel->UnUseGEM(); // then use the default channel of G4Evaporation. 
  //theQMDModel->UseFRAG(); // change criterion of a inelastic reaction
  //theQMDModel->SetTMAX(100); // will have maxTime-th time step, default = 100
  //theQMDModel->SetDT(1); // in fsec (c=1), default = 1
  //theQMDModel->SetEF(1.05); // 10% for Peripheral reactions, default = 1.05
  
  // FTFP
  G4FTFBuilder* theFTFBuilder= new G4FTFBuilder("FTFP",thePreCompound);
  G4HadronicInteraction* theFTFPModel= theFTFBuilder->GetModel();
  theFTFPModel->SetMinEnergy(emaxQMD-overlap);
  theFTFPModel->SetMaxEnergy(emaxFTFP);
  
  //Total Cross Section
  G4VComponentCrossSection* theGGNuclNuclXS = new G4ComponentGGNuclNuclXsc();
  G4VCrossSectionDataSet* theNuclNuclInelaData = new G4CrossSectionInelastic(theGGNuclNuclXS);
  G4VCrossSectionDataSet* theNuclNuclElaData = new G4CrossSectionElastic(theGGNuclNuclXS);
  //G4TripathiCrossSection* TripathiCrossSection= new G4TripathiCrossSection;
  //G4IonsShenCrossSection* aShen = new G4IonsShenCrossSection;

  // Elastic Process
  G4HadronElastic* theElasticModel = new G4HadronElastic();
  G4HadronElasticProcess* fElasticProcess = new G4HadronElasticProcess();
  fElasticProcess->RegisterMe(theElasticModel);
  fElasticProcess->AddDataSet(theNuclNuclElaData);

  // Generic Ion
  pManager = G4GenericIon::GenericIon()->GetProcessManager();
  pManager->AddDiscreteProcess(fElasticProcess);
  G4HadronInelasticProcess* theIPGenericIon = new G4HadronInelasticProcess("IonInelastic", G4GenericIon::GenericIon());
  //theIPGenericIon->AddDataSet(TripathiCrossSection);
  //theIPGenericIon->AddDataSet(aShen);
  theIPGenericIon->AddDataSet(theNuclNuclInelaData);
  theIPGenericIon->RegisterMe(theBLIModel);
  theIPGenericIon->RegisterMe(theQMDModel);
  theIPGenericIon->RegisterMe(theFTFPModel);
  pManager->AddDiscreteProcess(theIPGenericIon);
  pManager->AddProcess(fIonMultipleScattering, -1, 1, 1);
  pManager->AddProcess(fIonIonisation, -1, 2, 2);

  // Deuteron 
  pManager = G4Deuteron::Deuteron()->GetProcessManager();
  pManager->AddDiscreteProcess(fElasticProcess);
  G4DeuteronInelasticProcess*  fDeuteronProcess = new G4DeuteronInelasticProcess();
  //fDeuteronProcess->AddDataSet(TripathiCrossSection);
  //fDeuteronProcess->AddDataSet(aShen);
  fDeuteronProcess->AddDataSet(theNuclNuclInelaData);
  fDeuteronProcess->RegisterMe(theBLIModel);
  fDeuteronProcess->RegisterMe(theQMDModel);
  fDeuteronProcess->RegisterMe(theFTFPModel);
  pManager->AddDiscreteProcess(fElasticProcess);
  pManager->AddDiscreteProcess(fDeuteronProcess);
  pManager->AddProcess(fIonMultipleScattering, -1, 1, 1);
  pManager->AddProcess(fIonIonisation, -1, 2, 2);

  // Triton 
  pManager = G4Triton::Triton()->GetProcessManager();
  pManager->AddDiscreteProcess(fElasticProcess);
  G4TritonInelasticProcess* fTritonProcess = new G4TritonInelasticProcess();
  //fTritonProcess->AddDataSet(TripathiCrossSection);
  //fTritonProcess->AddDataSet(aShen);
  fTritonProcess->AddDataSet(theNuclNuclInelaData);
  fTritonProcess->RegisterMe(theBLIModel);
  fTritonProcess->RegisterMe(theQMDModel);
  fTritonProcess->RegisterMe(theFTFPModel);
  pManager->AddDiscreteProcess(fElasticProcess);
  pManager->AddDiscreteProcess(fTritonProcess);
  pManager->AddProcess(fIonMultipleScattering, -1, 1, 1);
  pManager->AddProcess(fIonIonisation, -1, 2, 2);
 
  // Alpha 
  pManager = G4Alpha::Alpha()->GetProcessManager();
  pManager->AddDiscreteProcess(fElasticProcess);
  G4AlphaInelasticProcess* fAlphaProcess = new G4AlphaInelasticProcess();
  //fAlphaProcess->AddDataSet(TripathiCrossSection);
  //fAlphaProcess->AddDataSet(aShen);
  fAlphaProcess->AddDataSet(theNuclNuclInelaData);
  fAlphaProcess->RegisterMe(theBLIModel);
  fAlphaProcess->RegisterMe(theQMDModel);
  fAlphaProcess->RegisterMe(theFTFPModel);
  pManager->AddDiscreteProcess(fElasticProcess);
  pManager->AddDiscreteProcess(fAlphaProcess);
  pManager->AddProcess(fIonMultipleScattering, -1, 1, 1);
  pManager->AddProcess(fIonIonisation, -1, 2, 2);

  // He3
  pManager = G4He3::He3()->GetProcessManager();
  pManager->AddDiscreteProcess(fElasticProcess);
  G4HadronInelasticProcess* fHe3Process = new G4HadronInelasticProcess("He3Inelastic", G4He3::He3());
  //fHe3Process->AddDataSet(TripathiCrossSection);
  //fHe3Process->AddDataSet(aShen);
  fHe3Process->AddDataSet(theNuclNuclInelaData);
  fHe3Process->RegisterMe(theBLIModel);
  fHe3Process->RegisterMe(theQMDModel);
  fHe3Process->RegisterMe(theFTFPModel);
  pManager->AddDiscreteProcess(fHe3Process); 
  pManager->AddProcess(fIonMultipleScattering, -1, 1, 1);
  pManager->AddProcess(fIonIonisation, -1, 2, 2);   
}

// -------------------------------------------------------
void PhysicsList::SetCuts()
{
  SetCutsWithDefault();
}