HadronPhysicsQGSPCMS_FTFP_BERT.cc

Go to the documentation of this file.

#include "SimG4Core/PhysicsLists/interface/HadronPhysicsQGSPCMS_FTFP_BERT.h"

#include "globals.hh"
#include "G4ios.hh"
#include <iomanip>   
#include "G4ParticleDefinition.hh"
#include "G4ParticleTable.hh"

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

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

#include "G4PhysListUtil.hh"
#include "G4SystemOfUnits.hh"

G4ThreadLocal HadronPhysicsQGSPCMS_FTFP_BERT::ThreadPrivate* 
HadronPhysicsQGSPCMS_FTFP_BERT::tpdata = nullptr;

HadronPhysicsQGSPCMS_FTFP_BERT::HadronPhysicsQGSPCMS_FTFP_BERT(G4int)
  :  G4VPhysicsConstructor("hInelasticQGSPCMS_FTFP_BERT")
{}

void HadronPhysicsQGSPCMS_FTFP_BERT::CreateModels()
{
  // First transition, between BERT and FTF/P
  G4double minFTFP= 6.0 * GeV;  
  G4double maxBERT= 8.0 * GeV;  
  // Second transition, between FTF/P and QGS/P
  G4double minQGSP= 12.0 * GeV;
  G4double maxFTFP= 25.0 * GeV; 

  G4bool   quasiElasFTF= false;   // Use built-in quasi-elastic (not add-on)
  G4bool   quasiElasQGS= true;    // For QGS, it must use it.

  G4cout << " New QGSPCMS_FTFP_BERT hadronic inealstic physics" << G4endl;
  G4cout << "   Thresholds: " << G4endl;
  G4cout << "     1) between BERT  and FTFP over the interval " 
     << minFTFP/GeV << " to " << maxBERT/GeV << " GeV. " << G4endl;
  G4cout << "     2) between FTFP and QGSP over the interval " 
     << minQGSP/GeV << " to " << maxFTFP/GeV << " GeV. " << G4endl;
  G4cout << "   QuasiElastic: " << quasiElasQGS << " for QGS "
     << " and " << quasiElasFTF << " for FTF " << G4endl;

  tpdata->theNeutrons=new G4NeutronBuilder;
  tpdata->theNeutrons->RegisterMe(
    tpdata->theQGSPNeutron=new G4QGSPNeutronBuilder(quasiElasQGS));
  tpdata->theQGSPNeutron->SetMinEnergy(minQGSP);   
  tpdata->theNeutrons->RegisterMe(
    tpdata->theFTFPNeutron=new G4FTFPNeutronBuilder(quasiElasFTF));
  tpdata->theFTFPNeutron->SetMinEnergy(minFTFP);   
  tpdata->theFTFPNeutron->SetMaxEnergy(maxFTFP);   

  tpdata->theNeutrons->RegisterMe(
    tpdata->theBertiniNeutron=new G4BertiniNeutronBuilder);
  tpdata->theBertiniNeutron->SetMinEnergy(0.0*GeV);
  tpdata->theBertiniNeutron->SetMaxEnergy(maxBERT);  

  tpdata->thePro=new G4ProtonBuilder;
  tpdata->thePro->RegisterMe(
    tpdata->theQGSPPro=new G4QGSPProtonBuilder(quasiElasQGS));
  tpdata->theQGSPPro->SetMinEnergy(minQGSP);   
  tpdata->thePro->RegisterMe(
    tpdata->theFTFPPro=new G4FTFPProtonBuilder(quasiElasFTF));
  tpdata->theFTFPPro->SetMinEnergy(minFTFP);   
  tpdata->theFTFPPro->SetMaxEnergy(maxFTFP);   
  tpdata->thePro->RegisterMe(
    tpdata->theBertiniPro=new G4BertiniProtonBuilder);
  tpdata->theBertiniPro->SetMaxEnergy(maxBERT); 
  
  tpdata->thePiK=new G4PiKBuilder;
  tpdata->thePiK->RegisterMe(
    tpdata->theQGSPPiK=new G4QGSPPiKBuilder(quasiElasQGS));
  tpdata->theQGSPPiK->SetMinEnergy(minQGSP);   
  tpdata->thePiK->RegisterMe(
    tpdata->theFTFPPiK=new G4FTFPPiKBuilder(quasiElasFTF));
  tpdata->theFTFPPiK->SetMaxEnergy(maxFTFP);  
  tpdata->theFTFPPiK->SetMinEnergy(minFTFP); 
  tpdata->thePiK->RegisterMe(tpdata->theBertiniPiK=new G4BertiniPiKBuilder);
  tpdata->theBertiniPiK->SetMaxEnergy(maxBERT); 
  
  // Hyperons use FTF
  tpdata->theHyperon=new G4HyperonFTFPBuilder;

  tpdata->theAntiBaryon=new G4AntiBarionBuilder;
  tpdata->theAntiBaryon->RegisterMe(
    tpdata->theFTFPAntiBaryon=new G4FTFPAntiBarionBuilder(quasiElasFTF));
}

HadronPhysicsQGSPCMS_FTFP_BERT::~HadronPhysicsQGSPCMS_FTFP_BERT()
{
  if(nullptr != tpdata) {
   delete tpdata->theQGSPNeutron;
   delete tpdata->theFTFPNeutron;
   delete tpdata->theBertiniNeutron;
   delete tpdata->theNeutrons;

   delete tpdata->theQGSPPro;
   delete tpdata->theFTFPPro;
   delete tpdata->thePro;
   delete tpdata->theBertiniPro;

   delete tpdata->theQGSPPiK;
   delete tpdata->theFTFPPiK;
   delete tpdata->theBertiniPiK;
   delete tpdata->thePiK;

   delete tpdata->theHyperon;
   delete tpdata->theAntiBaryon;
   delete tpdata->theFTFPAntiBaryon;
   
   delete tpdata;
   tpdata = nullptr;
  }
}

void HadronPhysicsQGSPCMS_FTFP_BERT::ConstructParticle()
{
  G4MesonConstructor pMesonConstructor;
  pMesonConstructor.ConstructParticle();

  G4BaryonConstructor pBaryonConstructor;
  pBaryonConstructor.ConstructParticle();

  G4ShortLivedConstructor pShortLivedConstructor;
  pShortLivedConstructor.ConstructParticle();  

  G4IonConstructor pIonConstructor;
  pIonConstructor.ConstructParticle();
}

#include "G4ProcessManager.hh"
void HadronPhysicsQGSPCMS_FTFP_BERT::ConstructProcess()
{
  if ( tpdata == nullptr ) { tpdata = new ThreadPrivate; }
  CreateModels();
  tpdata->theNeutrons->Build();
  tpdata->thePro->Build();
  tpdata->thePiK->Build();
  tpdata->theHyperon->Build(); 
  tpdata->theAntiBaryon->Build(); 

  // --- Neutrons ---
  G4PhysListUtil::FindInelasticProcess(G4Neutron::Neutron())
    ->AddDataSet(new G4NeutronInelasticXS());

  G4HadronicProcess* capture = nullptr;
  G4ProcessManager* pmanager = G4Neutron::Neutron()->GetProcessManager();
  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);
  }
  capture->AddDataSet(new G4NeutronCaptureXS());
  capture->RegisterMe(new G4NeutronRadCapture());
}