8 #include "G4ParticleDefinition.hh"
14 #include "G4DecayTable.hh"
15 #include "G4ParticleTable.hh"
16 #include "G4PhaseSpaceDecayChannel.hh"
28 static const std::set<G4ParticleDefinition *>
particles =
load();
41 std::map< std::string, std::vector<std::vector<std::string>* > > decayMap;
45 std::istringstream is(
line);
46 std::vector<std::string>* txtvec =
new std::vector<std::string>(std::istream_iterator<std::string>(is),std::istream_iterator<std::string>());
47 const std::string
name = (*txtvec)[0];
49 decayMap.at(
name).push_back(txtvec);
51 catch (const::std::out_of_range&
e) {
52 std::vector<std::vector<std::string>* >
decays;
60 G4ParticleTable *theParticleTable = G4ParticleTable::GetParticleTable();
62 std::set<G4ParticleDefinition *>
particles;
65 G4String pType=
"custom";
67 G4double spectatormass = 0;
68 G4ParticleDefinition* spectator = 0;
73 G4double lifetime{-1.0};
79 G4cout <<
"-------------------------------------------------" << G4endl;
80 std::string::size_type beg_idx,end_idx;
82 beg_idx =
line.find_first_not_of(
"\t #");
83 if (beg_idx > 0 &&
line[beg_idx-1] ==
'#')
continue;
84 end_idx =
line.find_first_of(
"\t ", beg_idx);
85 if (end_idx == std::string::npos)
continue;
87 const std::string pdgCodeString =
line.substr( beg_idx, end_idx - beg_idx );
88 pdgCode = strtol(pdgCodeString.c_str(), &endptr, 0);
89 if (endptr[0] !=
'\0') {
90 throw std::invalid_argument(
"Could not convert string to int: " +
line.substr( beg_idx, end_idx - beg_idx ));
93 beg_idx =
line.find_first_not_of(
"\t ",end_idx);
94 end_idx =
line.find_first_of(
"\t #", beg_idx);
95 if (end_idx == std::string::npos)
continue;
96 mass =
atof(
line.substr( beg_idx, end_idx - beg_idx ).c_str());
98 beg_idx =
line.find_first_not_of(
"\t# ",end_idx);
99 end_idx =
line.length();
100 if (beg_idx==std::string::npos) beg_idx = end_idx;
101 name =
line.substr( beg_idx, end_idx - beg_idx );
102 while (
name.c_str()[0] ==
' ') {
name.erase(0,1); }
105 if (abs(pdgCode) / 1000000 == 0){
106 G4cout <<
"Pdg code too low " << pdgCode <<
" "<<abs(pdgCode) / 1000000 << std::endl;
117 G4cout<<
"pType: "<<pType<<G4endl;
120 G4ParticleDefinition* previousDefinition = theParticleTable->FindParticle(pdgCode);
122 if (previousDefinition) {
123 G4cout <<
"Found an existing G4ParticleDefinition for " <<
name <<
"(" << pdgCode <<
") in G4ParticleTable. Assuming that we want the new definition, so removing previous definition! May cause issues elsewhere though..." << G4endl;
124 previousDefinition->DumpTable();
125 theParticleTable->Remove(previousDefinition);
126 delete previousDefinition;
139 pType, 0, +1, pdgCode,
140 stable, lifetime,
nullptr );
141 if (pType==
"custom") {
142 spectatormass =
mass;
147 spectatormass =
mass;
150 G4String cloudname =
name+
"cloud";
151 G4String cloudtype = pType+
"cloud";
152 G4double cloudmass =
mass-spectatormass;
153 const bool cloudstable = !(
parameters->DoDecays() == 1 || decayMap.contains(cloudname));
154 const G4double cloudlifetime = cloudstable ? -1.0 :
parameters->Lifetime();
156 std::unique_ptr<G4ParticleDefinition> tmpParticle = std::unique_ptr<CustomParticle>(
new CustomParticle(
161 cloudstable, cloudlifetime,
nullptr ) );
165 G4cout <<
name <<
" being assigned cloud " <<
particle->GetCloud()->GetParticleName() <<
" and spectator " <<
particle->GetSpectator()->GetParticleName() << G4endl;
176 G4cout <<
"-------------------------------------------------" << G4endl;
181 if ( decayMap.contains(
name) ) {
182 std::vector<std::vector<std::string>* >& mydecays = decayMap.at(
name);
184 if (mydecays.empty()) {
continue; }
185 int ndec=mydecays.size();
186 G4DecayTable*
table =
new G4DecayTable();
187 G4VDecayChannel**
mode =
new G4VDecayChannel*[ndec]{};
188 for (
int i=0;
i!=ndec;
i++) {
189 std::vector<std::string> thisdec=*(mydecays[
i]);
190 if (thisdec.empty())
continue;
191 const double branch = std::stod(thisdec.back());
193 for (
const auto & thisString:thisdec) G4cout<<thisString<<G4endl;
194 if (thisdec.size()==3) {
195 mode[
i] =
new G4PhaseSpaceDecayChannel(thisdec[0],
branch,2,thisdec[1],thisdec[2]);
196 }
else if (thisdec.size()==4) {
197 mode[
i] =
new G4PhaseSpaceDecayChannel(thisdec[0],
branch,3,thisdec[1],thisdec[2],thisdec[3]);
199 G4cout<<
"Decay chain invalid!!!"<<std::endl;