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
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00059
00060 #include "G4Event.hh"
00061 #include "Randomize.hh"
00062 #include "G4GeneralParticleSource.hh"
00063
00064 G4GeneralParticleSource::G4GeneralParticleSource()
00065 : multiple_vertex(false), flat_sampling(false)
00066 {
00067 sourceVector.clear();
00068 sourceIntensity.clear();
00069 sourceProbability.clear();
00070 currentSource = new G4SingleParticleSource();
00071 sourceVector.push_back(currentSource);
00072 sourceIntensity.push_back(1.);
00073 currentSourceIdx = G4int(sourceVector.size() - 1);
00074 theMessenger = new G4GeneralParticleSourceMessenger(this);
00075 theMessenger->SetParticleGun(currentSource);
00076 IntensityNormalization();
00077 }
00078
00079 G4GeneralParticleSource::~G4GeneralParticleSource()
00080 {
00081 delete theMessenger;
00082 }
00083
00084 void G4GeneralParticleSource::AddaSource(G4double aV)
00085 {
00086 currentSource = new G4SingleParticleSource();
00087 theMessenger->SetParticleGun(currentSource);
00088 sourceVector.push_back(currentSource);
00089 sourceIntensity.push_back(aV);
00090 currentSourceIdx = G4int(sourceVector.size() - 1);
00091 IntensityNormalization();
00092 }
00093
00094 void G4GeneralParticleSource::IntensityNormalization()
00095 {
00096 G4double total = 0.;
00097 size_t i = 0 ;
00098 for (i = 0; i < sourceIntensity.size(); i++)
00099 total += sourceIntensity[i] ;
00100
00101 sourceProbability.clear();
00102 std::vector <G4double> sourceNormalizedIntensity;
00103 sourceNormalizedIntensity.clear();
00104
00105 sourceNormalizedIntensity.push_back(sourceIntensity[0]/total);
00106 sourceProbability.push_back(sourceNormalizedIntensity[0]);
00107
00108 for ( i = 1 ; i < sourceIntensity.size(); i++) {
00109 sourceNormalizedIntensity.push_back(sourceIntensity[i]/total);
00110 sourceProbability.push_back(sourceNormalizedIntensity[i] + sourceProbability[i-1]);
00111 }
00112
00113
00114 for ( i = 0 ; i < sourceIntensity.size(); i++) {
00115 if (!flat_sampling) {
00116 sourceVector[i]->GetBiasRndm()->SetIntensityWeight(1.);
00117 } else {
00118 sourceVector[i]->GetBiasRndm()->SetIntensityWeight(sourceNormalizedIntensity[i]*sourceIntensity.size());
00119 }
00120 }
00121
00122 normalised = true;
00123 }
00124
00125 void G4GeneralParticleSource::ListSource()
00126 {
00127 G4cout << " The number of particle sources is " << sourceIntensity.size() << G4endl;
00128 for (size_t i = 0 ; i < sourceIntensity.size(); i++)
00129 G4cout << " source " << i << " intensity is " << sourceIntensity[i] << G4endl;
00130 }
00131
00132 void G4GeneralParticleSource::SetCurrentSourceto(G4int aV)
00133 {
00134 size_t id = size_t (aV) ;
00135 if ( id <= sourceIntensity.size() ) {
00136 currentSourceIdx = aV;
00137 currentSource = sourceVector[id];
00138 theMessenger->SetParticleGun(currentSource);
00139
00140 } else {
00141 G4cout << " source index is invalid " << G4endl;
00142 G4cout << " it shall be <= " << sourceIntensity.size() << G4endl;
00143 }
00144 }
00145
00146 void G4GeneralParticleSource::SetCurrentSourceIntensity(G4double aV)
00147 {
00148 sourceIntensity[currentSourceIdx] = aV;
00149 normalised = false;
00150 }
00151
00152 void G4GeneralParticleSource::ClearAll()
00153 {
00154 currentSourceIdx = -1;
00155 currentSource = 0;
00156 sourceVector.clear();
00157 sourceIntensity.clear();
00158 sourceProbability.clear();
00159 }
00160
00161 void G4GeneralParticleSource::DeleteaSource(G4int aV)
00162 {
00163 size_t id = size_t (aV) ;
00164 if ( id <= sourceIntensity.size() ) {
00165 sourceVector.erase(sourceVector.begin()+aV);
00166 sourceIntensity.erase(sourceIntensity.begin()+aV);
00167 normalised = false ;
00168 if (currentSourceIdx == aV ) {
00169 if ( sourceIntensity.size() > 0 ) {
00170 currentSource = sourceVector[0];
00171 currentSourceIdx = 1;
00172 } else {
00173 currentSource = 0;
00174 currentSourceIdx = -1;
00175 }
00176 }
00177 } else {
00178 G4cout << " source index is invalid " << G4endl;
00179 G4cout << " it shall be <= " << sourceIntensity.size() << G4endl;
00180 }
00181 }
00182
00183 void G4GeneralParticleSource::GeneratePrimaryVertex(G4Event* evt)
00184 {
00185 if (!multiple_vertex){
00186 if (sourceIntensity.size() > 1) {
00187 if (!normalised) IntensityNormalization();
00188 G4double rndm = G4UniformRand();
00189 size_t i = 0 ;
00190 if (!flat_sampling) {
00191 while ( rndm > sourceProbability[i] ) i++;
00192 (currentSource = sourceVector[i]);
00193 } else {
00194 i = size_t (sourceIntensity.size()*rndm);
00195 currentSource = sourceVector[i];
00196 }
00197 }
00198 currentSource-> GeneratePrimaryVertex(evt);
00199 }
00200 else {
00201 for (size_t i = 0; i < sourceIntensity.size(); i++) {
00202 sourceVector[i]->GeneratePrimaryVertex(evt);
00203 }
00204 }
00205 }