|
| 1 | +//--------------------------------------------------------------------------- |
| 2 | +// Author: Alberto Ribon |
| 3 | +// Date: April 2016 |
| 4 | +// |
| 5 | +// Hadron physics for the new physics list FTFP_BERT_ATL. |
| 6 | +// This is a modified version of the FTFP_BERT hadron physics for ATLAS. |
| 7 | +// The hadron physics of FTFP_BERT_ATL has the transition between Bertini |
| 8 | +// (BERT) intra-nuclear cascade model and Fritiof (FTF) string model in the |
| 9 | +// energy region [9, 12] GeV (instead of [4, 5] GeV as in FTFP_BERT). |
| 10 | +//--------------------------------------------------------------------------- |
| 11 | +// |
| 12 | +#include <iomanip> |
| 13 | + |
| 14 | +#include "SimG4Core/PhysicsLists/interface/CMSHadronPhysicsFTFP_BERT_ATL.h" |
| 15 | + |
| 16 | +#include "globals.hh" |
| 17 | +#include "G4ios.hh" |
| 18 | +#include "G4SystemOfUnits.hh" |
| 19 | +#include "G4ParticleDefinition.hh" |
| 20 | +#include "G4ParticleTable.hh" |
| 21 | + |
| 22 | +#include "G4MesonConstructor.hh" |
| 23 | +#include "G4BaryonConstructor.hh" |
| 24 | +#include "G4ShortLivedConstructor.hh" |
| 25 | + |
| 26 | +#include "G4ComponentGGHadronNucleusXsc.hh" |
| 27 | +#include "G4CrossSectionInelastic.hh" |
| 28 | +#include "G4HadronCaptureProcess.hh" |
| 29 | +#include "G4NeutronRadCapture.hh" |
| 30 | +#include "G4NeutronInelasticXS.hh" |
| 31 | +#include "G4NeutronCaptureXS.hh" |
| 32 | + |
| 33 | +#include "G4CrossSectionDataSetRegistry.hh" |
| 34 | + |
| 35 | +#include "G4PhysListUtil.hh" |
| 36 | + |
| 37 | +G4ThreadLocal CMSHadronPhysicsFTFP_BERT_ATL::ThreadPrivate* CMSHadronPhysicsFTFP_BERT_ATL::tpdata=nullptr; |
| 38 | + |
| 39 | +CMSHadronPhysicsFTFP_BERT_ATL::CMSHadronPhysicsFTFP_BERT_ATL(G4int) |
| 40 | + : G4VPhysicsConstructor("hInelastic FTFP_BERT_ATL") |
| 41 | + , QuasiElastic(false) |
| 42 | +{} |
| 43 | + |
| 44 | +void CMSHadronPhysicsFTFP_BERT_ATL::CreateModels() |
| 45 | +{ |
| 46 | + G4double minFTFP = 9.0 * GeV; |
| 47 | + G4double maxBERT = 12.0 * GeV; |
| 48 | + G4cout << " CMS_FTFP_BERT_ATL : new threshold between BERT and FTFP" |
| 49 | + << " is over the interval " << minFTFP/GeV << " to " << maxBERT/GeV |
| 50 | + << " GeV." << G4endl; |
| 51 | + QuasiElastic= false; |
| 52 | + |
| 53 | + tpdata->theNeutrons=new G4NeutronBuilder; |
| 54 | + tpdata->theNeutrons->RegisterMe(tpdata->theFTFPNeutron=new G4FTFPNeutronBuilder(QuasiElastic)); |
| 55 | + tpdata->theFTFPNeutron->SetMinEnergy(minFTFP); |
| 56 | + tpdata->theNeutrons->RegisterMe(tpdata->theBertiniNeutron=new G4BertiniNeutronBuilder); |
| 57 | + tpdata->theBertiniNeutron->SetMinEnergy(0.0*GeV); |
| 58 | + tpdata->theBertiniNeutron->SetMaxEnergy(maxBERT); |
| 59 | + |
| 60 | + tpdata->thePro=new G4ProtonBuilder; |
| 61 | + tpdata->thePro->RegisterMe(tpdata->theFTFPPro=new G4FTFPProtonBuilder(QuasiElastic)); |
| 62 | + tpdata->theFTFPPro->SetMinEnergy(minFTFP); |
| 63 | + tpdata->thePro->RegisterMe(tpdata->theBertiniPro=new G4BertiniProtonBuilder); |
| 64 | + tpdata->theBertiniPro->SetMaxEnergy(maxBERT); |
| 65 | + |
| 66 | + tpdata->thePiK=new G4PiKBuilder; |
| 67 | + tpdata->thePiK->RegisterMe(tpdata->theFTFPPiK=new G4FTFPPiKBuilder(QuasiElastic)); |
| 68 | + tpdata->theFTFPPiK->SetMinEnergy(minFTFP); |
| 69 | + tpdata->thePiK->RegisterMe(tpdata->theBertiniPiK=new G4BertiniPiKBuilder); |
| 70 | + tpdata->theBertiniPiK->SetMaxEnergy(maxBERT); |
| 71 | + |
| 72 | + tpdata->theHyperon=new G4HyperonFTFPBuilder; |
| 73 | + |
| 74 | + tpdata->theAntiBaryon=new G4AntiBarionBuilder; |
| 75 | + tpdata->theAntiBaryon->RegisterMe(tpdata->theFTFPAntiBaryon=new G4FTFPAntiBarionBuilder(QuasiElastic)); |
| 76 | +} |
| 77 | + |
| 78 | +CMSHadronPhysicsFTFP_BERT_ATL::~CMSHadronPhysicsFTFP_BERT_ATL() |
| 79 | +{ |
| 80 | + if (!tpdata) return; |
| 81 | + |
| 82 | + delete tpdata->theNeutrons; |
| 83 | + delete tpdata->theBertiniNeutron; |
| 84 | + delete tpdata->theFTFPNeutron; |
| 85 | + |
| 86 | + delete tpdata->thePiK; |
| 87 | + delete tpdata->theBertiniPiK; |
| 88 | + delete tpdata->theFTFPPiK; |
| 89 | + |
| 90 | + delete tpdata->thePro; |
| 91 | + delete tpdata->theBertiniPro; |
| 92 | + delete tpdata->theFTFPPro; |
| 93 | + |
| 94 | + delete tpdata->theHyperon; |
| 95 | + delete tpdata->theAntiBaryon; |
| 96 | + delete tpdata->theFTFPAntiBaryon; |
| 97 | + |
| 98 | + //Note that here we need to set to 0 the pointer |
| 99 | + //since tpdata is static and if thread are "reused" |
| 100 | + //it can be problematic |
| 101 | + delete tpdata; tpdata = nullptr; |
| 102 | +} |
| 103 | + |
| 104 | +void CMSHadronPhysicsFTFP_BERT_ATL::ConstructParticle() |
| 105 | +{ |
| 106 | + G4MesonConstructor pMesonConstructor; |
| 107 | + pMesonConstructor.ConstructParticle(); |
| 108 | + |
| 109 | + G4BaryonConstructor pBaryonConstructor; |
| 110 | + pBaryonConstructor.ConstructParticle(); |
| 111 | + |
| 112 | + G4ShortLivedConstructor pShortLivedConstructor; |
| 113 | + pShortLivedConstructor.ConstructParticle(); |
| 114 | +} |
| 115 | + |
| 116 | +#include "G4ProcessManager.hh" |
| 117 | +void CMSHadronPhysicsFTFP_BERT_ATL::ConstructProcess() |
| 118 | +{ |
| 119 | + if ( tpdata == 0 ) tpdata = new ThreadPrivate; |
| 120 | + CreateModels(); |
| 121 | + tpdata->theNeutrons->Build(); |
| 122 | + tpdata->thePro->Build(); |
| 123 | + tpdata->thePiK->Build(); |
| 124 | + |
| 125 | + // --- Kaons --- |
| 126 | + tpdata->xsKaon = new G4ComponentGGHadronNucleusXsc(); |
| 127 | + G4VCrossSectionDataSet * kaonxs = new G4CrossSectionInelastic(tpdata->xsKaon); |
| 128 | + G4PhysListUtil::FindInelasticProcess(G4KaonMinus::KaonMinus())->AddDataSet(kaonxs); |
| 129 | + G4PhysListUtil::FindInelasticProcess(G4KaonPlus::KaonPlus())->AddDataSet(kaonxs); |
| 130 | + G4PhysListUtil::FindInelasticProcess(G4KaonZeroShort::KaonZeroShort())->AddDataSet(kaonxs); |
| 131 | + G4PhysListUtil::FindInelasticProcess(G4KaonZeroLong::KaonZeroLong())->AddDataSet(kaonxs); |
| 132 | + |
| 133 | + tpdata->theHyperon->Build(); |
| 134 | + tpdata->theAntiBaryon->Build(); |
| 135 | + |
| 136 | + // --- Neutrons --- |
| 137 | + tpdata->xsNeutronInelasticXS = new G4NeutronInelasticXS(); |
| 138 | + |
| 139 | + G4HadronicProcess* capture = nullptr; |
| 140 | + G4ProcessManager* pmanager = G4Neutron::Neutron()->GetProcessManager(); |
| 141 | + G4ProcessVector* pv = pmanager->GetProcessList(); |
| 142 | + for ( size_t i=0; i < static_cast<size_t>(pv->size()); ++i ) { |
| 143 | + if ( fCapture == ((*pv)[i])->GetProcessSubType() ) { |
| 144 | + capture = static_cast<G4HadronicProcess*>((*pv)[i]); |
| 145 | + } |
| 146 | + } |
| 147 | + if ( ! capture ) { |
| 148 | + capture = new G4HadronCaptureProcess("nCapture"); |
| 149 | + pmanager->AddDiscreteProcess(capture); |
| 150 | + } |
| 151 | + tpdata->xsNeutronCaptureXS = new G4NeutronCaptureXS(); |
| 152 | + capture->AddDataSet(tpdata->xsNeutronCaptureXS); |
| 153 | + capture->RegisterMe(new G4NeutronRadCapture()); |
| 154 | +} |
0 commit comments