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
00034
00035
00036
00037
00038 #include "G4ios.hh"
00039
00040 #include "G4PiMinusAbsorptionAtRest.hh"
00041
00042 #include "G4SystemOfUnits.hh"
00043 #include "G4PiMinusStopLi.hh"
00044 #include "G4PiMinusStopC.hh"
00045 #include "G4PiMinusStopN.hh"
00046 #include "G4PiMinusStopO.hh"
00047 #include "G4PiMinusStopAl.hh"
00048 #include "G4PiMinusStopCu.hh"
00049 #include "G4PiMinusStopCo.hh"
00050 #include "G4PiMinusStopTa.hh"
00051 #include "G4PiMinusStopPb.hh"
00052 #include "G4StopTheoDeexcitation.hh"
00053 #include "G4StopDummyDeexcitation.hh"
00054 #include "G4DynamicParticle.hh"
00055 #include "G4DynamicParticleVector.hh"
00056 #include "Randomize.hh"
00057 #include "G4ThreeVector.hh"
00058 #include "G4LorentzVector.hh"
00059 #include "G4HadronicProcessStore.hh"
00060 #include "G4HadronicDeprecate.hh"
00061
00062
00063
00064
00065 G4PiMinusAbsorptionAtRest::G4PiMinusAbsorptionAtRest(const G4String& processName,
00066 G4ProcessType aType) :
00067 G4VRestProcess (processName, aType)
00068 {
00069 G4HadronicDeprecate("G4PiMinusAbsorptionAtRest");
00070
00071 SetProcessSubType(fHadronAtRest);
00072
00073 _indexDeexcitation = 0;
00074
00075 if (verboseLevel>0)
00076 { G4cout << GetProcessName() << " is created "<< G4endl; }
00077 G4HadronicProcessStore::Instance()->RegisterExtraProcess(this);
00078 }
00079
00080
00081
00082
00083 G4PiMinusAbsorptionAtRest::~G4PiMinusAbsorptionAtRest()
00084 {
00085 G4HadronicProcessStore::Instance()->DeRegisterExtraProcess(this);
00086 }
00087
00088 void G4PiMinusAbsorptionAtRest::PreparePhysicsTable(const G4ParticleDefinition& p)
00089 {
00090 G4HadronicProcessStore::Instance()->RegisterParticleForExtraProcess(this, &p);
00091 }
00092
00093 void G4PiMinusAbsorptionAtRest::BuildPhysicsTable(const G4ParticleDefinition& p)
00094 {
00095 G4HadronicProcessStore::Instance()->PrintInfo(&p);
00096 }
00097
00098 G4VParticleChange* G4PiMinusAbsorptionAtRest::AtRestDoIt(const G4Track& track, const G4Step& )
00099 {
00100 const G4DynamicParticle* stoppedHadron = track.GetDynamicParticle();
00101
00102
00103 if (! IsApplicable(*(stoppedHadron->GetDefinition())))
00104 {
00105 G4cerr
00106 << "G4PiMinusAbsorptionAtRest: ERROR, particle must be a pion minus!"
00107 << G4endl;
00108 return NULL;
00109 }
00110
00111
00112 const G4Material* material = track.GetMaterial();
00113
00114 G4double A=-1;
00115 G4double Z=-1;
00116 G4double random = G4UniformRand();
00117 const G4ElementVector* theElementVector = material->GetElementVector();
00118 unsigned int i;
00119 G4double sum = 0;
00120 G4double totalsum=0;
00121 for(i=0; i<material->GetNumberOfElements(); ++i)
00122 {
00123 if((*theElementVector)[i]->GetZ()!=1) totalsum+=material->GetFractionVector()[i];
00124 }
00125 for (i = 0; i<material->GetNumberOfElements(); ++i)
00126 {
00127 if((*theElementVector)[i]->GetZ()!=1) sum += material->GetFractionVector()[i];
00128 if ( sum/totalsum > random )
00129 {
00130 A = (*theElementVector)[i]->GetA()*mole/g;
00131 Z = (*theElementVector)[i]->GetZ();
00132 break;
00133 }
00134 }
00135
00136
00137
00138 G4PiMinusStopMaterial* algorithm = LoadAlgorithm(static_cast<G4int>(Z));
00139 G4PiMinusStopAbsorption stopAbsorption(algorithm,Z,A);
00140 stopAbsorption.SetVerboseLevel(verboseLevel);
00141
00142 G4DynamicParticleVector* absorptionProducts = stopAbsorption.DoAbsorption();
00143
00144
00145
00146 G4double pionEnergy = stoppedHadron->GetTotalEnergy();
00147 G4double excitation = pionEnergy - stopAbsorption.Energy();
00148 if (excitation < 0.)
00149 {
00150 G4Exception("G4PiMinusAbsorptionAtRest::AtRestDoIt()", "HAD_STOP_0000",
00151 FatalException, "Excitation energy < 0");
00152 }
00153 if (verboseLevel>0) { G4cout << " excitation " << excitation << G4endl; }
00154
00155 G4StopDeexcitationAlgorithm* nucleusAlgorithm = LoadNucleusAlgorithm();
00156 G4StopDeexcitation stopDeexcitation(nucleusAlgorithm);
00157
00158 G4double newZ = Z - stopAbsorption.NProtons();
00159 G4double newN = A - Z - stopAbsorption.NNeutrons();
00160 G4double newA = newZ + newN;
00161 G4ReactionProductVector* fragmentationProducts = stopDeexcitation.DoBreakUp(newA,newZ,excitation,stopAbsorption.RecoilMomentum());
00162
00163 unsigned int nAbsorptionProducts = 0;
00164 if (absorptionProducts != 0)
00165 { nAbsorptionProducts = absorptionProducts->size(); }
00166
00167 unsigned int nFragmentationProducts = 0;
00168 if (fragmentationProducts != 0)
00169 { nFragmentationProducts = fragmentationProducts->size(); }
00170
00171 if (verboseLevel>0)
00172 {
00173 G4cout << "nAbsorptionProducts = " << nAbsorptionProducts
00174 << " nFragmentationProducts = " << nFragmentationProducts
00175 << G4endl;
00176 }
00177
00178
00179
00180 aParticleChange.Initialize(track);
00181 aParticleChange.SetNumberOfSecondaries(G4int(nAbsorptionProducts + nFragmentationProducts));
00182
00183 for (i = 0; i<nAbsorptionProducts; i++)
00184 { aParticleChange.AddSecondary((*absorptionProducts)[i]); }
00185
00186
00187
00188 for(i=0; i<nFragmentationProducts; i++)
00189 {
00190 G4DynamicParticle * aNew =
00191 new G4DynamicParticle((*fragmentationProducts)[i]->GetDefinition(),
00192 (*fragmentationProducts)[i]->GetMomentum());
00193 G4double newTime = aParticleChange.GetGlobalTime((*fragmentationProducts)[i]->GetFormationTime());
00194 aParticleChange.AddSecondary(aNew, newTime);
00195 delete (*fragmentationProducts)[i];
00196 }
00197
00198 if (fragmentationProducts != 0) delete fragmentationProducts;
00199
00200 if (_indexDeexcitation == 1) aParticleChange.ProposeLocalEnergyDeposit(excitation);
00201
00202
00203 aParticleChange.ProposeTrackStatus(fStopAndKill);
00204
00205 return &aParticleChange;
00206
00207 }
00208
00209 G4PiMinusStopMaterial* G4PiMinusAbsorptionAtRest::LoadAlgorithm(int Z)
00210 {
00211 if (verboseLevel>0)
00212 {
00213 G4cout << "Load material algorithm " << Z << G4endl;
00214 }
00215
00216 G4int index = 0;
00217 if (Z > 0 && Z < 4) {index = 3;}
00218 if (Z > 3 && Z < 7) {index = 6;}
00219 if (Z == 7) {index = 7;}
00220 if (Z >= 8 && Z<= 11) {index = 8;}
00221 if (Z >= 12 && Z<= 18) {index = 13;}
00222 if (Z >=19 && Z<= 27) {index = 27;}
00223 if (Z >= 28 && Z<= 51) {index = 29;}
00224 if (Z >=52 ) {index = 73;}
00225
00226 switch (index)
00227 {
00228 case 3:
00229 if (verboseLevel>0)
00230 { G4cout << " =================== Load Li algorithm " << G4endl; }
00231 return new G4PiMinusStopLi();
00232 case 6:
00233 if (verboseLevel>0)
00234 { G4cout << " =================== Load C algorithm " << G4endl; }
00235 return new G4PiMinusStopC();
00236 case 7:
00237 if (verboseLevel>0)
00238 { G4cout << " =================== Load N algorithm " << G4endl; }
00239 return new G4PiMinusStopN();
00240 case 8:
00241 if (verboseLevel>0)
00242 { G4cout << " =================== Load O algorithm " << G4endl; }
00243 return new G4PiMinusStopO();
00244 case 13:
00245 if (verboseLevel>0)
00246 { G4cout << " =================== Load Al algorithm " << G4endl; }
00247 return new G4PiMinusStopAl();
00248 case 27:
00249 if (verboseLevel>0)
00250 { G4cout << " =================== Load Cu algorithm " << G4endl; }
00251 return new G4PiMinusStopCu();
00252 case 29:
00253 if (verboseLevel>0)
00254 { G4cout << " =================== Load Co algorithm " << G4endl; }
00255 return new G4PiMinusStopCo();
00256 case 73:
00257 if (verboseLevel>0)
00258 { G4cout << " =================== Load Ta algorithm " << G4endl; }
00259 return new G4PiMinusStopTa();
00260 default:
00261 if (verboseLevel>0)
00262 { G4cout << " =================== Load default material algorithm " << G4endl; }
00263 return new G4PiMinusStopC();
00264 }
00265 }
00266
00267 G4StopDeexcitationAlgorithm* G4PiMinusAbsorptionAtRest::LoadNucleusAlgorithm()
00268 {
00269
00270 switch (_indexDeexcitation)
00271 {
00272 case 0:
00273 if (verboseLevel>0)
00274 { G4cout << " =================== Load Theo deexcitation " << G4endl; }
00275 return new G4StopTheoDeexcitation();
00276 case 1:
00277 if (verboseLevel>0)
00278 { G4cout << " =================== Load Dummy deexcitation " << G4endl; }
00279 return new G4StopDummyDeexcitation();
00280 default:
00281 if (verboseLevel>0)
00282 { G4cout << " =================== Load default deexcitation " << G4endl; }
00283 return new G4StopTheoDeexcitation();
00284 }
00285 }
00286
00287 void G4PiMinusAbsorptionAtRest::SetDeexcitationAlgorithm(G4int index)
00288 {
00289 _indexDeexcitation = index;
00290 }