| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // INCL++ intra-nuclear cascade model
27 // Alain Boudard, CEA-Saclay, France
28 // Joseph Cugnon, University of Liege, Belgium
29 // Jean-Christophe David, CEA-Saclay, France
30 // Pekka Kaitaniemi, CEA-Saclay, France, and Helsinki Institute of Physics, Finland
31 // Sylvie Leray, CEA-Saclay, France
32 // Davide Mancusi, CEA-Saclay, France
33 //
34 #define INCLXX_IN_GEANT4_MODE 1
35
36 #include "globals.hh"
37
38 #include "G4INCLStore.hh"
39 #include <fstream>
40 #include "G4INCLLogger.hh"
41 #include "G4INCLCluster.hh"
42
43 namespace G4INCL {
44
45 Store::Store(Config const * const config) :
46 loadedA(0),
47 loadedZ(0),
48 loadedStoppingTime(0.),
49 theConfig(config)
50 {
51 }
52
53 Store::~Store() {
54 theBook.reset();
55 clear();
56 }
57
58 void Store::add(Particle *p) {
59 inside.push_back(p);
60 }
61
62 void Store::add(ParticleList const &pL) {
63 inside.insert(inside.end(), pL.begin(), pL.end());
64 }
65
66 void Store::addParticleEntryAvatar(IAvatar *a) {
67 // Add the avatar to the avatar map
68 avatarList.push_back(a);
69
70 ParticleList pList = a->getParticles();
71 for(ParticleIter i=pList.begin(), e=pList.end(); i!=e; ++i) {
72 addIncomingParticle((*i));
73 // Connect each particle to the avatar
74 connectAvatarToParticle(a, *i);
75 }
76 }
77
78 void Store::addParticleEntryAvatars(IAvatarList const &al) {
79 for(IAvatarIter a=al.begin(), e=al.end(); a!=e; ++a)
80 addParticleEntryAvatar(*a);
81 }
82
83 void Store::add(IAvatar *a) {
84 // Add the avatar to the avatar map
85 avatarList.push_back(a);
86
87 ParticleList pList = a->getParticles();
88 for(ParticleIter i=pList.begin(), e=pList.end(); i!=e; ++i) {
89 // Connect each particle to the avatar
90 connectAvatarToParticle(a, *i);
91 }
92
93 }
94
95 void Store::addIncomingParticle(Particle * const p) {
96 incoming.push_back(p);
97 }
98
99 void Store::connectAvatarToParticle(IAvatar * const a, Particle * const p) {
100 particleAvatarConnections.insert(PAPair(p,a));
101 }
102
103 void Store::disconnectAvatarFromParticle(IAvatar * const a, Particle * const p) {
104 PAIterPair iterPair = particleAvatarConnections.equal_range(p);
105 for(PAIter i=iterPair.first, last=iterPair.second; i!=last; ++i) {
106 if(i->second==a) {
107 particleAvatarConnections.erase(i);
108 return;
109 }
110 }
111 INCL_WARN("Loop in Store::disconnectAvatarFromParticle fell through." << std::endl
112 << "This indicates an inconsistent state of the particleAvatarConnections map." << std::endl);
113 }
114
115 void Store::removeAvatar(IAvatar * const avatar) {
116 // Disconnect the avatar from particles
117 ParticleList particlesRelatedToAvatar = avatar->getParticles();
118 for(ParticleIter particleIter = particlesRelatedToAvatar.begin(), e = particlesRelatedToAvatar.end();
119 particleIter != e; ++particleIter) {
120 disconnectAvatarFromParticle(avatar, *particleIter);
121 }
122
123 // Remove the avatar itself
124 avatarList.remove(avatar);
125 }
126
127 void Store::particleHasBeenUpdated(Particle * const particle) {
128 PAIterPair iterPair = particleAvatarConnections.equal_range(particle);
129 for(PAIter i=iterPair.first, last=iterPair.second; i!=last; ++i) {
130 avatarsToBeRemoved.insert(i->second);
131 }
132 }
133
134 void Store::removeScheduledAvatars() {
135 for(ASIter a=avatarsToBeRemoved.begin(), e=avatarsToBeRemoved.end(); a!=e; ++a) {
136 removeAvatar(*a);
137 delete *a;
138 }
139 avatarsToBeRemoved.clear();
140 }
141
142 IAvatar* Store::findSmallestTime() {
143 if(avatarList.empty()) return NULL;
144
145 #ifdef INCL_AVATAR_SEARCH_FullSort
146
147 /* Full sort algorithm.
148 *
149 * Simple, but guaranteed to work.
150 */
151 avatarList.sort(Store::avatarComparisonPredicate);
152 IAvatar *avatar = avatarList.front();
153
154 #elif defined(INCL_AVATAR_SEARCH_MinElement)
155
156 /* Algorithm provided by the C++ stdlib. */
157 IAvatar *avatar = *(std::min_element(avatarList.begin(), avatarList.end(),
158 Store::avatarComparisonPredicate));
159
160 #else
161 #error Unrecognized INCL_AVATAR_SEARCH. Allowed values are: FullSort, MinElement.
162 #endif
163
164 removeAvatar(avatar);
165 return avatar;
166 }
167
168 void Store::timeStep(G4double step) {
169 for(ParticleIter particleIter = inside.begin(), particleEnd=inside.end();
170 particleIter != particleEnd; ++particleIter) {
171 (*particleIter)->propagate(step);
172 }
173 }
174
175 void Store::particleHasBeenEjected(Particle * const p) {
176 particleHasBeenUpdated(p);
177 // The particle will be destroyed when destroying the Store
178 inside.remove(p);
179 }
180
181 void Store::particleHasBeenDestroyed(Particle * const p) {
182 particleHasBeenUpdated(p);
183 // Have to destroy the particle here, the Store will forget about it
184 inside.remove(p);
185 delete p;
186 }
187
188 void Store::particleHasEntered(Particle * const particle) {
189 removeFromIncoming(particle);
190 add(particle);
191 }
192
193 void Store::clearAvatars() {
194 for(IAvatarIter iter = avatarList.begin(), e = avatarList.end(); iter != e; ++iter) {
195 delete *iter;
196 }
197
198 particleAvatarConnections.clear();
199 avatarList.clear();
200 avatarsToBeRemoved.clear();
201 }
202
203 void Store::clear() {
204 clearAvatars();
205
206 clearInside();
207 clearOutgoing();
208
209 if( incoming.size() != 0 ) {
210 INCL_WARN("Incoming list is not empty when Store::clear() is called" << '\n');
211 }
212 incoming.clear();
213
214 }
215
216 void Store::clearInside() {
217 for(ParticleIter iter=inside.begin(), e=inside.end(); iter!=e; ++iter) {
218 delete *iter;
219 }
220 inside.clear();
221 }
222
223 void Store::clearOutgoing() {
224 for(ParticleIter iter=outgoing.begin(), e=outgoing.end(); iter!=e; ++iter) {
225 if((*iter)->isCluster()) {
226 Cluster *c = dynamic_cast<Cluster *>(*iter);
227 // assert(c);
228 #ifdef INCLXX_IN_GEANT4_MODE
229 if(!c)
230 continue;
231 #endif
232 c->deleteParticles();
233 }
234 delete (*iter);
235 }
236 outgoing.clear();
237 }
238
239 void Store::loadParticles(std::string const &filename) {
240 clear();
241 G4int projectileA, projectileZ, A, Z;
242 G4double stoppingTime, cutNN;
243 G4int ID, type, isParticipant;
244 G4double x, y, z;
245 G4double px, py, pz, E, v;
246
247 std::ifstream in(filename.c_str());
248 in >> projectileA >> projectileZ >> A >> Z >> stoppingTime >> cutNN;
249 loadedA = A;
250 loadedZ = Z;
251 loadedStoppingTime = stoppingTime;
252
253 while(1) { /* Loop checking, 10.07.2015, D.Mancusi */
254 in >> ID >> type >> isParticipant >> x >> y >> z >> px >> py >> pz >> E >> v;
255 if(!in.good()) break;
256 ParticleType t;
257 if(type == 1) {
258 t = Proton;
259 }
260 else if(type == -1) {
261 t = Neutron;
262 }
263 else {
264 INCL_FATAL("Unrecognized particle type while loading particles; type=" << type << '\n');
265 t = UnknownParticle;
266 }
267
268 Particle *p = new Particle(t, E, ThreeVector(px, py, pz),
269 ThreeVector(x, y, z));
270 p->setPotentialEnergy(v);
271 if(isParticipant == 1) {
272 p->makeParticipant();
273 theBook.incrementCascading();
274 }
275 add(p);
276 }
277
278 in.close();
279 }
280
281 std::string Store::printParticleConfiguration() {
282 std::stringstream ss;
283 G4int A = 0, Z = 0;
284 for(ParticleIter i=inside.begin(), e=inside.end(); i!=e; ++i) {
285 if((*i)->getType() == Proton) {
286 A++;
287 Z++;
288 }
289 if((*i)->getType() == Neutron) {
290 A++;
291 }
292 }
293 // Note: Projectile A and Z are set to 0 (we don't really know
294 // anything about them at this point).
295 ss << "0 0 " << A << " " << Z << " "
296 << "100.0" << " "
297 << "0.0" << '\n';
298
299 for(ParticleIter i=inside.begin(), e=inside.end(); i!=e; ++i) {
300 G4int ID = (G4int)(*i)->getID();
301 G4int type = 0;
302 if((*i)->getType() == Proton) {
303 type = 1;
304 }
305 if((*i)->getType() == Neutron) {
306 type = -1;
307 }
308
309 G4int isParticipant = 0;
310 if((*i)->isParticipant()) {
311 isParticipant = 1;
312 }
313
314 G4double x = (*i)->getPosition().getX();
315 G4double y = (*i)->getPosition().getY();
316 G4double z = (*i)->getPosition().getZ();
317 G4double E = (*i)->getEnergy();
318 G4double px = (*i)->getMomentum().getX();
319 G4double py = (*i)->getMomentum().getY();
320 G4double pz = (*i)->getMomentum().getZ();
321 G4double V = (*i)->getPotentialEnergy();
322
323 ss << ID << " " << type << " " << isParticipant << " "
324 << x << " " << y << " " << z << " "
325 << px << " " << py << " " << pz << " "
326 << E << " " << V << '\n';
327 }
328
329 return ss.str();
330 }
331
332 void Store::writeParticles(std::string const &filename) {
333 std::ofstream out(filename.c_str());
334 out << printParticleConfiguration();
335 out.close();
336 }
337
338 std::string Store::printAvatars() {
339 std::stringstream ss;
340 for(IAvatarIter i = avatarList.begin(), e = avatarList.end(); i != e; ++i) {
341 ss << (*i)->toString() << '\n';
342 }
343 return ss.str();
344 }
345
346 G4bool Store::containsCollisions() const {
347 for(IAvatarIter i = avatarList.begin(), e = avatarList.end(); i != e; ++i)
348 if((*i)->getType()==CollisionAvatarType) return true;
349 return false;
350 }
351
352 }
353