include/Geant4/G4INCLClusteringModelIntercomparison.hh

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 // INCL++ intra-nuclear cascade model
0027 // Alain Boudard, CEA-Saclay, France
0028 // Joseph Cugnon, University of Liege, Belgium
0029 // Jean-Christophe David, CEA-Saclay, France
0030 // Pekka Kaitaniemi, CEA-Saclay, France, and Helsinki Institute of Physics, Finland
0031 // Sylvie Leray, CEA-Saclay, France
0032 // Davide Mancusi, CEA-Saclay, France
0033 //
0034 #define INCLXX_IN_GEANT4_MODE 1
0035
0036 #include "globals.hh"
0037
0038 #ifndef G4INCLClusteringModelIntercomparison_hh
0039 #define G4INCLClusteringModelIntercomparison_hh 1
0040
0041 #ifdef INCLXX_IN_GEANT4_MODE
0042 #define INCL_CACHING_CLUSTERING_MODEL_INTERCOMPARISON_Set 1
0043 #endif // INCLXX_IN_GEANT4_MODE
0044
0045 #include "G4INCLIClusteringModel.hh"
0046 #include "G4INCLParticle.hh"
0047 #include "G4INCLParticleTable.hh"
0048 #include "G4INCLCluster.hh"
0049 #include "G4INCLNucleus.hh"
0050 #include "G4INCLKinematicsUtils.hh"
0051 #include "G4INCLHashing.hh"
0052
0053 #include <set>
0054 #include <algorithm>
0055
0056 namespace G4INCL {
0057
0058   /** \brief Container for the relevant information
0059    *
0060    * This struct contains all the information that is relevant for the
0061    * clustering algorithm. It is probably more compact than the Particles it
0062    * feeds on, hopefully improving cache performance.
0063    */
0064   struct ConsideredPartner {
0065     Particle *particle;
0066     G4bool isTargetSpectator;
0067     G4int Z;
0068     G4int S;
0069     ThreeVector position;
0070     ThreeVector momentum;
0071     G4double energy;
0072     G4double potentialEnergy;
0073
0074     ConsideredPartner() :
0075       particle(NULL),
0076       isTargetSpectator(false),
0077       Z(0),
0078       S(0),
0079       energy(0.),
0080       potentialEnergy(0.)
0081     {}
0082
0083     ConsideredPartner(Particle * const p) :
0084       particle(p),
0085       isTargetSpectator(particle->isTargetSpectator()),
0086       Z(particle->getZ()),
0087       S(particle->getS()),
0088       position(particle->getPosition()),
0089       momentum(particle->getMomentum()),
0090       energy(particle->getEnergy()),
0091       potentialEnergy(particle->getPotentialEnergy())
0092     {}
0093   };
0094
0095   /// \brief Cluster coalescence algorithm used in the IAEA intercomparison
0096   class ClusteringModelIntercomparison : public IClusteringModel {
0097   public:
0098     ClusteringModelIntercomparison(Config const * const theConfig) :
0099       theNucleus(NULL),
0100       selectedA(0),
0101       selectedZ(0),
0102       selectedS(0),
0103       sqtot(0.),
0104       cascadingEnergyPool(0.),
0105       protonMass(ParticleTable::getRealMass(Proton)),
0106       neutronMass(ParticleTable::getRealMass(Neutron)),
0107       lambdaMass(ParticleTable::getRealMass(Lambda)),
0108       runningMaxClusterAlgorithmMass(theConfig->getClusterMaxMass()),
0109       nConsideredMax(0),
0110       nConsidered(0),
0111       consideredPartners(NULL),
0112       isInRunningConfiguration(NULL),
0113       maxMassConfigurationSkipping(ParticleTable::maxClusterMass)
0114     {
0115       // Set up the maximum charge and neutron number for clusters
0116       clusterZMaxAll = 0;
0117       clusterNMaxAll = 0;
0118       for(G4int A=0; A<=runningMaxClusterAlgorithmMass; ++A) {
0119         if(clusterZMax[A]>clusterZMaxAll)
0120           clusterZMaxAll = clusterZMax[A];
0121         if(A-clusterZMin[A]>clusterNMaxAll)
0122           clusterNMaxAll = A-clusterZMin[A];
0123       }
0124       std::fill(candidateConfiguration,
0125                 candidateConfiguration + ParticleTable::maxClusterMass,
0126                 static_cast<Particle*>(NULL));
0127
0128       std::fill(runningEnergies,
0129                 runningEnergies + ParticleTable::maxClusterMass,
0130                 0.0);
0131
0132       std::fill(runningPotentials,
0133                 runningPotentials + ParticleTable::maxClusterMass,
0134                 0.0);
0135
0136       std::fill(runningConfiguration,
0137                 runningConfiguration + ParticleTable::maxClusterMass,
0138                 -1);
0139
0140     }
0141
0142     virtual ~ClusteringModelIntercomparison() {
0143       delete [] consideredPartners;
0144       delete [] isInRunningConfiguration;
0145     }
0146
0147     virtual Cluster* getCluster(Nucleus*, Particle*);
0148     virtual G4bool clusterCanEscape(Nucleus const * const, Cluster const * const);
0149
0150   private:
0151     void findClusterStartingFrom(const G4int oldA, const G4int oldZ, const G4int oldS);
0152     G4double getPhaseSpace(const G4int oldA, ConsideredPartner const &p);
0153
0154     Nucleus *theNucleus;
0155
0156     G4double runningEnergies[ParticleTable::maxClusterMass+1];
0157     ThreeVector runningMomenta[ParticleTable::maxClusterMass+1];
0158     ThreeVector runningPositions[ParticleTable::maxClusterMass+1];
0159     G4double runningPotentials[ParticleTable::maxClusterMass+1];
0160 #if defined(INCL_CACHING_CLUSTERING_MODEL_INTERCOMPARISON_HashMask)
0161     Hashing::NucleonItem runningConfiguration[ParticleTable::maxClusterMass];
0162 #elif defined(INCL_CACHING_CLUSTERING_MODEL_INTERCOMPARISON_Set) || defined(INCL_CACHING_CLUSTERING_MODEL_INTERCOMPARISON_None)
0163     G4int runningConfiguration[ParticleTable::maxClusterMass];
0164 #else
0165 #error Unrecognized INCL_CACHING_CLUSTERING_MODEL_INTERCOMPARISON. Allowed values are: Set, HashMask, None.
0166 #endif
0167
0168     G4int selectedA, selectedZ, selectedS;
0169     G4double sqtot;
0170
0171     G4int clusterZMaxAll, clusterNMaxAll;
0172
0173     G4double cascadingEnergyPool;
0174
0175     /// \brief Lower limit of Z for cluster of mass A
0176     static const G4int clusterZMin[ParticleTable::maxClusterMass+1];
0177     /// \brief Upper limit of Z for cluster of mass A
0178     static const G4int clusterZMax[ParticleTable::maxClusterMass+1];
0179
0180     /// \brief Precomputed factor 1.0/A
0181     static const G4double clusterPosFact[ParticleTable::maxClusterMass+1];
0182
0183     /// \brief Precomputed factor (1.0/A)^2
0184     static const G4double clusterPosFact2[ParticleTable::maxClusterMass+1];
0185
0186     /// \brief Phase-space parameters for cluster formation
0187     static const G4double clusterPhaseSpaceCut[ParticleTable::maxClusterMass+1];
0188
0189     static const G4double limitCosEscapeAngle;
0190
0191     const G4double protonMass;
0192     const G4double neutronMass;
0193     const G4double lambdaMass;
0194
0195     G4int runningMaxClusterAlgorithmMass;
0196
0197     G4int nConsideredMax;
0198     G4int nConsidered;
0199
0200     /** \brief Array of considered cluster partners
0201      *
0202      * A dynamical array of ConsideredPartner objects is allocated on this
0203      * variable and filled with pointers to nucleons which are eligible for
0204      * clustering. We used to use a ParticleList for this purpose, but this
0205      * made it very cumbersome to check whether nucleons had already been
0206      * included in the running configuration. Using an array of Particle*
0207      * coupled with a boolean mask (\see{isInRunningConfiguration}) reduces the
0208      * overhead by a large amount.  Running times for 1-GeV p+Pb208 went down
0209      * by almost 30% (!).
0210      *
0211      * Lesson learnt: when you need speed, nothing beats a good ol' array.
0212      */
0213     ConsideredPartner *consideredPartners;
0214
0215     /** \brief Array of flags for nucleons in the running configuration
0216      *
0217      * Clustering partners that are already used in the running cluster
0218      * configuration are flagged as "true" in this array.
0219      */
0220     G4bool *isInRunningConfiguration;
0221
0222     /** \brief Best cluster configuration
0223      *
0224      * This array contains pointers to the nucleons which make up the best
0225      * cluster configuration that has been found so far.
0226      */
0227     Particle *candidateConfiguration[ParticleTable::maxClusterMass];
0228
0229 #if defined(INCL_CACHING_CLUSTERING_MODEL_INTERCOMPARISON_HashMask)
0230     typedef std::set<Hashing::HashType> HashContainer;
0231     typedef HashContainer::iterator HashIterator;
0232
0233     /// \brief Array of containers for configurations that have already been checked
0234     HashContainer checkedConfigurations[ParticleTable::maxClusterMass-2];
0235 #elif defined(INCL_CACHING_CLUSTERING_MODEL_INTERCOMPARISON_Set)
0236     /** \brief Class for storing and comparing sorted nucleon configurations
0237      *
0238      * This class is actually just a wrapper around an array of Particle*
0239      * pointers. It provides a lexicographical comparison operator
0240      * (SortedNucleonConfiguration::operator<) for inclusion in std::set
0241      * containers.
0242      */
0243     class SortedNucleonConfiguration {
0244       public:
0245         // Use Particle* as nucleon identifiers
0246         typedef G4int NucleonItem;
0247
0248         /// \brief Constructor
0249         SortedNucleonConfiguration() : theSize(0), nucleons(NULL) {}
0250
0251         /// \brief Copy constructor
0252         SortedNucleonConfiguration(const SortedNucleonConfiguration &rhs) :
0253           theSize(rhs.theSize),
0254           nucleons(new NucleonItem[theSize])
0255       {
0256         std::copy(rhs.nucleons, rhs.nucleons+theSize, nucleons);
0257       }
0258
0259         /// \brief Destructor
0260         ~SortedNucleonConfiguration() {
0261           delete [] nucleons;
0262         }
0263
0264         /// \brief Helper method for the assignment operator
0265         void swap(SortedNucleonConfiguration &rhs) {
0266           std::swap(theSize, rhs.theSize);
0267           std::swap(nucleons, rhs.nucleons);
0268         }
0269
0270         /// \brief Assignment operator
0271         SortedNucleonConfiguration &operator=(const SortedNucleonConfiguration &rhs) {
0272           SortedNucleonConfiguration tempConfig(rhs);
0273           swap(tempConfig);
0274           return *this;
0275         }
0276
0277         /** \brief Order operator for SortedNucleonConfiguration
0278          *
0279          * The comparison is done lexicographically (i.e. from the first
0280          * element to the last).
0281          */
0282         G4bool operator<(const SortedNucleonConfiguration &rhs) const {
0283 // assert(theSize==rhs.theSize);
0284           return std::lexicographical_compare(nucleons, nucleons+theSize, rhs.nucleons, rhs.nucleons+theSize);
0285         }
0286
0287         /// \brief Fill configuration with array of NucleonItem
0288         void fill(NucleonItem *config, size_t n) {
0289           theSize = n;
0290           nucleons = new NucleonItem[theSize];
0291           std::copy(config, config+theSize, nucleons);
0292           std::sort(nucleons, nucleons+theSize);
0293         }
0294
0295       private:
0296         /// \brief Size of the array
0297         size_t theSize;
0298
0299         /// \brief The real array
0300         NucleonItem *nucleons;
0301     };
0302
0303     typedef std::set<SortedNucleonConfiguration> SortedNucleonConfigurationContainer;
0304     typedef SortedNucleonConfigurationContainer::iterator SortedNucleonConfigurationIterator;
0305
0306     /// \brief Array of containers for configurations that have already been checked
0307     SortedNucleonConfigurationContainer checkedConfigurations[ParticleTable::maxClusterMass-2];
0308 #elif !defined(INCL_CACHING_CLUSTERING_MODEL_INTERCOMPARISON_None)
0309 #error Unrecognized INCL_CACHING_CLUSTERING_MODEL_INTERCOMPARISON. Allowed values are: Set, HashMask, None.
0310 #endif
0311
0312     /** \brief Maximum mass for configuration storage
0313      *
0314      * Skipping configurations becomes inefficient above this mass.
0315      */
0316     G4int maxMassConfigurationSkipping;
0317   };
0318
0319 }
0320
0321 #endif