//this
#include "ELYSE/PhysicsList.hh"

//Geant 4
#include "G4EmStandardPhysics.hh"
#include "G4EmStandardPhysics71.hh"
#include "G4LossTableManager.hh"
#include "G4UnitsTable.hh"

#include "G4Cerenkov.hh"
#include "G4Scintillation.hh"
#include "G4OpAbsorption.hh"
#include "G4OpRayleigh.hh"
#include "G4OpBoundaryProcess.hh"


// Bosons
#include "G4ChargedGeantino.hh"
#include "G4Geantino.hh"
#include "G4Gamma.hh"
#include "G4OpticalPhoton.hh"

// leptons
#include "G4MuonPlus.hh"
#include "G4MuonMinus.hh"
#include "G4NeutrinoMu.hh"
#include "G4AntiNeutrinoMu.hh"

#include "G4Electron.hh"
#include "G4Positron.hh"
#include "G4NeutrinoE.hh"
#include "G4AntiNeutrinoE.hh"

// Mesons
#include "G4PionPlus.hh"
#include "G4PionMinus.hh"
#include "G4PionZero.hh"
#include "G4Eta.hh"
#include "G4EtaPrime.hh"

#include "G4KaonPlus.hh"
#include "G4KaonMinus.hh"
#include "G4KaonZero.hh"
#include "G4AntiKaonZero.hh"
#include "G4KaonZeroLong.hh"
#include "G4KaonZeroShort.hh"

// Baryons
#include "G4Proton.hh"
#include "G4AntiProton.hh"
#include "G4Neutron.hh"
#include "G4AntiNeutron.hh"

// Nuclei
#include "G4Deuteron.hh"
#include "G4Triton.hh"
#include "G4Alpha.hh"
#include "G4GenericIon.hh"

#include "G4Gamma.hh"
#include "G4Electron.hh"
#include "G4Positron.hh"


//ELYSE
#include "ELYSE/PhysicsListMessenger.hh"
#include "ELYSE/PhysListEmStandard.hh"
#include "ELYSE/PhysListEmLivermore.hh"
#include "ELYSE/PhysListEmPenelope.hh"


//----------------------------------------------------------------------

ELYSE::PhysicsList::PhysicsList() : G4VModularPhysicsList(),
  theCerenkovProcess(0),
  theScintillationProcess(0),
  theAbsorptionProcess(0),
  theRayleighScatteringProcess(0),
  theBoundaryProcess(0)
 {
  G4LossTableManager::Instance();
  currentDefaultCut   = 1.0*mm;
  cutForGamma         = currentDefaultCut;
  cutForElectron      = currentDefaultCut;
  cutForPositron      = currentDefaultCut;

  pMessenger = new PhysicsListMessenger(this);

  SetVerboseLevel(1);

  // EM physics
  emName = G4String("standard");
  emPhysicsList = new PhysListEmStandard(emName);

}

//----------------------------------------------------------------------

ELYSE::PhysicsList::~PhysicsList() {
  delete pMessenger;
}

//----------------------------------------------------------------------


void ELYSE::PhysicsList::ConstructParticle() {
// pseudo-particles
  G4Geantino::GeantinoDefinition();
  G4ChargedGeantino::ChargedGeantinoDefinition();
  
// gamma
  G4Gamma::GammaDefinition();
  
// optical photon
  G4OpticalPhoton::OpticalPhotonDefinition();

// leptons
  G4Electron::ElectronDefinition();
  G4Positron::PositronDefinition();
  G4MuonPlus::MuonPlusDefinition();
  G4MuonMinus::MuonMinusDefinition();

  G4NeutrinoE::NeutrinoEDefinition();
  G4AntiNeutrinoE::AntiNeutrinoEDefinition();
  G4NeutrinoMu::NeutrinoMuDefinition();
  G4AntiNeutrinoMu::AntiNeutrinoMuDefinition();  

// mesons
  G4PionPlus::PionPlusDefinition();
  G4PionMinus::PionMinusDefinition();
  G4PionZero::PionZeroDefinition();
  G4Eta::EtaDefinition();
  G4EtaPrime::EtaPrimeDefinition();
  G4KaonPlus::KaonPlusDefinition();
  G4KaonMinus::KaonMinusDefinition();
  G4KaonZero::KaonZeroDefinition();
  G4AntiKaonZero::AntiKaonZeroDefinition();
  G4KaonZeroLong::KaonZeroLongDefinition();
  G4KaonZeroShort::KaonZeroShortDefinition();

// barions
  G4Proton::ProtonDefinition();
  G4AntiProton::AntiProtonDefinition();
  G4Neutron::NeutronDefinition();
  G4AntiNeutron::AntiNeutronDefinition();

// ions
  G4Deuteron::DeuteronDefinition();
  G4Triton::TritonDefinition();
  G4Alpha::AlphaDefinition();
  G4GenericIon::GenericIonDefinition();
}

