| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
>> 26 //
26 // Hadronic Process: Nuclear De-excitations 27 // Hadronic Process: Nuclear De-excitations
27 // by V. Lara (May 1998) 28 // by V. Lara (May 1998)
28 // 29 //
29 // 30 //
30 // Modified: 31 // Modified:
31 // 30 June 1998 by V. Lara: 32 // 30 June 1998 by V. Lara:
32 // -Modified the Transform method for use 33 // -Modified the Transform method for use G4ParticleTable and
33 // therefore G4IonTable. It makes possib 34 // therefore G4IonTable. It makes possible to convert all kind
34 // of fragments (G4Fragment) produced in 35 // of fragments (G4Fragment) produced in deexcitation to
35 // G4DynamicParticle 36 // G4DynamicParticle
36 // -It uses default algorithms for: 37 // -It uses default algorithms for:
37 // Evaporation: G4Evaporation 38 // Evaporation: G4Evaporation
38 // MultiFragmentation: G4StatMF 39 // MultiFragmentation: G4StatMF
39 // Fermi Breakup model: G4FermiBr 40 // Fermi Breakup model: G4FermiBreakUp
40 // 24 Jul 2008 by M. A. Cortes Giraldo: 41 // 24 Jul 2008 by M. A. Cortes Giraldo:
41 // -Max Z,A for Fermi Break-Up turns to 9 42 // -Max Z,A for Fermi Break-Up turns to 9,17 by default
42 // -BreakItUp() reorganised and bug in Ev 43 // -BreakItUp() reorganised and bug in Evaporation loop fixed
43 // -Transform() optimised 44 // -Transform() optimised
44 // (September 2008) by J. M. Quesada. External 45 // (September 2008) by J. M. Quesada. External choices have been added for :
45 // -inverse cross section option (default 46 // -inverse cross section option (default OPTxs=3)
46 // -superimposed Coulomb barrier (if useS 47 // -superimposed Coulomb barrier (if useSICB is set true, by default it is false)
47 // September 2009 by J. M. Quesada: 48 // September 2009 by J. M. Quesada:
48 // -according to Igor Pshenichnov, SMM wi 49 // -according to Igor Pshenichnov, SMM will be applied (just in case) only once.
49 // 27 Nov 2009 by V.Ivanchenko: 50 // 27 Nov 2009 by V.Ivanchenko:
50 // -cleanup the logic, reduce number inte 51 // -cleanup the logic, reduce number internal vectors, fixed memory leak.
51 // 11 May 2010 by V.Ivanchenko: 52 // 11 May 2010 by V.Ivanchenko:
52 // -FermiBreakUp activated, used integer 53 // -FermiBreakUp activated, used integer Z and A, used BreakUpFragment method for
53 // final photon deexcitation; used check 54 // final photon deexcitation; used check on adundance of a fragment, decay
54 // unstable fragments with A <5 55 // unstable fragments with A <5
55 // 22 March 2011 by V.Ivanchenko: general clea 56 // 22 March 2011 by V.Ivanchenko: general cleanup and addition of a condition:
56 // products of Fermi Break Up cannot be 57 // products of Fermi Break Up cannot be further deexcited by this model
57 // 30 March 2011 by V.Ivanchenko removed priva 58 // 30 March 2011 by V.Ivanchenko removed private inline methods, moved Set methods
58 // to the source 59 // to the source
59 // 23 January 2012 by V.Ivanchenko general cle 60 // 23 January 2012 by V.Ivanchenko general cleanup including destruction of
60 // objects, propagate G4PhotonEvaporation p 61 // objects, propagate G4PhotonEvaporation pointer to G4Evaporation class and
61 // not delete it here 62 // not delete it here
62 63
63 #include "G4ExcitationHandler.hh" 64 #include "G4ExcitationHandler.hh"
64 #include "G4SystemOfUnits.hh" 65 #include "G4SystemOfUnits.hh"
65 #include "G4LorentzVector.hh" 66 #include "G4LorentzVector.hh"
66 #include "G4ThreeVector.hh" << 67 #include "G4NistManager.hh"
67 #include "G4ParticleTable.hh" 68 #include "G4ParticleTable.hh"
68 #include "G4ParticleTypes.hh" 69 #include "G4ParticleTypes.hh"
69 #include "G4Ions.hh" 70 #include "G4Ions.hh"
70 #include "G4Electron.hh" 71 #include "G4Electron.hh"
71 #include "G4Lambda.hh" <<
72 72
73 #include "G4VMultiFragmentation.hh" 73 #include "G4VMultiFragmentation.hh"
74 #include "G4VFermiBreakUp.hh" 74 #include "G4VFermiBreakUp.hh"
75 #include "G4Element.hh" <<
76 #include "G4ElementTable.hh" <<
77 75
78 #include "G4VEvaporation.hh" 76 #include "G4VEvaporation.hh"
79 #include "G4VEvaporationChannel.hh" 77 #include "G4VEvaporationChannel.hh"
80 #include "G4Evaporation.hh" 78 #include "G4Evaporation.hh"
81 #include "G4PhotonEvaporation.hh" 79 #include "G4PhotonEvaporation.hh"
82 #include "G4StatMF.hh" 80 #include "G4StatMF.hh"
83 #include "G4FermiBreakUpVI.hh" 81 #include "G4FermiBreakUpVI.hh"
84 #include "G4NuclearLevelData.hh" 82 #include "G4NuclearLevelData.hh"
85 #include "G4PhysicsModelCatalog.hh" << 83 #include "G4Pow.hh"
86 84
87 G4ExcitationHandler::G4ExcitationHandler() 85 G4ExcitationHandler::G4ExcitationHandler()
88 : minEForMultiFrag(1.*CLHEP::TeV), minExcita << 86 : maxZForFermiBreakUp(9),maxAForFermiBreakUp(17),
89 maxExcitation(100.*CLHEP::MeV) << 87 fVerbose(0),isInitialised(false),isEvapLocal(true)
90 { << 88 {
91 thePartTable = G4ParticleTable::GetParticleT << 89 theTableOfIons = G4ParticleTable::GetParticleTable()->GetIonTable();
92 theTableOfIons = thePartTable->GetIonTable() <<
93 nist = G4NistManager::Instance(); 90 nist = G4NistManager::Instance();
94 91
95 theMultiFragmentation = new G4StatMF(); << 92 theMultiFragmentation = nullptr;
96 theFermiModel = new G4FermiBreakUpVI(); << 93 theFermiModel = nullptr;
97 thePhotonEvaporation = new G4PhotonEvaporati << 94 G4Pow::GetInstance();
98 SetEvaporation(new G4Evaporation(thePhotonEv << 95 theEvaporation = new G4Evaporation();
>> 96 thePhotonEvaporation = theEvaporation->GetPhotonEvaporation();
99 theResults.reserve(60); 97 theResults.reserve(60);
100 results.reserve(30); 98 results.reserve(30);
101 theEvapList.reserve(30); 99 theEvapList.reserve(30);
>> 100 thePhotoEvapList.reserve(10);
102 101
103 theElectron = G4Electron::Electron(); << 102 SetParameters();
104 theNeutron = G4Neutron::NeutronDefinition(); << 103 electron = G4Electron::Electron();
105 theProton = G4Proton::ProtonDefinition(); << 104
106 theDeuteron = G4Deuteron::DeuteronDefinition << 105 if(fVerbose > 0) { G4cout << "### New handler " << this << G4endl; }
107 theTriton = G4Triton::TritonDefinition(); <<
108 theHe3 = G4He3::He3Definition(); <<
109 theAlpha = G4Alpha::AlphaDefinition(); <<
110 theLambda = G4Lambda::Lambda(); <<
111 <<
112 fLambdaMass = theLambda->GetPDGMass(); <<
113 <<
114 if(fVerbose > 1) { G4cout << "### New handle <<
115 } 106 }
116 107
117 G4ExcitationHandler::~G4ExcitationHandler() 108 G4ExcitationHandler::~G4ExcitationHandler()
118 { 109 {
>> 110 //G4cout << "### Delete handler " << this << G4endl;
119 delete theMultiFragmentation; 111 delete theMultiFragmentation;
120 delete theFermiModel; 112 delete theFermiModel;
121 if(isEvapLocal) { delete theEvaporation; } 113 if(isEvapLocal) { delete theEvaporation; }
122 } 114 }
123 115
124 void G4ExcitationHandler::SetParameters() 116 void G4ExcitationHandler::SetParameters()
125 { 117 {
126 G4NuclearLevelData* ndata = G4NuclearLevelDa << 118 G4DeexPrecoParameters* param =
127 auto param = ndata->GetParameters(); << 119 G4NuclearLevelData::GetInstance()->GetParameters();
128 isActive = true; 120 isActive = true;
129 // check if de-excitation is needed << 121 if(fDummy == param->GetDeexChannelsType()) { isActive = false; }
130 if (fDummy == param->GetDeexChannelsType()) <<
131 isActive = false; <<
132 } else { <<
133 // upload data for elements used in geomet <<
134 G4int Zmax = 20; <<
135 const G4ElementTable* table = G4Element::G <<
136 for (auto const & elm : *table) { Zmax = s <<
137 ndata->UploadNuclearLevelData(Zmax+1); <<
138 } <<
139 minEForMultiFrag = param->GetMinExPerNucleou 122 minEForMultiFrag = param->GetMinExPerNucleounForMF();
140 minExcitation = param->GetMinExcitation(); 123 minExcitation = param->GetMinExcitation();
141 maxExcitation = param->GetPrecoHighEnergy(); << 124 icID = param->GetInternalConversionID();
142 << 125 if(isActive) {
143 // allowing local debug printout << 126 if(!thePhotonEvaporation) {
144 fVerbose = std::max(fVerbose, param->GetVerb << 127 SetPhotonEvaporation(new G4PhotonEvaporation());
145 if (isActive) { << 128 }
146 if (nullptr == thePhotonEvaporation) { << 129 if(!theFermiModel) { SetFermiModel(new G4FermiBreakUpVI()); }
147 SetPhotonEvaporation(new G4PhotonEvapora << 130 if(!theMultiFragmentation) { SetMultiFragmentation(new G4StatMF()); }
148 } <<
149 if (nullptr == theFermiModel) { <<
150 SetFermiModel(new G4FermiBreakUpVI()); <<
151 } <<
152 if (nullptr == theMultiFragmentation) { <<
153 SetMultiFragmentation(new G4StatMF()); <<
154 } <<
155 if (nullptr == theEvaporation) { <<
156 SetEvaporation(new G4Evaporation(thePhot <<
157 } <<
158 } <<
159 theFermiModel->SetVerbose(fVerbose); <<
160 if(fVerbose > 1) { <<
161 G4cout << "G4ExcitationHandler::SetParamet <<
162 } 131 }
163 } 132 }
164 133
165 void G4ExcitationHandler::Initialise() 134 void G4ExcitationHandler::Initialise()
166 { 135 {
167 if(isInitialised) { return; } 136 if(isInitialised) { return; }
168 if(fVerbose > 1) { << 137 if(fVerbose > 0) {
169 G4cout << "G4ExcitationHandler::Initialise 138 G4cout << "G4ExcitationHandler::Initialise() started " << this << G4endl;
170 } 139 }
171 G4DeexPrecoParameters* param = 140 G4DeexPrecoParameters* param =
172 G4NuclearLevelData::GetInstance()->GetPara 141 G4NuclearLevelData::GetInstance()->GetParameters();
173 isInitialised = true; 142 isInitialised = true;
174 SetParameters(); 143 SetParameters();
175 if(isActive) { 144 if(isActive) {
176 theFermiModel->Initialise(); 145 theFermiModel->Initialise();
177 theEvaporation->InitialiseChannels(); 146 theEvaporation->InitialiseChannels();
178 } 147 }
179 // dump level is controlled by parameter cla << 148 if(fVerbose > 0) { param->Dump(); }
180 param->Dump(); <<
181 } 149 }
182 150
183 void G4ExcitationHandler::SetEvaporation(G4VEv 151 void G4ExcitationHandler::SetEvaporation(G4VEvaporation* ptr, G4bool flag)
184 { 152 {
185 if(nullptr != ptr && ptr != theEvaporation) << 153 if(ptr && ptr != theEvaporation) {
>> 154 delete theEvaporation;
186 theEvaporation = ptr; 155 theEvaporation = ptr;
187 theEvaporation->SetPhotonEvaporation(thePh << 156 thePhotonEvaporation = ptr->GetPhotonEvaporation();
188 theEvaporation->SetFermiBreakUp(theFermiMo 157 theEvaporation->SetFermiBreakUp(theFermiModel);
189 isEvapLocal = flag; 158 isEvapLocal = flag;
190 if(fVerbose > 1) { <<
191 G4cout << "G4ExcitationHandler::SetEvapo <<
192 } <<
193 } 159 }
194 } 160 }
195 161
196 void 162 void
197 G4ExcitationHandler::SetMultiFragmentation(G4V 163 G4ExcitationHandler::SetMultiFragmentation(G4VMultiFragmentation* ptr)
198 { 164 {
199 if(nullptr != ptr && ptr != theMultiFragment << 165 if(ptr && ptr != theMultiFragmentation) {
200 delete theMultiFragmentation; 166 delete theMultiFragmentation;
201 theMultiFragmentation = ptr; 167 theMultiFragmentation = ptr;
202 } 168 }
203 } 169 }
204 170
205 void G4ExcitationHandler::SetFermiModel(G4VFer 171 void G4ExcitationHandler::SetFermiModel(G4VFermiBreakUp* ptr)
206 { 172 {
207 if(nullptr != ptr && ptr != theFermiModel) { << 173 if(ptr && ptr != theFermiModel) {
208 delete theFermiModel; 174 delete theFermiModel;
209 theFermiModel = ptr; 175 theFermiModel = ptr;
210 if(nullptr != theEvaporation) { << 176 theEvaporation->SetFermiBreakUp(theFermiModel);
211 theEvaporation->SetFermiBreakUp(theFermi <<
212 } <<
213 } 177 }
214 } 178 }
215 179
216 void 180 void
217 G4ExcitationHandler::SetPhotonEvaporation(G4VE 181 G4ExcitationHandler::SetPhotonEvaporation(G4VEvaporationChannel* ptr)
218 { 182 {
219 if(nullptr != ptr && ptr != thePhotonEvapora << 183 if(ptr && ptr != thePhotonEvaporation) {
220 delete thePhotonEvaporation; <<
221 thePhotonEvaporation = ptr; 184 thePhotonEvaporation = ptr;
222 if(nullptr != theEvaporation) { << 185 theEvaporation->SetPhotonEvaporation(ptr);
223 theEvaporation->SetPhotonEvaporation(ptr <<
224 } <<
225 if(fVerbose > 1) { <<
226 G4cout << "G4ExcitationHandler::SetPhoto <<
227 << " for handler " << this << G4e <<
228 } <<
229 } 186 }
230 } 187 }
231 188
232 void G4ExcitationHandler::SetDeexChannelsType( 189 void G4ExcitationHandler::SetDeexChannelsType(G4DeexChannelType val)
233 { 190 {
234 G4Evaporation* evap = static_cast<G4Evaporat 191 G4Evaporation* evap = static_cast<G4Evaporation*>(theEvaporation);
235 if(fVerbose > 1) { <<
236 G4cout << "G4ExcitationHandler::SetDeexCha <<
237 << " for " << this << G4endl; <<
238 } <<
239 if(val == fDummy) { 192 if(val == fDummy) {
240 isActive = false; 193 isActive = false;
241 return; 194 return;
242 } 195 }
243 if (nullptr == evap) { return; } << 196 if(!evap) { return; }
244 if (val == fEvaporation) { << 197 if(val == fEvaporation) {
245 evap->SetDefaultChannel(); 198 evap->SetDefaultChannel();
246 } else if (val == fCombined) { << 199 } else if(val == fCombined) {
247 evap->SetCombinedChannel(); 200 evap->SetCombinedChannel();
248 } else if (val == fGEM) { << 201 } else if(val == fGEM) {
249 evap->SetGEMChannel(); 202 evap->SetGEMChannel();
250 } else if (val == fGEMVI) { <<
251 evap->SetGEMVIChannel(); <<
252 } 203 }
253 evap->InitialiseChannels(); 204 evap->InitialiseChannels();
254 if (fVerbose > 1) { << 205 if(fVerbose > 0) {
255 if (G4Threading::IsMasterThread()) { << 206 if(G4Threading::IsMasterThread()) {
256 G4cout << "Number of de-excitation chann 207 G4cout << "Number of de-excitation channels is changed to: "
257 << theEvaporation->GetNumberOfChannels( 208 << theEvaporation->GetNumberOfChannels();
258 G4cout << " " << this; 209 G4cout << " " << this;
259 } 210 }
260 G4cout << G4endl; 211 G4cout << G4endl;
261 } 212 }
262 } 213 }
263 214
264 G4VEvaporation* G4ExcitationHandler::GetEvapor <<
265 { <<
266 if (nullptr != theEvaporation) { SetParamete <<
267 return theEvaporation; <<
268 } <<
269 <<
270 G4VMultiFragmentation* G4ExcitationHandler::Ge <<
271 { <<
272 if (nullptr != theMultiFragmentation) { SetP <<
273 return theMultiFragmentation; <<
274 } <<
275 <<
276 G4VFermiBreakUp* G4ExcitationHandler::GetFermi <<
277 { <<
278 if (nullptr != theFermiModel) { SetParameter <<
279 return theFermiModel; <<
280 } <<
281 <<
282 G4VEvaporationChannel* G4ExcitationHandler::Ge <<
283 { <<
284 if(nullptr != thePhotonEvaporation) { SetPar <<
285 return thePhotonEvaporation; <<
286 } <<
287 <<
288 G4ReactionProductVector * 215 G4ReactionProductVector *
289 G4ExcitationHandler::BreakItUp(const G4Fragmen 216 G4ExcitationHandler::BreakItUp(const G4Fragment & theInitialState)
290 { 217 {
291 // Variables existing until end of method 218 // Variables existing until end of method
292 G4Fragment * theInitialStatePtr = new G4Frag 219 G4Fragment * theInitialStatePtr = new G4Fragment(theInitialState);
293 if (fVerbose > 1) { << 220 if(fVerbose > 1) {
294 G4cout << "@@@@@@@@@@ Start G4Excitation H 221 G4cout << "@@@@@@@@@@ Start G4Excitation Handler @@@@@@@@@@@@@ " << G4endl;
295 G4cout << theInitialState << G4endl; 222 G4cout << theInitialState << G4endl;
296 } 223 }
297 if (!isInitialised) { Initialise(); } << 224 if(!isInitialised) { Initialise(); }
298 225
299 // pointer to fragment vector which receives 226 // pointer to fragment vector which receives temporal results
300 G4FragmentVector * theTempResult = nullptr; 227 G4FragmentVector * theTempResult = nullptr;
301 228
302 theResults.clear(); 229 theResults.clear();
>> 230 thePhotoEvapList.clear();
303 theEvapList.clear(); 231 theEvapList.clear();
304 232
305 // Variables to describe the excited configu 233 // Variables to describe the excited configuration
306 G4double exEnergy = theInitialState.GetExcit 234 G4double exEnergy = theInitialState.GetExcitationEnergy();
307 G4int A = theInitialState.GetA_asInt(); 235 G4int A = theInitialState.GetA_asInt();
308 G4int Z = theInitialState.GetZ_asInt(); 236 G4int Z = theInitialState.GetZ_asInt();
309 G4int nL = theInitialState.GetNumberOfLambda << 237
310 <<
311 // too much excitation <<
312 if (exEnergy > A*maxExcitation && A > 0) { <<
313 ++fWarnings; <<
314 if(fWarnings < 0) { <<
315 G4ExceptionDescription ed; <<
316 ed << "High excitation Fragment Z= " << <<
317 << " Eex/A(MeV)= " << exEnergy/A; <<
318 G4Exception("G4ExcitationHandler::BreakI <<
319 } <<
320 } <<
321 <<
322 // for hyper-nuclei subtract lambdas from th <<
323 G4double lambdaF = 0.0; <<
324 G4LorentzVector lambdaLV = theInitialStatePt <<
325 if (0 < nL) { <<
326 <<
327 // is it a stable hyper-nuclei? <<
328 if(A >= 3 && A <= 5 && nL <= 2) { <<
329 G4int pdg = 0; <<
330 if(3 == A && 1 == nL) { <<
331 pdg = 1010010030; <<
332 } else if(5 == A && 2 == Z && 1 == nL) { <<
333 pdg = 1010020050; <<
334 } else if(4 == A) { <<
335 if(1 == Z && 1 == nL) { <<
336 pdg = 1010010040; <<
337 } else if(2 == Z && 1 == nL) { <<
338 pdg = 1010020040; <<
339 } else if(0 == Z && 2 == nL) { <<
340 pdg = 1020000040; <<
341 } else if(1 == Z && 2 == nL) { <<
342 pdg = 1020010040; <<
343 } <<
344 } <<
345 // initial state is one of hyper-nuclei <<
346 if (0 < pdg) { <<
347 const G4ParticleDefinition* part = thePartTa <<
348 if(nullptr != part) { <<
349 G4ReactionProduct* theNew = new G4Reaction <<
350 G4ThreeVector dir = G4ThreeVector( 0.0, 0. <<
351 if ( lambdaLV.vect().mag() > CLHEP:: <<
352 dir = lambdaLV.vect().unit(); <<
353 } <<
354 G4double mass = part->GetPDGMass(); <<
355 G4double etot = std::max(lambdaLV.e(), mas <<
356 dir *= std::sqrt((etot - mass)*(etot <<
357 theNew->SetMomentum(dir); <<
358 theNew->SetTotalEnergy(etot); <<
359 theNew->SetFormationTime(theInitialState.G <<
360 theNew->SetCreatorModelID(theInitialState. <<
361 G4ReactionProductVector* v = new G4Reactio <<
362 v->push_back(theNew); <<
363 return v; <<
364 } <<
365 } <<
366 } <<
367 G4double mass = theInitialStatePtr->GetGro <<
368 lambdaF = nL*(fLambdaMass - CLHEP::neutron <<
369 <<
370 // de-excitation with neutrons instead of <<
371 theInitialStatePtr->SetZAandMomentum(lambd <<
372 <<
373 // 4-momentum not used in de-excitation <<
374 lambdaLV *= lambdaF; <<
375 } else if (0 > nL) { <<
376 ++fWarnings; <<
377 if(fWarnings < 0) { <<
378 G4ExceptionDescription ed; <<
379 ed << "Fragment with negative L: Z=" << <<
380 << " Eex/A(MeV)= " << exEnergy/A; <<
381 G4Exception("G4ExcitationHandler::BreakI <<
382 } <<
383 } <<
384 <<
385 // In case A <= 1 the fragment will not perf 238 // In case A <= 1 the fragment will not perform any nucleon emission
386 if (A <= 1 || !isActive) { 239 if (A <= 1 || !isActive) {
387 theResults.push_back( theInitialStatePtr ) 240 theResults.push_back( theInitialStatePtr );
388 241
389 // check if a fragment is stable 242 // check if a fragment is stable
390 } else if (exEnergy < minExcitation && nist- << 243 } else if(exEnergy < minExcitation && nist->GetIsotopeAbundance(Z, A) > 0.0) {
391 theResults.push_back( theInitialStatePtr ) 244 theResults.push_back( theInitialStatePtr );
392 245
393 // JMQ 150909: first step in de-excitation 246 // JMQ 150909: first step in de-excitation is treated separately
394 // Fragments after the first step are stor 247 // Fragments after the first step are stored in theEvapList
395 } else { << 248 } else {
396 if ((A<maxAForFermiBreakUp && Z<maxZForFer << 249 if((A<maxAForFermiBreakUp && Z<maxZForFermiBreakUp)
397 || exEnergy <= minEForMultiFrag*A) { << 250 || exEnergy <= minEForMultiFrag*A) {
398 theEvapList.push_back(theInitialStatePtr 251 theEvapList.push_back(theInitialStatePtr);
399 252
400 // Statistical Multifragmentation will tak 253 // Statistical Multifragmentation will take place only once
401 } else { 254 } else {
402 theTempResult = theMultiFragmentation->B 255 theTempResult = theMultiFragmentation->BreakItUp(theInitialState);
403 if (nullptr == theTempResult) { << 256 if(!theTempResult) {
404 theEvapList.push_back(theInitialStatePtr); 257 theEvapList.push_back(theInitialStatePtr);
405 } else { 258 } else {
406 std::size_t nsec = theTempResult->size(); << 259 size_t nsec = theTempResult->size();
407 260
408 // no fragmentation 261 // no fragmentation
409 if (0 == nsec) { << 262 if(0 == nsec) {
410 theEvapList.push_back(theInitialStatePtr); 263 theEvapList.push_back(theInitialStatePtr);
411 264
412 // secondary are produced - sort out secon 265 // secondary are produced - sort out secondary fragments
413 } else { 266 } else {
414 G4bool deletePrimary = true; 267 G4bool deletePrimary = true;
415 for (auto const & ptr : *theTempResult) { << 268 G4FragmentVector::iterator j;
416 if (ptr == theInitialStatePtr) { deleteP << 269 for (j = theTempResult->begin(); j != theTempResult->end(); ++j) {
417 SortSecondaryFragment(ptr); << 270
>> 271 if((*j) == theInitialStatePtr) { deletePrimary = false; }
>> 272 A = (*j)->GetA_asInt();
>> 273
>> 274 // gamma, p, n
>> 275 if(A <= 1) {
>> 276 theResults.push_back(*j);
>> 277
>> 278 // Analyse fragment A > 1
>> 279 } else {
>> 280 G4double exEnergy1 = (*j)->GetExcitationEnergy();
>> 281
>> 282 // cold fragments
>> 283 if(exEnergy1 < minExcitation) {
>> 284 Z = (*j)->GetZ_asInt();
>> 285 if(nist->GetIsotopeAbundance(Z, A) > 0.0) {
>> 286 theResults.push_back(*j); // stable fragment
>> 287 } else {
>> 288 theEvapList.push_back(*j);
>> 289 }
>> 290 // hot fragments are unstable
>> 291 } else {
>> 292 theEvapList.push_back(*j);
>> 293 }
>> 294 }
418 } 295 }
419 if (deletePrimary) { delete theInitialStat << 296 if( deletePrimary ) { delete theInitialStatePtr; }
420 } 297 }
421 delete theTempResult; // end multifragmentat 298 delete theTempResult; // end multifragmentation
422 } 299 }
423 } 300 }
424 } 301 }
425 if (fVerbose > 2) { << 302 if(fVerbose > 2) {
426 G4cout << "## After first step of handler << 303 G4cout << "## After first step " << theEvapList.size() << " for evap; "
427 << " for evap; " << 304 << thePhotoEvapList.size() << " for photo-evap; "
428 << theResults.size() << " results. " << G 305 << theResults.size() << " results. " << G4endl;
429 } 306 }
430 // ----------------------------------- 307 // -----------------------------------
431 // FermiBreakUp and De-excitation loop 308 // FermiBreakUp and De-excitation loop
432 // ----------------------------------- 309 // -----------------------------------
433 310
434 static const G4int countmax = 1000; 311 static const G4int countmax = 1000;
435 std::size_t kk; << 312 G4Fragment* frag;
>> 313 size_t kk;
436 for (kk=0; kk<theEvapList.size(); ++kk) { 314 for (kk=0; kk<theEvapList.size(); ++kk) {
437 G4Fragment* frag = theEvapList[kk]; << 315 frag = theEvapList[kk];
438 if (fVerbose > 3) { << 316 if(fVerbose > 2) {
439 G4cout << "Next evaporate: " << G4endl; 317 G4cout << "Next evaporate: " << G4endl;
440 G4cout << *frag << G4endl; 318 G4cout << *frag << G4endl;
441 } 319 }
442 if (kk >= countmax) { << 320 if(kk >= countmax) {
443 G4ExceptionDescription ed; 321 G4ExceptionDescription ed;
444 ed << "Infinite loop in the de-excitatio 322 ed << "Infinite loop in the de-excitation module: " << kk
445 << " iterations \n" 323 << " iterations \n"
446 << " Initial fragment: \n" << theIniti 324 << " Initial fragment: \n" << theInitialState
447 << "\n Current fragment: \n" << *frag; 325 << "\n Current fragment: \n" << *frag;
448 G4Exception("G4ExcitationHandler::BreakI 326 G4Exception("G4ExcitationHandler::BreakItUp","had0333",FatalException,
449 ed,"Stop execution"); 327 ed,"Stop execution");
450 328
451 } 329 }
452 A = frag->GetA_asInt(); 330 A = frag->GetA_asInt();
453 Z = frag->GetZ_asInt(); 331 Z = frag->GetZ_asInt();
454 results.clear(); 332 results.clear();
455 if (fVerbose > 2) { << 333
456 G4cout << "G4ExcitationHandler# " << kk <<
457 << " Eex(MeV)= " << frag->GetExci <<
458 } <<
459 // Fermi Break-Up 334 // Fermi Break-Up
460 if (theFermiModel->IsApplicable(Z, A, frag << 335 if(theFermiModel->IsApplicable(Z, A, frag->GetExcitationEnergy())) {
461 theFermiModel->BreakFragment(&results, f 336 theFermiModel->BreakFragment(&results, frag);
462 std::size_t nsec = results.size(); << 337 size_t nsec = results.size();
463 if (fVerbose > 2) { G4cout << "FermiBrea << 338 if(fVerbose > 2) { G4cout << "FermiBreakUp Nsec= " << nsec << G4endl; }
464 339
465 // FBU takes care to delete input fragme 340 // FBU takes care to delete input fragment or add it to the results
466 // The secondary may be excited - photo- 341 // The secondary may be excited - photo-evaporation should be applied
467 if (1 < nsec) { << 342 for(size_t j=0; j<nsec; ++j) {
468 for (auto const & res : results) { << 343 exEnergy = results[j]->GetExcitationEnergy();
469 SortSecondaryFragment(res); << 344 if(exEnergy < minExcitation) { theResults.push_back(results[j]); }
470 } << 345 else { thePhotoEvapList.push_back(results[j]); }
471 continue; <<
472 } 346 }
473 // evaporation will be applied << 347 continue;
474 } 348 }
475 // apply Evaporation, residual nucleus is 349 // apply Evaporation, residual nucleus is always added to the results
476 // photon evaporation is possible <<
477 theEvaporation->BreakFragment(&results, fr 350 theEvaporation->BreakFragment(&results, frag);
478 if (fVerbose > 3) { << 351 size_t nsec = results.size();
479 G4cout << "Evaporation Nsec= " << result << 352 if(fVerbose > 2) { G4cout << "Evaporation Nsec= " << nsec << G4endl; }
480 } << 353
481 if (0 == results.size()) { << 354 // no evaporation
>> 355 if(1 >= nsec) {
482 theResults.push_back(frag); 356 theResults.push_back(frag);
483 } else { << 357 continue;
484 SortSecondaryFragment(frag); <<
485 } 358 }
486 359
487 // Sort out secondary fragments 360 // Sort out secondary fragments
488 for (auto const & res : results) { << 361 for (size_t j = 0; j<nsec; ++j) {
489 if(fVerbose > 4) { << 362 if(fVerbose > 3) {
490 G4cout << "Evaporated product #" << *res << << 363 G4cout << "Evaporated product #" << j << G4endl;
>> 364 G4cout << results[j] << G4endl;
>> 365 }
>> 366 A = results[j]->GetA_asInt();
>> 367 if(A <= 1) {
>> 368 theResults.push_back(results[j]); // gamma, p, n
>> 369 continue;
>> 370 }
>> 371 exEnergy = results[j]->GetExcitationEnergy();
>> 372
>> 373 // hot fragment
>> 374 if(exEnergy >= minExcitation) {
>> 375 theEvapList.push_back(results[j]);
>> 376
>> 377 // cold fragment
>> 378 } else {
>> 379 Z = results[j]->GetZ_asInt();
>> 380
>> 381 // natural isotope
>> 382 if(nist->GetIsotopeAbundance(Z, A) > 0.0) {
>> 383 theResults.push_back(results[j]); // stable fragment
>> 384
>> 385 } else {
>> 386 theEvapList.push_back(results[j]);
>> 387 }
491 } 388 }
492 SortSecondaryFragment(res); <<
493 } // end of loop on secondary 389 } // end of loop on secondary
494 } // end of the loop over theEvapList 390 } // end of the loop over theEvapList
495 if (fVerbose > 2) { << 391 if(fVerbose > 2) {
496 G4cout << "## After 2nd step of handler " << 392 G4cout << "## After 2nd step " << theEvapList.size() << " was evap; "
497 << " was evap; " << 393 << thePhotoEvapList.size() << " for photo-evap; "
>> 394 << theResults.size() << " results. " << G4endl;
>> 395 }
>> 396 // -----------------------
>> 397 // Photon-Evaporation loop
>> 398 // -----------------------
>> 399
>> 400 // at this point only photon evaporation is possible
>> 401 size_t kkmax = thePhotoEvapList.size();
>> 402 for (kk=0; kk<kkmax; ++kk) {
>> 403 frag = thePhotoEvapList[kk];
>> 404 if(fVerbose > 2) {
>> 405 G4cout << "Next photon evaporate: " << thePhotonEvaporation << G4endl;
>> 406 G4cout << *frag << G4endl;
>> 407 }
>> 408 exEnergy = frag->GetExcitationEnergy();
>> 409
>> 410 // photon de-excitation only for hot fragments
>> 411 if(exEnergy > minExcitation) {
>> 412 thePhotonEvaporation->BreakUpChain(&theResults, frag);
>> 413 }
>> 414
>> 415 // primary fragment is kept
>> 416 theResults.push_back(frag);
>> 417
>> 418 } // end of photon-evaporation loop
>> 419 if(fVerbose > 2) {
>> 420 G4cout << "## After 3d step " << theEvapList.size() << " was evap; "
>> 421 << thePhotoEvapList.size() << " was photo-evap; "
498 << theResults.size() << " results. " << G 422 << theResults.size() << " results. " << G4endl;
499 } 423 }
500 G4ReactionProductVector * theReactionProduct 424 G4ReactionProductVector * theReactionProductVector =
501 new G4ReactionProductVector(); 425 new G4ReactionProductVector();
502 426
503 // MAC (24/07/08) 427 // MAC (24/07/08)
504 // To optimise the storing speed, we reserve << 428 // To optimise the storing speed, we reserve space in memory for the vector
505 // in memory for the vector <<
506 theReactionProductVector->reserve( theResult 429 theReactionProductVector->reserve( theResults.size() );
507 430
508 if (fVerbose > 1) { << 431 G4int theFragmentA, theFragmentZ;
>> 432
>> 433 if(fVerbose > 1) {
509 G4cout << "### ExcitationHandler provides 434 G4cout << "### ExcitationHandler provides " << theResults.size()
510 << " evaporated products:" << G4endl; 435 << " evaporated products:" << G4endl;
511 } 436 }
512 G4LorentzVector partOfLambdaLV; << 437 kkmax = theResults.size();
513 if ( nL > 0 ) partOfLambdaLV = lambdaLV/(G4d << 438 for (kk=0; kk<kkmax; ++kk) {
514 for (auto const & frag : theResults) { << 439 frag = theResults[kk];
515 G4LorentzVector lv0 = frag->GetMomentum(); <<
516 G4double etot = lv0.e(); <<
517 440
518 // in the case of dummy de-excitation, exc 441 // in the case of dummy de-excitation, excitation energy is transfered
519 // into kinetic energy of output ion << 442 // into kinetic energy
520 if (!isActive) { << 443 if(!isActive && 0 == kk) {
521 G4double mass = frag->GetGroundStateMass 444 G4double mass = frag->GetGroundStateMass();
522 G4double ptot = lv0.vect().mag(); << 445 G4double ptot = (frag->GetMomentum()).vect().mag();
>> 446 G4double etot = (frag->GetMomentum()).e();
523 G4double fac = (etot <= mass || 0.0 == 447 G4double fac = (etot <= mass || 0.0 == ptot) ? 0.0
524 : std::sqrt((etot - mass)*(etot + mass))/pto 448 : std::sqrt((etot - mass)*(etot + mass))/ptot;
525 G4LorentzVector lv((frag->GetMomentum()) 449 G4LorentzVector lv((frag->GetMomentum()).px()*fac,
526 (frag->GetMomentum()).py()*fac, 450 (frag->GetMomentum()).py()*fac,
527 (frag->GetMomentum()).pz()*fac, etot); 451 (frag->GetMomentum()).pz()*fac, etot);
528 frag->SetMomentum(lv); 452 frag->SetMomentum(lv);
529 } 453 }
530 if (fVerbose > 3) { << 454 if(fVerbose > 1) {
531 G4cout << *frag; << 455 G4cout << kk << "-th fragment " << frag;
532 if (frag->NuclearPolarization()) { << 456 if(frag->NuclearPolarization()) {
533 G4cout << " " << frag->NuclearPolarization( 457 G4cout << " " << frag->NuclearPolarization();
534 } 458 }
535 G4cout << G4endl; 459 G4cout << G4endl;
>> 460 G4cout << *frag << G4endl;
536 } 461 }
537 462
538 G4int fragmentA = frag->GetA_asInt(); << 463 theFragmentA = frag->GetA_asInt();
539 G4int fragmentZ = frag->GetZ_asInt(); << 464 theFragmentZ = frag->GetZ_asInt();
>> 465 G4double etot= frag->GetMomentum().e();
540 G4double eexc = 0.0; 466 G4double eexc = 0.0;
541 const G4ParticleDefinition* theKindOfFragm 467 const G4ParticleDefinition* theKindOfFragment = nullptr;
542 G4bool isHyperN = false; << 468 if (theFragmentA == 0) { // photon or e-
543 if (fragmentA == 0) { // photon or e <<
544 theKindOfFragment = frag->GetParticleDef 469 theKindOfFragment = frag->GetParticleDefinition();
545 } else if (fragmentA == 1 && fragmentZ == << 470 } else if (theFragmentA == 1 && theFragmentZ == 0) { // neutron
546 theKindOfFragment = theNeutron; << 471 theKindOfFragment = G4Neutron::NeutronDefinition();
547 } else if (fragmentA == 1 && fragmentZ == << 472 } else if (theFragmentA == 1 && theFragmentZ == 1) { // proton
548 theKindOfFragment = theProton; << 473 theKindOfFragment = G4Proton::ProtonDefinition();
549 } else if (fragmentA == 2 && fragmentZ == << 474 } else if (theFragmentA == 2 && theFragmentZ == 1) { // deuteron
550 theKindOfFragment = theDeuteron; << 475 theKindOfFragment = G4Deuteron::DeuteronDefinition();
551 } else if (fragmentA == 3 && fragmentZ == << 476 } else if (theFragmentA == 3 && theFragmentZ == 1) { // triton
552 theKindOfFragment = theTriton; << 477 theKindOfFragment = G4Triton::TritonDefinition();
553 if(0 < nL) { << 478 } else if (theFragmentA == 3 && theFragmentZ == 2) { // helium3
554 const G4ParticleDefinition* p = thePar << 479 theKindOfFragment = G4He3::He3Definition();
555 if(nullptr != p) { << 480 } else if (theFragmentA == 4 && theFragmentZ == 2) { // alpha
556 theKindOfFragment = p; << 481 theKindOfFragment = G4Alpha::AlphaDefinition();;
557 isHyperN = true; <<
558 --nL; <<
559 } <<
560 } <<
561 } else if (fragmentA == 3 && fragmentZ == <<
562 theKindOfFragment = theHe3; <<
563 } else if (fragmentA == 4 && fragmentZ == <<
564 theKindOfFragment = theAlpha; <<
565 if (0 < nL) { <<
566 const G4ParticleDefinition* p = thePar <<
567 if(nullptr != p) { <<
568 theKindOfFragment = p; <<
569 isHyperN = true; <<
570 --nL; <<
571 } <<
572 } <<
573 } else { 482 } else {
574 483
575 // fragment 484 // fragment
576 eexc = frag->GetExcitationEnergy(); 485 eexc = frag->GetExcitationEnergy();
577 G4int idxf = frag->GetFloatingLevelNumbe 486 G4int idxf = frag->GetFloatingLevelNumber();
578 if (eexc < minExcitation) { << 487 if(eexc < minExcitation) {
579 eexc = 0.0; 488 eexc = 0.0;
580 idxf = 0; 489 idxf = 0;
581 } 490 }
582 491
583 theKindOfFragment = theTableOfIons->GetI << 492 theKindOfFragment = theTableOfIons->GetIon(theFragmentZ,theFragmentA,eexc,
584 493 G4Ions::FloatLevelBase(idxf));
585 if (fVerbose > 3) { << 494 if(fVerbose > 1) {
586 G4cout << "### EXCH: Find ion Z= " << fragme << 495 G4cout << "### EXCH: Find ion Z= " << theFragmentZ << " A= " << theFragmentA
587 << " A= " << fragmentA <<
588 << " Eexc(MeV)= " << eexc/MeV << " id 496 << " Eexc(MeV)= " << eexc/MeV << " idx= " << idxf
589 << " " << theKindOfFragment->GetParti << 497 << " " << theKindOfFragment << G4endl;
590 << G4endl; <<
591 } 498 }
592 } 499 }
593 // fragment identified 500 // fragment identified
594 if (nullptr != theKindOfFragment) { << 501 if(theKindOfFragment) {
595 G4ReactionProduct * theNew = new G4React 502 G4ReactionProduct * theNew = new G4ReactionProduct(theKindOfFragment);
596 if (isHyperN) { << 503 theNew->SetMomentum(frag->GetMomentum().vect());
597 G4LorentzVector lv = lv0 + partOfLambd <<
598 G4ThreeVector dir = lv.vect().unit(); <<
599 G4double mass = theKindOfFragment->Get <<
600 etot = std::max(lv.e(), mass); <<
601 G4double ptot = std::sqrt((etot - mass <<
602 dir *= ptot; <<
603 theNew->SetMomentum(dir); <<
604 // remaining not compensated 4-momentum <<
605 lambdaLV += (lv0 - G4LorentzVector(dir <<
606 } else { <<
607 theNew->SetMomentum(lv0.vect()); <<
608 } <<
609 theNew->SetTotalEnergy(etot); 504 theNew->SetTotalEnergy(etot);
610 theNew->SetFormationTime(frag->GetCreati 505 theNew->SetFormationTime(frag->GetCreationTime());
611 theNew->SetCreatorModelID(frag->GetCreat << 506 if(theKindOfFragment == electron) { theNew->SetCreatorModel(icID); }
612 theReactionProductVector->push_back(theN 507 theReactionProductVector->push_back(theNew);
613 508
614 // fragment not found out ground state i 509 // fragment not found out ground state is created
615 } else { 510 } else {
616 theKindOfFragment = 511 theKindOfFragment =
617 theTableOfIons->GetIon(fragmentZ,fragmentA,0 << 512 theTableOfIons->GetIon(theFragmentZ,theFragmentA,0.0,noFloat,0);
618 if (theKindOfFragment) { << 513 if(theKindOfFragment) {
619 G4ThreeVector mom(0.0,0.0,0.0); 514 G4ThreeVector mom(0.0,0.0,0.0);
620 G4double ionmass = theKindOfFragment->GetPDG 515 G4double ionmass = theKindOfFragment->GetPDGMass();
621 if (etot <= ionmass) { << 516 if(etot <= ionmass) {
622 etot = ionmass; 517 etot = ionmass;
623 } else { 518 } else {
624 G4double ptot = std::sqrt((etot - ionmass) 519 G4double ptot = std::sqrt((etot - ionmass)*(etot + ionmass));
625 mom = (frag->GetMomentum().vect().unit())* 520 mom = (frag->GetMomentum().vect().unit())*ptot;
626 } 521 }
627 G4ReactionProduct * theNew = new G4ReactionP 522 G4ReactionProduct * theNew = new G4ReactionProduct(theKindOfFragment);
628 theNew->SetMomentum(mom); 523 theNew->SetMomentum(mom);
629 theNew->SetTotalEnergy(etot); 524 theNew->SetTotalEnergy(etot);
630 theNew->SetFormationTime(frag->GetCreationTi 525 theNew->SetFormationTime(frag->GetCreationTime());
631 theNew->SetCreatorModelID(frag->GetCreatorMo <<
632 theReactionProductVector->push_back(theNew); 526 theReactionProductVector->push_back(theNew);
633 if (fVerbose > 3) { << 527 if(fVerbose > 2) {
634 G4cout << " ground state, energy << 528 G4cout << "### Find ion Z= " << theFragmentZ << " A= " << theFragmentA
635 << etot << G4endl; << 529 << " ground state, energy corrected E(MeV)= " << etot << G4endl;
636 } 530 }
637 } 531 }
638 } 532 }
639 delete frag; 533 delete frag;
>> 534 if(fVerbose > 1) { G4cout << "G4Fragment #" << kk << " is deleted" << G4endl; }
640 } 535 }
641 // remaining lambdas are free; conserve quan << 536 if(fVerbose > 2) {
642 // not 4-momentum <<
643 if (0 < nL) { <<
644 G4ThreeVector dir = G4ThreeVector(0.0, 0.0 <<
645 if (lambdaLV.vect().mag() > CLHEP::eV) { <<
646 dir = lambdaLV.vect().unit(); <<
647 } <<
648 G4double etot = std::max(lambdaLV.e()/(G4d <<
649 dir *= std::sqrt((etot - fLambdaMass)*(eto <<
650 for (G4int i=0; i<nL; ++i) { <<
651 G4ReactionProduct* theNew = new G4Reacti <<
652 theNew->SetMomentum(dir); <<
653 theNew->SetTotalEnergy(etot); <<
654 theNew->SetFormationTime(theInitialState <<
655 theNew->SetCreatorModelID(theInitialStat <<
656 theReactionProductVector->push_back(theN <<
657 } <<
658 } <<
659 if (fVerbose > 3) { <<
660 G4cout << "@@@@@@@@@@ End G4Excitation Han 537 G4cout << "@@@@@@@@@@ End G4Excitation Handler "<< G4endl;
661 } 538 }
662 return theReactionProductVector; 539 return theReactionProductVector;
663 } 540 }
664 541
665 void G4ExcitationHandler::ModelDescription(std 542 void G4ExcitationHandler::ModelDescription(std::ostream& outFile) const
666 { 543 {
667 outFile << "G4ExcitationHandler description\ 544 outFile << "G4ExcitationHandler description\n"
668 << "This class samples de-excitation of ex 545 << "This class samples de-excitation of excited nucleus using\n"
669 << "Fermi Break-up model for light fragmen 546 << "Fermi Break-up model for light fragments (Z < 9, A < 17), "
670 << "evaporation, fission, and photo-evapor 547 << "evaporation, fission, and photo-evaporation models. Evaporated\n"
671 << "particle may be proton, neutron, and o 548 << "particle may be proton, neutron, and other light fragment \n"
672 << "(Z < 13, A < 29). During photon evapor 549 << "(Z < 13, A < 29). During photon evaporation produced gamma \n"
673 << "or electrons due to internal conversio 550 << "or electrons due to internal conversion \n";
674 } 551 }
>> 552
>> 553
>> 554
675 555
676 556
677 557
678 558