| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 26
27 #include "G4LENDCapture.hh" 27 #include "G4LENDCapture.hh"
28 #include "G4Fragment.hh" 28 #include "G4Fragment.hh"
29 #include "G4PhotonEvaporation.hh" 29 #include "G4PhotonEvaporation.hh"
30 #include "G4SystemOfUnits.hh" 30 #include "G4SystemOfUnits.hh"
31 #include "G4Nucleus.hh" 31 #include "G4Nucleus.hh"
32 #include "G4ParticleTable.hh" 32 #include "G4ParticleTable.hh"
33 #include "G4IonTable.hh" 33 #include "G4IonTable.hh"
34 34
35 G4HadFinalState * G4LENDCapture::ApplyYourself 35 G4HadFinalState * G4LENDCapture::ApplyYourself(const G4HadProjectile& aTrack, G4Nucleus& aTarg )
36 { 36 {
37 37
38 G4double temp = aTrack.GetMaterial()->GetTe 38 G4double temp = aTrack.GetMaterial()->GetTemperature();
39 39
40 //G4int iZ = int ( aTarg.GetZ() ); 40 //G4int iZ = int ( aTarg.GetZ() );
41 //G4int iA = int ( aTarg.GetN() ); 41 //G4int iA = int ( aTarg.GetN() );
42 //migrate to integer A and Z (GetN_asInt re 42 //migrate to integer A and Z (GetN_asInt returns number of neutrons in the nucleus since this)
43 G4int iZ = aTarg.GetZ_asInt(); 43 G4int iZ = aTarg.GetZ_asInt();
44 G4int iA = aTarg.GetA_asInt(); 44 G4int iA = aTarg.GetA_asInt();
45 G4int iM = 0; 45 G4int iM = 0;
46 if ( aTarg.GetIsotope() != NULL ) { 46 if ( aTarg.GetIsotope() != NULL ) {
47 iM = aTarg.GetIsotope()->Getm(); 47 iM = aTarg.GetIsotope()->Getm();
48 } 48 }
49 49
50 G4double ke = aTrack.GetKineticEnergy(); 50 G4double ke = aTrack.GetKineticEnergy();
51 51
52 G4HadFinalState* theResult = &theParticleCh 52 G4HadFinalState* theResult = &theParticleChange;
53 theResult->Clear(); 53 theResult->Clear();
54 54
55 G4GIDI_target* aTarget = get_target_from_ma 55 G4GIDI_target* aTarget = get_target_from_map( lend_manager->GetNucleusEncoding( iZ , iA , iM ) );
56 if ( aTarget == NULL ) return returnUnchang 56 if ( aTarget == NULL ) return returnUnchanged( aTrack , theResult );
57 std::vector<G4GIDI_Product>* products = aTa 57 std::vector<G4GIDI_Product>* products = aTarget->getCaptureFinalState( ke*MeV, temp, MyRNG, NULL );
58 58
59 G4int ipZ = aTrack.GetDefinition()->GetAtom 59 G4int ipZ = aTrack.GetDefinition()->GetAtomicNumber();
60 G4int ipA = aTrack.GetDefinition()->GetAtom 60 G4int ipA = aTrack.GetDefinition()->GetAtomicMass();
61 61
62 G4bool needResidual=true; 62 G4bool needResidual=true;
63 63
64 G4ThreeVector p(0,0,0); 64 G4ThreeVector p(0,0,0);
65 if ( products != NULL ) 65 if ( products != NULL )
66 { 66 {
67 67
68 G4int totN = 0; 68 G4int totN = 0;
69 69
70 for ( G4int j = 0; j < int( products->si 70 for ( G4int j = 0; j < int( products->size() ); j++ )
71 { 71 {
72 G4int jZ = (*products)[j].Z; 72 G4int jZ = (*products)[j].Z;
73 G4int jA = (*products)[j].A; 73 G4int jA = (*products)[j].A;
74 74
75 //G4cout << "ZA = " << 1000 * (*produ 75 //G4cout << "ZA = " << 1000 * (*products)[j].Z + (*products)[j].A << " EK = "
76 // << (*products)[j].kineticEnerg 76 // << (*products)[j].kineticEnergy
77 // << " px " << (*products)[j]. 77 // << " px " << (*products)[j].px
78 // << " py " << (*products)[j]. 78 // << " py " << (*products)[j].py
79 // << " pz " << (*products)[j]. 79 // << " pz " << (*products)[j].pz
80 // << G4endl; 80 // << G4endl;
81 81
82 if ( jZ == iZ + ipZ && jA == iA + ipA 82 if ( jZ == iZ + ipZ && jA == iA + ipA ) needResidual = false;
83 83
84 G4ThreeVector dp((*products)[j].px,(* 84 G4ThreeVector dp((*products)[j].px,(*products)[j].py,(*products)[j].pz);
85 p += dp; 85 p += dp;
86 86
87 G4DynamicParticle* theSec = new G4Dyn 87 G4DynamicParticle* theSec = new G4DynamicParticle;
88 88
89 if ( jA == 1 && jZ == 1 ) { 89 if ( jA == 1 && jZ == 1 ) {
90 theSec->SetDefinition( G4Proton::P 90 theSec->SetDefinition( G4Proton::Proton() );
91 totN += 1; 91 totN += 1;
92 } 92 }
93 else if ( jA == 1 && jZ == 0 ) 93 else if ( jA == 1 && jZ == 0 )
94 { 94 {
95 theSec->SetDefinition( G4Neutron:: 95 theSec->SetDefinition( G4Neutron::Neutron() );
96 totN += 1; 96 totN += 1;
97 } 97 }
98 else if ( jZ > 0 ) { 98 else if ( jZ > 0 ) {
99 if ( jA != 0 ) 99 if ( jA != 0 )
100 { 100 {
101 theSec->SetDefinition( G4IonTab 101 theSec->SetDefinition( G4IonTable::GetIonTable()->GetIon( jZ , jA , iM ) );
102 totN += jA; 102 totN += jA;
103 } 103 }
104 else 104 else
105 { 105 {
106 theSec->SetDefinition( G4IonTab 106 theSec->SetDefinition( G4IonTable::GetIonTable()->GetIon( jZ , iA+1-totN , iM ) );
107 } 107 }
108 } 108 }
109 else { 109 else {
110 theSec->SetDefinition( G4Gamma::Ga 110 theSec->SetDefinition( G4Gamma::Gamma() );
111 } 111 }
112 112
113 theSec->SetMomentum( G4ThreeVector( ( 113 theSec->SetMomentum( G4ThreeVector( (*products)[j].px*MeV , (*products)[j].py*MeV , (*products)[j].pz*MeV ) );
114 114
115 /* 115 /*
116 if ( dp.mag() == 0 ) 116 if ( dp.mag() == 0 )
117 { 117 {
118 //theSec->SetMomentum( -p*MeV ); 118 //theSec->SetMomentum( -p*MeV );
119 } 119 }
120 */ 120 */
121 121
122 theResult->AddSecondary( theSec, secI 122 theResult->AddSecondary( theSec, secID );
123 } 123 }
124 } 124 }
125 else 125 else
126 { 126 {
127 127
128 //For the case data does not provide fin 128 //For the case data does not provide final states
129 129
130 //G4cout << "products != NULL; iZ = " << 130 //G4cout << "products != NULL; iZ = " << iZ << ", iA = " << iA << G4endl;
131 131
132 // TK comment 132 // TK comment
133 // aTarg->ReturnTargetParticle()->Get4Mo 133 // aTarg->ReturnTargetParticle()->Get4Momentum has trouble, thus we use following
134 G4Fragment nucleus( iA + ipA , iZ + ipZ 134 G4Fragment nucleus( iA + ipA , iZ + ipZ , aTrack.Get4Momentum() + G4LorentzVector( G4ThreeVector(0,0,0) , G4IonTable::GetIonTable()->GetIon( iZ + ipZ , iA )->GetPDGMass() ) );
135 G4PhotonEvaporation photonEvaporation; 135 G4PhotonEvaporation photonEvaporation;
136 photonEvaporation.SetICM( TRUE ); 136 photonEvaporation.SetICM( TRUE );
137 G4FragmentVector* products_from_PE = pho 137 G4FragmentVector* products_from_PE = photonEvaporation.BreakItUp(nucleus);
138 G4FragmentVector::iterator it; 138 G4FragmentVector::iterator it;
139 139
140 for (it = products_from_PE->begin(); it 140 for (it = products_from_PE->begin(); it != products_from_PE->end(); it++)
141 { 141 {
142 if ( (*it)->GetZ_asInt() == iZ + ipZ 142 if ( (*it)->GetZ_asInt() == iZ + ipZ && (*it)->GetA_asInt() == iA + ipA ) needResidual = false;
143 G4DynamicParticle* theSec = new G4Dyn 143 G4DynamicParticle* theSec = new G4DynamicParticle;
144 if ( (*it)->GetParticleDefinition() ! 144 if ( (*it)->GetParticleDefinition() != NULL ) {
145 //G4cout << (*it)->GetParticleDefi 145 //G4cout << (*it)->GetParticleDefinition()->GetParticleName() << G4endl;
146 theSec->SetDefinition( (*it)->GetP 146 theSec->SetDefinition( (*it)->GetParticleDefinition() );
147 theSec->Set4Momentum( (*it)->GetMo 147 theSec->Set4Momentum( (*it)->GetMomentum() );
148 } else { 148 } else {
149 //G4cout << (*it)->GetZ_asInt() << 149 //G4cout << (*it)->GetZ_asInt() << " " << (*it)->GetA_asInt() << G4endl;
150 theSec->SetDefinition( G4IonTable: 150 theSec->SetDefinition( G4IonTable::GetIonTable()->GetIon( (*it)->GetZ_asInt() , (*it)->GetA_asInt() ) );
151 theSec->Set4Momentum( (*it)->GetMo 151 theSec->Set4Momentum( (*it)->GetMomentum() );
152 } 152 }
153 theResult->AddSecondary( theSec, secI 153 theResult->AddSecondary( theSec, secID );
154 } 154 }
155 delete products_from_PE; 155 delete products_from_PE;
156 } 156 }
157 157
158 //if necessary, generate residual nucleus 158 //if necessary, generate residual nucleus
159 if ( needResidual ) { 159 if ( needResidual ) {
160 G4DynamicParticle* residual = new G4Dyna 160 G4DynamicParticle* residual = new G4DynamicParticle;
161 residual->SetDefinition( G4IonTable::Get 161 residual->SetDefinition( G4IonTable::GetIonTable()->GetIon( iZ + ipZ , iA + ipA ) );
162 residual->SetMomentum( -p*MeV ); 162 residual->SetMomentum( -p*MeV );
163 theResult->AddSecondary( residual, secID 163 theResult->AddSecondary( residual, secID );
164 } 164 }
165 165
166 delete products; 166 delete products;
167 167
168 theResult->SetStatusChange( stopAndKill ); 168 theResult->SetStatusChange( stopAndKill );
169 169
170 return theResult; 170 return theResult;
171 171
172 } 172 }
173 173