//----------------------------------------------------------------------

#include "G4ProcessManager.hh"
#include "G4Decay.hh"

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

  // electromagnetic Physics List
  //
  emPhysicsList->ConstructProcess();

  //Decay
  ConstructGeneral();

  //Optical Photons
  ConstructOp();
  //JEC 20/2/07 not useful? ConstructHad();
  SetVerbose(0);
    
}//ConstructProcess

//--------------------------------------------------------------------------------------------------

void ELYSE::PhysicsList::ConstructGeneral() {
  // Add Decay Process
  G4Decay* theDecayProcess = new G4Decay();
  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);
    }
  }
}//ConstructGeneral

//-----------------------------------------------------------------------------------------------
void ELYSE::PhysicsList::ConstructOp(){

  theCerenkovProcess           = new G4Cerenkov("Cerenkov");
  theScintillationProcess      = new G4Scintillation("Scintillation");
  theAbsorptionProcess         = new G4OpAbsorption();
  theRayleighScatteringProcess = new G4OpRayleigh();
  theBoundaryProcess           = new G4OpBoundaryProcess();

  /*
    theCerenkovProcess->DumpPhysicsTable();
    theScintillationProcess->DumpPhysicsTable();
    theAbsorptionProcess->DumpPhysicsTable();
    theRayleighScatteringProcess->DumpPhysicsTable();
  */


  theCerenkovProcess->SetTrackSecondariesFirst(false); //NV: if "true" generates as many trajectory objs of  
  theCerenkovProcess->SetMaxNumPhotonsPerStep(300);    //op parents as the # of parents' steps. 
                                                       //The number of photons may be tuned by the messenger

  //JEC 19/2/07  theScintillationProcess->SetScintillationYieldFactor(1.);
  theScintillationProcess->SetScintillationYieldFactor(0.);
  theScintillationProcess->SetTrackSecondariesFirst(false);

  G4OpticalSurfaceModel themodel = unified; //JEC where does it come from?
  theBoundaryProcess->SetModel(themodel);

  theParticleIterator->reset();
  while( (*theParticleIterator)() )
  {
    G4ParticleDefinition* particle     = theParticleIterator->value();
    G4ProcessManager*     pmanager     = particle->GetProcessManager();
    G4String              particleName = particle->GetParticleName();

    if (theCerenkovProcess->IsApplicable(*particle)) {
      pmanager->AddContinuousProcess(theCerenkovProcess);
    }
    
    if (theScintillationProcess->IsApplicable(*particle)) {
      pmanager->AddProcess(theScintillationProcess);
      pmanager->SetProcessOrderingToLast(theScintillationProcess, idxAtRest);
      pmanager->SetProcessOrderingToLast(theScintillationProcess, idxPostStep);
    }

    if (particleName == "opticalphoton") {
      pmanager->AddDiscreteProcess(theAbsorptionProcess);
      pmanager->AddDiscreteProcess(theRayleighScatteringProcess);
      pmanager->AddDiscreteProcess(theBoundaryProcess);
    }
  }
}//ConstructOp

//----------------------------------------------------------------------

