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