ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4LENDCapture.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4LENDCapture.cc
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 
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 
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 
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 );
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 );
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 );
164  }
165 
166  delete products;
167 
168  theResult->SetStatusChange( stopAndKill );
169 
170  return theResult;
171 
172 }