void ELYSE::PhysicsList::AddPhysicsList(const G4String& name) {
  if (verboseLevel>1) {
    G4cout << "ELYSE::PhysicsList::AddPhysicsList: <" << name << ">" << G4endl;
  }

  if (name == emName) return;

  if (name == "standard") {

    emName = name;
    delete emPhysicsList;
    emPhysicsList = new PhysListEmStandard(name);

  } else if (name == "G4standard") {

    emName = name;
    delete emPhysicsList;
    emPhysicsList = new G4EmStandardPhysics(name);

  } else if (name == "G4standard_fast") {
 
    emName = name;
    delete emPhysicsList;
    emPhysicsList = new G4EmStandardPhysics71(name);

  } else if (name == "Livermore") {

    emName = name;
    delete emPhysicsList;
    emPhysicsList = new PhysListEmLivermore(name);

  } else if (name == "Penelope") {

    emName = name;
    delete emPhysicsList;
    emPhysicsList = new PhysListEmPenelope(name);

  } else {

    G4cout << "ELYSE::PhysicsList::AddPhysicsList: <" << name << ">"
           << " is not defined"
           << G4endl;
  }
}

//----------------------------------------------------------------------
void ELYSE::PhysicsList::SetCuts()
{
     
  if (verboseLevel >0){
    G4cout << "PhysicsList::SetCuts:";
    G4cout << "CutLength : " << G4BestUnit(defaultCutValue,"Length") << G4endl;
  }  

  // set cut values for gamma at first and for e- second and next for e+,
  // because some processes for e+/e- need cut values for gamma
  SetCutValue(cutForGamma, "gamma");
  SetCutValue(cutForElectron, "e-");
  SetCutValue(cutForPositron, "e+");   
    
  if (verboseLevel>0) DumpCutValuesTable();
}

//----------------------------------------------------------------------

void ELYSE::PhysicsList::SetCutForGamma(G4double cut) {
  cutForGamma = cut;
  SetParticleCuts(cutForGamma, G4Gamma::Gamma());
}

//----------------------------------------------------------------------

void ELYSE::PhysicsList::SetCutForElectron(G4double cut) {
  cutForElectron = cut;
  SetParticleCuts(cutForElectron, G4Electron::Electron());
}

//----------------------------------------------------------------------

void ELYSE::PhysicsList::SetCutForPositron(G4double cut) {
  cutForPositron = cut;
  SetParticleCuts(cutForPositron, G4Positron::Positron());
}

//-------------------------------------------------------------------
void ELYSE::PhysicsList::SetVerbose(G4int verbose) {

  theCerenkovProcess->SetVerboseLevel(verbose);
  theScintillationProcess->SetVerboseLevel(verbose);
  theAbsorptionProcess->SetVerboseLevel(verbose);
  theRayleighScatteringProcess->SetVerboseLevel(verbose);
  theBoundaryProcess->SetVerboseLevel(verbose);  

}//SetVerbose
//-------------------------------------------------------------------

void ELYSE::PhysicsList::SetNbOfPhotonsCerenkov(G4int MaxNumber) {  
  theCerenkovProcess->SetMaxNumPhotonsPerStep(MaxNumber);
}//SetNbOfPhotonsCerenkov


//-----------------------------------------------------------------------------------------------

// // Elastic processes:
// #include "G4HadronElasticProcess.hh"

// // Inelastic processes:
// #include "G4PionPlusInelasticProcess.hh"
// #include "G4PionMinusInelasticProcess.hh"
// #include "G4KaonPlusInelasticProcess.hh"
// #include "G4KaonZeroSInelasticProcess.hh"
// #include "G4KaonZeroLInelasticProcess.hh"
// #include "G4KaonMinusInelasticProcess.hh"
// #include "G4ProtonInelasticProcess.hh"
// #include "G4AntiProtonInelasticProcess.hh"
// #include "G4NeutronInelasticProcess.hh"
// #include "G4AntiNeutronInelasticProcess.hh"
// #include "G4DeuteronInelasticProcess.hh"
// #include "G4TritonInelasticProcess.hh"
// #include "G4AlphaInelasticProcess.hh"

