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 #include "G4StatMFMicroManager.hh"
00034 #include "G4HadronicException.hh"
00035
00036
00037
00038 G4StatMFMicroManager::G4StatMFMicroManager(const G4StatMFMicroManager & )
00039 {
00040 throw G4HadronicException(__FILE__, __LINE__, "G4StatMFMicroManager::copy_constructor meant to not be accessable");
00041 }
00042
00043
00044
00045 G4StatMFMicroManager & G4StatMFMicroManager::
00046 operator=(const G4StatMFMicroManager & )
00047 {
00048 throw G4HadronicException(__FILE__, __LINE__, "G4StatMFMicroManager::operator= meant to not be accessable");
00049 return *this;
00050 }
00051
00052
00053 G4bool G4StatMFMicroManager::operator==(const G4StatMFMicroManager & ) const
00054 {
00055 return false;
00056 }
00057
00058
00059 G4bool G4StatMFMicroManager::operator!=(const G4StatMFMicroManager & ) const
00060 {
00061 return true;
00062 }
00063
00064
00065
00066
00067 G4StatMFMicroManager::G4StatMFMicroManager(const G4Fragment & theFragment, const G4int multiplicity,
00068 const G4double FreeIntE, const G4double SCompNuc) :
00069 _Normalization(0.0)
00070 {
00071
00072 Initialize(theFragment,multiplicity,FreeIntE,SCompNuc);
00073 }
00074
00075
00076
00077 G4StatMFMicroManager::~G4StatMFMicroManager()
00078 {
00079 if (!_Partition.empty())
00080 {
00081 std::for_each(_Partition.begin(),_Partition.end(),
00082 DeleteFragment());
00083 }
00084 }
00085
00086
00087
00088
00089
00090 void G4StatMFMicroManager::Initialize(const G4Fragment & theFragment, const G4int im,
00091 const G4double FreeIntE, const G4double SCompNuc)
00092 {
00093 G4int i;
00094
00095 G4double U = theFragment.GetExcitationEnergy();
00096
00097 G4double A = theFragment.GetA();
00098 G4double Z = theFragment.GetZ();
00099
00100
00101 _WW = 0.0;
00102
00103
00104 _MeanMultiplicity = 0.0;
00105
00106
00107 _MeanTemperature = 0.0;
00108
00109
00110 _MeanEntropy = 0.0;
00111
00112
00113
00114
00115 G4int FragmentAtomicNumbers[4];
00116
00117
00118
00119
00120
00121 FragmentAtomicNumbers[im-1] = static_cast<G4int>(A);
00122 for (i = 0; i < (im - 1); i++) FragmentAtomicNumbers[i] = 0;
00123
00124
00125
00126
00127 while (MakePartition(im,FragmentAtomicNumbers)) {
00128
00129
00130 G4StatMFMicroPartition * aPartition = new G4StatMFMicroPartition(static_cast<G4int>(A),
00131 static_cast<G4int>(Z));
00132 G4double PartitionProbability = 0.0;
00133
00134 for (i = im-1; i >= 0; i--) aPartition->SetPartitionFragment(FragmentAtomicNumbers[i]);
00135 PartitionProbability = aPartition->CalcPartitionProbability(U,FreeIntE,SCompNuc);
00136 _Partition.push_back(aPartition);
00137
00138 _WW += PartitionProbability;
00139 _MeanMultiplicity += im*PartitionProbability;
00140 _MeanTemperature += aPartition->GetTemperature() * PartitionProbability;
00141 if (PartitionProbability > 0.0)
00142 _MeanEntropy += PartitionProbability * aPartition->GetEntropy();
00143
00144 }
00145
00146
00147
00148
00149
00150 }
00151
00152
00153 G4bool G4StatMFMicroManager::MakePartition(const G4int k, G4int * ANumbers)
00154
00155
00156
00157 {
00158 G4int l = 1;
00159 while (l < k) {
00160 G4int tmp = ANumbers[l-1] + ANumbers[k-1];
00161 ANumbers[l-1] += 1;
00162 ANumbers[k-1] -= 1;
00163 if (ANumbers[l-1] > ANumbers[l] || ANumbers[k-2] > ANumbers[k-1]) {
00164 ANumbers[l-1] = 1;
00165 ANumbers[k-1] = tmp - 1;
00166 l++;
00167 } else return true;
00168 }
00169 return false;
00170 }
00171
00172
00173
00174 void G4StatMFMicroManager::Normalize(const G4double Norm)
00175 {
00176 _Normalization = Norm;
00177 _WW /= Norm;
00178 _MeanMultiplicity /= Norm;
00179 _MeanTemperature /= Norm;
00180 _MeanEntropy /= Norm;
00181
00182 return;
00183 }
00184
00185 G4StatMFChannel * G4StatMFMicroManager::ChooseChannel(const G4double A0, const G4double Z0,
00186 const G4double MeanT)
00187 {
00188 G4double RandNumber = _Normalization * _WW * G4UniformRand();
00189 G4double AccumWeight = 0.0;
00190
00191 for (std::vector<G4StatMFMicroPartition*>::iterator i = _Partition.begin();
00192 i != _Partition.end(); ++i)
00193 {
00194 AccumWeight += (*i)->GetProbability();
00195 if (RandNumber < AccumWeight)
00196 return (*i)->ChooseZ(A0,Z0,MeanT);
00197 }
00198
00199 throw G4HadronicException(__FILE__, __LINE__,
00200 "G4StatMFMicroCanonical::ChooseChannel: Couldn't find a channel.");
00201 return 0;
00202 }