00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033 #define INCLXX_IN_GEANT4_MODE 1
00034
00035 #include "globals.hh"
00036
00037 #include "G4INCLStore.hh"
00038 #include <fstream>
00039 #include "G4INCLLogger.hh"
00040 #include "G4INCLCluster.hh"
00041
00042 namespace G4INCL {
00043
00044 Store::Store(Config const * const config) :
00045 theBook(new Book),
00046 loadedA(0),
00047 loadedZ(0),
00048 loadedStoppingTime(0.),
00049 theConfig(config)
00050 {
00051 }
00052
00053 Store::~Store() {
00054 theBook->reset();
00055 delete theBook;
00056 theBook = 0;
00057 clear();
00058 }
00059
00060 void Store::add(Particle *p) {
00061 const long ID = p->getID();
00062 particles[ID]=p;
00063 inside.push_back(p);
00064
00065 if(particleAvatarConnections.find(ID)==particleAvatarConnections.end()) {
00066 std::vector<long> *aIDs = new std::vector<long>;
00067 particleAvatarConnections[ID] = aIDs;
00068 }
00069 }
00070
00071 void Store::addParticleEntryAvatar(IAvatar *a) {
00072
00073 avatars[a->getID()]=a;
00074 avatarList.push_back(a);
00075
00076 ParticleList pList = a->getParticles();
00077 for(ParticleIter i = pList.begin(); i != pList.end(); ++i) {
00078 addIncomingParticle((*i));
00079
00080 connectParticleAndAvatar((*i)->getID(), a->getID());
00081 }
00082 }
00083
00084 void Store::addParticleEntryAvatars(IAvatarList const &al) {
00085 for(IAvatarIter a=al.begin(); a!=al.end(); ++a)
00086 addParticleEntryAvatar(*a);
00087 }
00088
00089 void Store::add(IAvatar *a) {
00090
00091 avatars[a->getID()]=a;
00092 avatarList.push_back(a);
00093
00094 ParticleList pList = a->getParticles();
00095 for(ParticleIter i = pList.begin(); i != pList.end(); ++i) {
00096
00097
00098
00099
00100 if(particles.find((*i)->getID()) == particles.end()) {
00101 ERROR("Avatar was added before related particles. This is probably a bug." << std::endl);
00102 add((*i));
00103 }
00104
00105 connectParticleAndAvatar((*i)->getID(), a->getID());
00106 }
00107
00108 }
00109
00110 void Store::addIncomingParticle(Particle * const p) {
00111 incoming.push_back(p);
00112 }
00113
00114 void Store::connectParticleAndAvatar(long particleID, long avatarID) {
00115 std::map<long, std::vector<long>* >::const_iterator iter = particleAvatarConnections.find(particleID);
00116
00117 if(iter!=particleAvatarConnections.end()) {
00118 std::vector<long> *p = iter->second;
00119 p->push_back(avatarID);
00120 } else {
00121 std::vector<long> *p = new std::vector<long>;
00122 p->push_back(avatarID);
00123 particleAvatarConnections[particleID]=p;
00124 }
00125 }
00126
00127 void Store::removeAvatarFromParticle(long particleID, long avatarID) {
00128 std::vector<long>* theseAvatars = particleAvatarConnections.find(particleID)->second;
00129 std::vector<long>* newAvatars = new std::vector<long>();
00130 for(std::vector<long>::const_iterator iter = theseAvatars->begin();
00131 iter != theseAvatars->end(); ++iter) {
00132 if((*iter) != avatarID) {
00133 newAvatars->push_back((*iter));
00134 }
00135 }
00136 delete theseAvatars;
00137 particleAvatarConnections[particleID] = newAvatars;
00138 }
00139
00140 void Store::removeAvatarByID(long ID) {
00141
00142 IAvatar *avatar = avatars.find(ID)->second;
00143 ParticleList particlesRelatedToAvatar = avatar->getParticles();
00144 for(ParticleIter particleIDiter
00145 = particlesRelatedToAvatar.begin();
00146 particleIDiter != particlesRelatedToAvatar.end(); ++particleIDiter) {
00147 removeAvatarFromParticle((*particleIDiter)->getID(), ID);
00148 }
00149
00150 #ifdef INCL_AVATAR_SEARCH_INCLSort
00151
00152 std::list<IAvatarIter>::iterator it=binaryIterSearch(avatars.find(ID)->second);
00153 if(it != avatarIterList.end())
00154 avatarIterList.erase(it);
00155 #endif
00156
00157
00158 avatarList.remove(avatar);
00159 avatars.erase(ID);
00160 }
00161
00162 void Store::particleHasBeenUpdated(long particleID) {
00163 std::vector<long> temp_aIDs;
00164 std::vector<long> *aIDs = particleAvatarConnections.find(particleID)->second;
00165 for(std::vector<long>::iterator i = aIDs->begin();
00166 i != aIDs->end(); ++i) {
00167 temp_aIDs.push_back((*i));
00168 }
00169
00170 for(std::vector<long>::iterator i = temp_aIDs.begin();
00171 i != temp_aIDs.end(); ++i) {
00172 IAvatar *tmpAvatar = avatars.find(*i)->second;
00173 removeAvatarByID((*i));
00174 delete tmpAvatar;
00175 }
00176 }
00177
00178 #ifdef INCL_AVATAR_SEARCH_INCLSort
00179 std::list<IAvatarIter>::iterator Store::binaryIterSearch(IAvatar const * const avatar) {
00180 std::list<IAvatarIter>::iterator it;
00181 std::iterator_traits<std::list<IAvatarIter>::iterator>::difference_type count, step;
00182 std::list<IAvatarIter>::iterator first = avatarIterList.begin();
00183 std::list<IAvatarIter>::iterator last = avatarIterList.end();
00184 const G4double avatarTime = avatar->getTime();
00185 count = distance(first,last);
00186 while (count>0)
00187 {
00188 it = first; step=count/2; advance(it,step);
00189 if ((**it)->getTime()>avatarTime)
00190 { first=++it; count-=step+1; }
00191 else count=step;
00192 }
00193 if(first!=last && (**first)->getID()==avatar->getID())
00194 return first;
00195 else
00196 return last;
00197 }
00198 #endif
00199
00200 void Store::removeAvatarsFromParticle(long ID) {
00201 std::vector<long> *relatedAvatars = particleAvatarConnections.find(ID)->second;
00202 for(std::vector<long>::const_iterator i = relatedAvatars->begin();
00203 i != relatedAvatars->end(); ++i) {
00204 removeAvatarByID((*i));
00205 }
00206 delete relatedAvatars;
00207 particleAvatarConnections.erase(ID);
00208 }
00209
00210 IAvatar* Store::findSmallestTime() {
00211 if(avatarList.empty()) return NULL;
00212
00213 #ifdef INCL_AVATAR_SEARCH_FullSort
00214
00215
00216
00217
00218
00219 avatarList.sort(Store::avatarComparisonPredicate);
00220 IAvatar *avatar = avatarList.front();
00221
00222 #elif defined(INCL_AVATAR_SEARCH_INCLSort)
00223
00224
00225
00226
00227
00228
00229
00230
00231 IAvatarIter best;
00232 if(avatarIterList.empty())
00233 best = avatarList.begin();
00234 else
00235 best = avatarIterList.back();
00236 G4double bestTime = (*best)->getTime();
00237 IAvatarIter a = best;
00238
00239 for(++a; a!=avatarList.end(); ++a)
00240 if((*a)->getTime() < bestTime) {
00241 best = a;
00242 bestTime = (*best)->getTime();
00243 avatarIterList.push_back(best);
00244 }
00245 IAvatar *avatar = *best;
00246
00247 #elif defined(INCL_AVATAR_SEARCH_MinElement)
00248
00249
00250 IAvatar *avatar = *(std::min_element(avatarList.begin(), avatarList.end(),
00251 Store::avatarComparisonPredicate));
00252
00253 #else
00254 #error Unrecognized INCL_AVATAR_SEARCH. Allowed values are: FullSort, INCLSort, MinElement.
00255 #endif
00256
00257 removeAvatarByID(avatar->getID());
00258 return avatar;
00259 }
00260
00261 void Store::timeStep(G4double step) {
00262 for(std::map<long, Particle*>::iterator particleIter
00263 = particles.begin();
00264 particleIter != particles.end(); ++particleIter) {
00265 (*particleIter).second->propagate(step);
00266 }
00267 }
00268
00269 void Store::particleHasBeenEjected(long ID) {
00270 particleHasBeenUpdated(ID);
00271
00272 inside.remove(particles.find(ID)->second);
00273 particles.erase(ID);
00274 delete particleAvatarConnections.find(ID)->second;
00275 particleAvatarConnections.erase(ID);
00276 }
00277
00278 void Store::particleHasBeenDestroyed(long ID) {
00279 particleHasBeenUpdated(ID);
00280
00281 Particle * const toDelete = particles.find(ID)->second;
00282 inside.remove(toDelete);
00283 delete toDelete;
00284 particles.erase(ID);
00285 }
00286
00287 void Store::particleHasEntered(Particle * const particle) {
00288 removeFromIncoming(particle);
00289 add(particle);
00290 }
00291
00292 ParticleList Store::getParticipants() {
00293 WARN("Store::getParticipants is probably slow..." << std::endl);
00294 ParticleList result;
00295 for(std::map<long, Particle*>::iterator iter = particles.begin();
00296 iter != particles.end(); ++iter) {
00297 if((*iter).second->isParticipant()) {
00298 result.push_back((*iter).second);
00299 }
00300 }
00301 return result;
00302 }
00303
00304 ParticleList Store::getSpectators() {
00305 WARN("Store::getSpectators is probably slow..." << std::endl);
00306 ParticleList result;
00307 for(std::map<long, Particle*>::iterator iter = particles.begin();
00308 iter != particles.end(); ++iter) {
00309 if(!(*iter).second->isParticipant()) {
00310 result.push_back((*iter).second);
00311 }
00312 }
00313 return result;
00314 }
00315
00316 void Store::clearAvatars() {
00317 for(std::map<long, IAvatar*>::iterator iter = avatars.begin();
00318 iter != avatars.end(); ++iter) {
00319 delete (*iter).second;
00320 }
00321
00322 for(std::map<long, std::vector<long>*>::iterator iter = particleAvatarConnections.begin();
00323 iter != particleAvatarConnections.end(); ++iter) {
00324 delete (*iter).second;
00325 }
00326
00327 particleAvatarConnections.clear();
00328 avatars.clear();
00329 avatarList.clear();
00330
00331 }
00332
00333 void Store::initialiseParticleAvatarConnections() {
00334 for(ParticleIter ip=inside.begin(); ip!=inside.end(); ++ip) {
00335 std::vector<long> *aIDs = new std::vector<long>;
00336 particleAvatarConnections[(*ip)->getID()] = aIDs;
00337 }
00338 }
00339
00340 void Store::clear() {
00341 clearAvatars();
00342 inside.clear();
00343
00344 for(std::map<long, Particle*>::iterator iter = particles.begin();
00345 iter != particles.end(); ++iter) {
00346 delete (*iter).second;
00347 }
00348 particles.clear();
00349
00350 clearOutgoing();
00351
00352 if( incoming.size() != 0 ) {
00353 WARN("Incoming list is not empty when Store::clear() is called" << std::endl);
00354 }
00355 incoming.clear();
00356
00357 #ifdef INCL_AVATAR_SEARCH_INCLSort
00358 avatarIterList.clear();
00359 #endif
00360
00361 }
00362
00363 void Store::clearOutgoing() {
00364 for(ParticleIter iter = outgoing.begin(); iter != outgoing.end(); ++iter) {
00365 if((*iter)->isCluster()) {
00366 Cluster *c = dynamic_cast<Cluster *>(*iter);
00367
00368 #ifdef INCLXX_IN_GEANT4_MODE
00369 if(!c)
00370 continue;
00371 #endif
00372 c->deleteParticles();
00373 }
00374 delete (*iter);
00375 }
00376 outgoing.clear();
00377 }
00378
00379 void Store::loadParticles(std::string filename) {
00380 clear();
00381 G4int projectileA, projectileZ, A, Z;
00382 G4double stoppingTime, cutNN;
00383 G4int ID, type, isParticipant;
00384 G4double x, y, z;
00385 G4double px, py, pz, E, v;
00386
00387 std::ifstream in(filename.c_str());
00388 in >> projectileA >> projectileZ >> A >> Z >> stoppingTime >> cutNN;
00389 loadedA = A;
00390 loadedZ = Z;
00391 loadedStoppingTime = stoppingTime;
00392
00393 G4int readA = 0;
00394 G4int readZ = 0;
00395 while(1) {
00396 in >> ID >> type >> isParticipant >> x >> y >> z >> px >> py >> pz >> E >> v;
00397 if(!in.good()) break;
00398 ParticleType t;
00399 if(type == 1) {
00400 t = Proton;
00401 readZ++;
00402 readA++;
00403 }
00404 else if(type == -1) {
00405 t = Neutron;
00406 readA++;
00407 }
00408 else {
00409 FATAL("Unrecognized particle type while loading particles; type=" << type << std::endl);
00410 t = UnknownParticle;
00411 }
00412
00413 Particle *p = new Particle(t, E, ThreeVector(px, py, pz),
00414 ThreeVector(x, y, z));
00415 p->setPotentialEnergy(v);
00416 if(isParticipant == 1) {
00417 p->makeParticipant();
00418 theBook->incrementCascading();
00419 }
00420 add(p);
00421 }
00422
00423 in.close();
00424 }
00425
00426 std::string Store::printParticleConfiguration() {
00427 std::stringstream ss;
00428 G4int A = 0, Z = 0;
00429 for(ParticleIter i = inside.begin(); i != inside.end(); ++i) {
00430 if((*i)->getType() == Proton) {
00431 A++;
00432 Z++;
00433 }
00434 if((*i)->getType() == Neutron) {
00435 A++;
00436 }
00437 }
00438
00439
00440 ss << "0 0 " << A << " " << Z << " "
00441 << "100.0" << " "
00442 << "0.0" << std::endl;
00443
00444 for(ParticleIter i = inside.begin(); i != inside.end(); ++i) {
00445 G4int ID = (*i)->getID();
00446 G4int type = 0;
00447 if((*i)->getType() == Proton) {
00448 type = 1;
00449 }
00450 if((*i)->getType() == Neutron) {
00451 type = -1;
00452 }
00453
00454 G4int isParticipant = 0;
00455 if((*i)->isParticipant()) {
00456 isParticipant = 1;
00457 }
00458
00459 G4double x = (*i)->getPosition().getX();
00460 G4double y = (*i)->getPosition().getY();
00461 G4double z = (*i)->getPosition().getZ();
00462 G4double E = (*i)->getEnergy();
00463 G4double px = (*i)->getMomentum().getX();
00464 G4double py = (*i)->getMomentum().getY();
00465 G4double pz = (*i)->getMomentum().getZ();
00466 G4double V = (*i)->getPotentialEnergy();
00467
00468 ss << ID << " " << type << " " << isParticipant << " "
00469 << x << " " << y << " " << z << " "
00470 << px << " " << py << " " << pz << " "
00471 << E << " " << V << std::endl;
00472 }
00473
00474 return ss.str();
00475 }
00476
00477 void Store::writeParticles(std::string filename) {
00478 std::ofstream out(filename.c_str());
00479 out << printParticleConfiguration();
00480 out.close();
00481 }
00482
00483 std::string Store::printAvatars() {
00484 std::stringstream ss;
00485 IAvatarIter i;
00486 for(i = avatarList.begin(); i != avatarList.end(); ++i) {
00487 ss << (*i)->toString() << std::endl;
00488 }
00489 return ss.str();
00490 }
00491
00492 G4bool Store::containsCollisions() const {
00493 IAvatarIter i;
00494 for(i = avatarList.begin(); i != avatarList.end(); ++i)
00495 if((*i)->getType()==CollisionAvatarType) return true;
00496 return false;
00497 }
00498
00499 }