// // Low-energy Models: < 20GeV
// #include "G4LElastic.hh"
// #include "G4LEPionPlusInelastic.hh"
// #include "G4LEPionMinusInelastic.hh"
// #include "G4LEKaonPlusInelastic.hh"
// #include "G4LEKaonZeroSInelastic.hh"
// #include "G4LEKaonZeroLInelastic.hh"
// #include "G4LEKaonMinusInelastic.hh"
// #include "G4LEProtonInelastic.hh"
// #include "G4LEAntiProtonInelastic.hh"
// #include "G4LENeutronInelastic.hh"
// #include "G4LEAntiNeutronInelastic.hh"
// #include "G4LEDeuteronInelastic.hh"
// #include "G4LETritonInelastic.hh"
// #include "G4LEAlphaInelastic.hh"

// // High-energy Models: >20 GeV
// #include "G4HEPionPlusInelastic.hh"
// #include "G4HEPionMinusInelastic.hh"
// #include "G4HEKaonPlusInelastic.hh"
// #include "G4HEKaonZeroInelastic.hh"
// #include "G4HEKaonZeroInelastic.hh"
// #include "G4HEKaonMinusInelastic.hh"
// #include "G4HEProtonInelastic.hh"
// #include "G4HEAntiProtonInelastic.hh"
// #include "G4HENeutronInelastic.hh"
// #include "G4HEAntiNeutronInelastic.hh"

// // Neutron high-precision models: <20 MeV
// #include "G4NeutronHPElastic.hh"
// #include "G4NeutronHPElasticData.hh"
// #include "G4NeutronHPCapture.hh"
// #include "G4NeutronHPCaptureData.hh"
// #include "G4NeutronHPInelastic.hh"
// #include "G4NeutronHPInelasticData.hh"
// #include "G4LCapture.hh"

// // Secondary hadronic interaction models
// #include "G4CascadeInterface.hh"
// #include "G4BinaryCascade.hh"

// // Stopping processes
// #include "G4PiMinusAbsorptionAtRest.hh"
// #include "G4KaonMinusAbsorptionAtRest.hh"
// #include "G4AntiProtonAnnihilationAtRest.hh"
// #include "G4AntiNeutronAnnihilationAtRest.hh"

// void ELYSE::PhysicsList::ConstructHad() {

// // This code stolen from:
// // examples/advanced/underground_physics/src/DMXPhysicsList.cc
// // CWW 2/23/05
// //
  
//   G4HadronElasticProcess* theElasticProcess = new G4HadronElasticProcess;
//   G4LElastic* theElasticModel = new G4LElastic;
//   theElasticProcess->RegisterMe(theElasticModel);
  
//   theParticleIterator->reset();
//   while ((*theParticleIterator)()) 
//     {
//       G4ParticleDefinition* particle = theParticleIterator->value();
//       G4ProcessManager* pmanager = particle->GetProcessManager();
//       G4String particleName = particle->GetParticleName();

//       if (particleName == "pi+") 
// 	{
// 	  pmanager->AddDiscreteProcess(theElasticProcess);
// 	  G4PionPlusInelasticProcess* theInelasticProcess = 
// 	    new G4PionPlusInelasticProcess();
// 	  G4LEPionPlusInelastic* theLEInelasticModel = 
// 	    new G4LEPionPlusInelastic;
// 	  theInelasticProcess->RegisterMe(theLEInelasticModel);
// 	  G4HEPionPlusInelastic* theHEInelasticModel = 
// 	    new G4HEPionPlusInelastic;
// 	  theInelasticProcess->RegisterMe(theHEInelasticModel);
// 	  pmanager->AddDiscreteProcess(theInelasticProcess);
// 	} 

//       else if (particleName == "pi-") 
// 	{
// 	  pmanager->AddDiscreteProcess(theElasticProcess);
// 	  G4PionMinusInelasticProcess* theInelasticProcess = 
// 	    new G4PionMinusInelasticProcess();
// 	  G4LEPionMinusInelastic* theLEInelasticModel = 
// 	    new G4LEPionMinusInelastic;
// 	  theInelasticProcess->RegisterMe(theLEInelasticModel);
// 	  G4HEPionMinusInelastic* theHEInelasticModel = 
// 	    new G4HEPionMinusInelastic;
// 	  theInelasticProcess->RegisterMe(theHEInelasticModel);
// 	  pmanager->AddDiscreteProcess(theInelasticProcess);
// 	  G4String prcNam;
// 	  pmanager->AddRestProcess(new G4PiMinusAbsorptionAtRest, ordDefault);
// 	}
      
//       else if (particleName == "kaon+") 
// 	{
// 	  pmanager->AddDiscreteProcess(theElasticProcess);
// 	  G4KaonPlusInelasticProcess* theInelasticProcess = 
// 	    new G4KaonPlusInelasticProcess();
// 	  G4LEKaonPlusInelastic* theLEInelasticModel = 
// 	    new G4LEKaonPlusInelastic;
// 	  theInelasticProcess->RegisterMe(theLEInelasticModel);
// 	  G4HEKaonPlusInelastic* theHEInelasticModel = 
// 	    new G4HEKaonPlusInelastic;
// 	  theInelasticProcess->RegisterMe(theHEInelasticModel);
// 	  pmanager->AddDiscreteProcess(theInelasticProcess);
// 	}
      
//       else if (particleName == "kaon0S") 
// 	{
// 	  pmanager->AddDiscreteProcess(theElasticProcess);
// 	  G4KaonZeroSInelasticProcess* theInelasticProcess = 
// 	    new G4KaonZeroSInelasticProcess();
// 	  G4LEKaonZeroSInelastic* theLEInelasticModel = 
// 	    new G4LEKaonZeroSInelastic;
// 	  theInelasticProcess->RegisterMe(theLEInelasticModel);
// 	  G4HEKaonZeroInelastic* theHEInelasticModel = 
// 	    new G4HEKaonZeroInelastic;
// 	  theInelasticProcess->RegisterMe(theHEInelasticModel);
// 	  pmanager->AddDiscreteProcess(theInelasticProcess);
// 	}

//       else if (particleName == "kaon0L") 
// 	{
// 	  pmanager->AddDiscreteProcess(theElasticProcess);
// 	  G4KaonZeroLInelasticProcess* theInelasticProcess = 
// 	    new G4KaonZeroLInelasticProcess();
// 	  G4LEKaonZeroLInelastic* theLEInelasticModel = 
// 	    new G4LEKaonZeroLInelastic;
// 	  theInelasticProcess->RegisterMe(theLEInelasticModel);
// 	  G4HEKaonZeroInelastic* theHEInelasticModel = 
// 	    new G4HEKaonZeroInelastic;
// 	  theInelasticProcess->RegisterMe(theHEInelasticModel);
// 	  pmanager->AddDiscreteProcess(theInelasticProcess);
// 	}

//       else if (particleName == "kaon-") 
// 	{
// 	  pmanager->AddDiscreteProcess(theElasticProcess);
// 	  G4KaonMinusInelasticProcess* theInelasticProcess = 
// 	    new G4KaonMinusInelasticProcess();
// 	  G4LEKaonMinusInelastic* theLEInelasticModel = 
// 	    new G4LEKaonMinusInelastic;
// 	  theInelasticProcess->RegisterMe(theLEInelasticModel);
// 	  G4HEKaonMinusInelastic* theHEInelasticModel = 
// 	    new G4HEKaonMinusInelastic;
// 	  theInelasticProcess->RegisterMe(theHEInelasticModel);
// 	  pmanager->AddDiscreteProcess(theInelasticProcess);
// 	  pmanager->AddRestProcess(new G4KaonMinusAbsorptionAtRest, ordDefault);
// 	}

//       else if (particleName == "proton") 
// 	{
// 	  pmanager->AddDiscreteProcess(theElasticProcess);
// 	  G4ProtonInelasticProcess* theInelasticProcess = 
// 	    new G4ProtonInelasticProcess();


// 	  //JEC 12/2/07 USE ONE SPECIFIC hadronic model... 
// // 	  if (gheishahad) {
// 	    G4LEProtonInelastic* theLEInelasticModel = new G4LEProtonInelastic;
// 	    theInelasticProcess->RegisterMe(theLEInelasticModel);
// // 	  }
// // 	  else if (bertinihad) {
// // 	    G4CascadeInterface* theBertiniModel = new G4CascadeInterface;
// // 	    theInelasticProcess->RegisterMe(theBertiniModel);
// // 	  }
// // 	  else if (binaryhad) {
// // 	    G4BinaryCascade* theBinaryModel = new G4BinaryCascade();
// // 	    theInelasticProcess->RegisterMe(theBinaryModel);
// // 	  }
// // 	  else {
// // 	    G4cout << "No secondary interaction model chosen! Using G4 BINARY." << G4endl;
// // 	    G4BinaryCascade* theBinaryModel = new G4BinaryCascade();
// // 	    theInelasticProcess->RegisterMe(theBinaryModel);
// // 	  }

// 	  G4HEProtonInelastic* theHEInelasticModel = new G4HEProtonInelastic;
// 	  theInelasticProcess->RegisterMe(theHEInelasticModel);
// 	  pmanager->AddDiscreteProcess(theInelasticProcess);
// 	}

//       else if (particleName == "anti_proton") 
// 	{
// 	  pmanager->AddDiscreteProcess(theElasticProcess);
// 	  G4AntiProtonInelasticProcess* theInelasticProcess = 
// 	    new G4AntiProtonInelasticProcess();
// 	  G4LEAntiProtonInelastic* theLEInelasticModel = 
// 	    new G4LEAntiProtonInelastic;
// 	  theInelasticProcess->RegisterMe(theLEInelasticModel);
// 	  G4HEAntiProtonInelastic* theHEInelasticModel = 
// 	    new G4HEAntiProtonInelastic;
// 	  theInelasticProcess->RegisterMe(theHEInelasticModel);
// 	  pmanager->AddDiscreteProcess(theInelasticProcess);
// 	}

//       else if (particleName == "neutron") 
// 	{
// 	  // elastic scattering
// 	  G4HadronElasticProcess* theNeutronElasticProcess = 
// 	    new G4HadronElasticProcess;
// 	  G4LElastic* theElasticModel1 = new G4LElastic;
// 	  G4NeutronHPElastic * theElasticNeutron = new G4NeutronHPElastic;
// 	  theNeutronElasticProcess->RegisterMe(theElasticModel1);
// 	  theElasticModel1->SetMinEnergy(19*MeV);
// 	  theNeutronElasticProcess->RegisterMe(theElasticNeutron);
// 	  G4NeutronHPElasticData * theNeutronData = new G4NeutronHPElasticData;
// 	  theNeutronElasticProcess->AddDataSet(theNeutronData);
// 	  pmanager->AddDiscreteProcess(theNeutronElasticProcess);
	  
// 	  // inelastic scattering
// 	  G4NeutronInelasticProcess* theInelasticProcess =
// 	    new G4NeutronInelasticProcess();
	  
// 	  //JEC 12/2/07 USE ONE SPECIFIC hadronic model... 
// // 	  if (gheishahad) {
// 	    G4LENeutronInelastic* theInelasticModel = new G4LENeutronInelastic;
// 	    theInelasticModel->SetMinEnergy(19*MeV);
// 	    theInelasticProcess->RegisterMe(theInelasticModel);
// // 	  }
// // 	  else if (bertinihad) {
// // 	    G4CascadeInterface* theBertiniModel = new G4CascadeInterface;
// // 	    theBertiniModel->SetMinEnergy(19*MeV);
// // 	    theInelasticProcess->RegisterMe(theBertiniModel);
// // 	  }
// // 	  else if (binaryhad) {
// // 	    G4BinaryCascade* theBinaryModel = new G4BinaryCascade();
// // 	    theBinaryModel->SetMinEnergy(19*MeV);
// // 	    theInelasticProcess->RegisterMe(theBinaryModel);
// // 	  }
// // 	  else {
// // 	    G4cout << "No secondary interaction model chosen! Using G4 BINARY." << G4endl;
// // 	    G4BinaryCascade* theBinaryModel = new G4BinaryCascade();
// // 	    theBinaryModel->SetMinEnergy(19*MeV);
// // 	    theInelasticProcess->RegisterMe(theBinaryModel);
// // 	  }
	  
// 	  G4HENeutronInelastic* theHEInelasticModel = new G4HENeutronInelastic;
// 	  theInelasticProcess->RegisterMe(theHEInelasticModel);
	  
// 	  G4NeutronHPInelastic * theLENeutronInelasticModel =
// 	    new G4NeutronHPInelastic;
// 	  theInelasticProcess->RegisterMe(theLENeutronInelasticModel);
	  
// 	  G4NeutronHPInelasticData * theNeutronData1 = 
// 	    new G4NeutronHPInelasticData;
// 	  theInelasticProcess->AddDataSet(theNeutronData1);
// 	  pmanager->AddDiscreteProcess(theInelasticProcess);

// 	  // capture
// 	  G4HadronCaptureProcess* theCaptureProcess =
// 	    new G4HadronCaptureProcess;
// 	  G4LCapture* theCaptureModel = new G4LCapture;
// 	  theCaptureModel->SetMinEnergy(19*MeV);
// 	  theCaptureProcess->RegisterMe(theCaptureModel);
// 	  G4NeutronHPCapture * theLENeutronCaptureModel = new G4NeutronHPCapture;
// 	  theCaptureProcess->RegisterMe(theLENeutronCaptureModel);
// 	  G4NeutronHPCaptureData * theNeutronData3 = new G4NeutronHPCaptureData;
// 	  theCaptureProcess->AddDataSet(theNeutronData3);
// 	  pmanager->AddDiscreteProcess(theCaptureProcess);
// 	  //  G4ProcessManager* pmanager = G4Neutron::Neutron->GetProcessManager();
// 	  //  pmanager->AddProcess(new G4UserSpecialCuts(),-1,-1,1);
// 	}

//       else if (particleName == "anti_neutron") 
// 	{
// 	  pmanager->AddDiscreteProcess(theElasticProcess);
// 	  G4AntiNeutronInelasticProcess* theInelasticProcess = 
// 	    new G4AntiNeutronInelasticProcess();
// 	  G4LEAntiNeutronInelastic* theLEInelasticModel = 
// 	    new G4LEAntiNeutronInelastic;
// 	  theInelasticProcess->RegisterMe(theLEInelasticModel);
// 	  G4HEAntiNeutronInelastic* theHEInelasticModel = 
// 	    new G4HEAntiNeutronInelastic;
// 	  theInelasticProcess->RegisterMe(theHEInelasticModel);
// 	  pmanager->AddDiscreteProcess(theInelasticProcess);
// 	}

//       else if (particleName == "deuteron") 
// 	{
// 	  pmanager->AddDiscreteProcess(theElasticProcess);
// 	  G4DeuteronInelasticProcess* theInelasticProcess = 
// 	    new G4DeuteronInelasticProcess();
// 	  G4LEDeuteronInelastic* theLEInelasticModel = 
// 	    new G4LEDeuteronInelastic;
// 	  theInelasticProcess->RegisterMe(theLEInelasticModel);
// 	  pmanager->AddDiscreteProcess(theInelasticProcess);
// 	}
      
//       else if (particleName == "triton") 
// 	{
// 	  pmanager->AddDiscreteProcess(theElasticProcess);
// 	  G4TritonInelasticProcess* theInelasticProcess = 
// 	    new G4TritonInelasticProcess();
// 	  G4LETritonInelastic* theLEInelasticModel = 
// 	    new G4LETritonInelastic;
// 	  theInelasticProcess->RegisterMe(theLEInelasticModel);
// 	  pmanager->AddDiscreteProcess(theInelasticProcess);
// 	}

//       else if (particleName == "alpha") 
// 	{
// 	  pmanager->AddDiscreteProcess(theElasticProcess);
// 	  G4AlphaInelasticProcess* theInelasticProcess = 
// 	    new G4AlphaInelasticProcess();
// 	  G4LEAlphaInelastic* theLEInelasticModel = 
// 	    new G4LEAlphaInelastic;
// 	  theInelasticProcess->RegisterMe(theLEInelasticModel);
// 	  pmanager->AddDiscreteProcess(theInelasticProcess);
// 	}

//     }
// }//ConstructHad