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 // $Id: G4QuasiElasticChannel.cc,v 1.4 2008/04/24 13:26:19 gunter Exp $ >> 28 // GEANT4 tag $Name: geant4-09-02-patch-01 $ 27 // 29 // 28 30 29 // Author : Gunter Folger March 2007 31 // Author : Gunter Folger March 2007 30 // Modified by Mikhail Kossov. Apr2009, E/M co << 31 // Class Description 32 // Class Description 32 // Final state production model for theoretica 33 // Final state production model for theoretical models of hadron inelastic 33 // quasi elastic scattering in geant4; 34 // quasi elastic scattering in geant4; 34 // Class Description - End 35 // Class Description - End 35 // 36 // 36 // Modified: 37 // Modified: 37 // 20110805 M. Kelsey -- Follow change to G4V << 38 // 20110808 M. Kelsey -- Move #includes from << 39 38 40 #include "G4QuasiElasticChannel.hh" 39 #include "G4QuasiElasticChannel.hh" 41 40 42 #include "G4Fancy3DNucleus.hh" 41 #include "G4Fancy3DNucleus.hh" 43 #include "G4DynamicParticle.hh" << 44 #include "G4HadTmpUtil.hh" /* lrint */ << 45 #include "G4KineticTrack.hh" << 46 #include "G4KineticTrackVector.hh" << 47 #include "G4LorentzVector.hh" << 48 #include "G4Neutron.hh" << 49 #include "G4Nucleon.hh" << 50 #include "G4Nucleus.hh" << 51 #include "G4ParticleDefinition.hh" << 52 #include "G4ParticleTable.hh" << 53 #include "G4IonTable.hh" << 54 #include "G4QuasiElRatios.hh" << 55 #include "globals.hh" << 56 #include <vector> << 57 #include "G4PhysicsModelCatalog.hh" << 58 42 59 //#define debug_scatter << 60 43 >> 44 #include "G4HadTmpUtil.hh" //lrint >> 45 >> 46 //#define debug_scatter 61 47 62 G4QuasiElasticChannel::G4QuasiElasticChannel() 48 G4QuasiElasticChannel::G4QuasiElasticChannel() 63 : G4HadronicInteraction("QuasiElastic"), << 49 { 64 theQuasiElastic(new G4QuasiElRatios), << 50 theQuasiElastic=G4QuasiFreeRatios::GetPointer(); 65 the3DNucleus(new G4Fancy3DNucleus), << 66 secID(-1) { << 67 secID = G4PhysicsModelCatalog::GetModelID( " << 68 } 51 } 69 52 70 G4QuasiElasticChannel::~G4QuasiElasticChannel( 53 G4QuasiElasticChannel::~G4QuasiElasticChannel() 71 { << 54 {} 72 delete the3DNucleus; << 73 delete theQuasiElastic; << 74 } << 75 55 76 G4double G4QuasiElasticChannel::GetFraction(G4 56 G4double G4QuasiElasticChannel::GetFraction(G4Nucleus &theNucleus, 77 const G4DynamicParticle & thePrimary) << 57 const G4DynamicParticle & thePrimary) 78 { 58 { 79 #ifdef debug_scatter << 59 std::pair<G4double,G4double> ratios; 80 G4cout << "G4QuasiElasticChannel:: P=" < << 60 ratios=theQuasiElastic->GetRatios(thePrimary.GetTotalMomentum(), 81 << ", pPDG=" << thePrimary.GetDef << 61 thePrimary.GetDefinition()->GetPDGEncoding(), 82 << ", Z = " << theNucleus.GetZ_a << 62 G4lrint(theNucleus.GetZ()), 83 << ", N = " << theNucleus.GetN_a << 63 G4lrint(theNucleus.GetN()-theNucleus.GetZ())); 84 << ", A = " << theNucleus.GetA_a << 64 #ifdef debug_getFraction 85 #endif << 65 G4cout << "G4QuasiElasticChannel::ratios " << ratios.first << " x " <<ratios.second 86 << 66 << " = " << ratios.first*ratios.second << G4endl; 87 std::pair<G4double,G4double> ratios; << 67 #endif 88 ratios=theQuasiElastic->GetRatios(thePrimary << 68 89 thePrimary << 69 return ratios.first*ratios.second; 90 theNucleus << 91 theNucleus << 92 #ifdef debug_scatter << 93 G4cout << "G4QuasiElasticChannel::ratios << 94 << " = " << ratios.first*ratios. << 95 #endif << 96 << 97 return ratios.first*ratios.second; << 98 } 70 } 99 71 100 G4KineticTrackVector * G4QuasiElasticChannel:: 72 G4KineticTrackVector * G4QuasiElasticChannel::Scatter(G4Nucleus &theNucleus, 101 << 73 const G4DynamicParticle & thePrimary) 102 { 74 { 103 G4int A=theNucleus.GetA_asInt(); << 75 104 G4int Z=theNucleus.GetZ_asInt(); << 76 105 // build Nucleus and choose random nucleon << 77 G4int A=G4lrint(theNucleus.GetN()); 106 the3DNucleus->Init(theNucleus.GetA_asInt(),t << 78 // G4int Z=G4lrint(theNucleus.GetZ()); 107 const std::vector<G4Nucleon>& nucleons=the3D << 79 // build Nucleus and choose random nucleon to scatter with 108 G4double targetNucleusMass=the3DNucleus->Get << 80 the3DNucleus.Init(theNucleus.GetN(),theNucleus.GetZ()); 109 G4LorentzVector targetNucleus4Mom(0.,0.,0.,t << 81 const std::vector<G4Nucleon *> nucleons=the3DNucleus.GetNucleons(); 110 G4int index; << 82 111 do { << 83 G4int index; 112 index=G4lrint((A-1)*G4UniformRand()); << 84 do { 113 } while (index < 0 || index >= static_cast<G << 85 index=G4lrint((A-1)*G4UniformRand()); 114 << 86 115 const G4ParticleDefinition * pDef= nucleons[ << 87 } while (index < 0 || index >= static_cast<G4int>(nucleons.size())); 116 << 88 117 G4int resA=A - 1; << 89 // G4double an=G4UniformRand()*A; 118 G4int resZ=Z - static_cast<int>(pDef->GetPDG << 90 G4ParticleDefinition * pDef= nucleons[index]->GetDefinition(); 119 const G4ParticleDefinition* resDef; << 91 120 G4double residualNucleusMass; << 92 #ifdef debug_scatter 121 if(resZ) << 93 G4cout << " neutron - proton? A, Z, an, pdg" <<" " 122 { << 94 << A <<" "<<G4lrint(theNucleus.GetZ()) 123 resDef=G4ParticleTable::GetParticleTable() << 95 << " "<<an <<" " << pDef->GetParticleName()<< G4endl; 124 residualNucleusMass=resDef->GetPDGMass(); << 96 #endif 125 } << 97 // G4LorentzVector pNucleon(G4ThreeVector(0,0,0),pDef->GetPDGMass()); 126 else { << 98 G4LorentzVector pNucleon=nucleons[index]->Get4Momentum(); 127 resDef=G4Neutron::Neutron(); << 99 128 residualNucleusMass=resA * G4Neutron::Neut << 100 std::pair<G4LorentzVector,G4LorentzVector> result; 129 } << 101 130 #ifdef debug_scatter << 102 result=theQuasiElastic->Scatter(pDef->GetPDGEncoding(),pNucleon, 131 G4cout<<"G4QElChan::Scatter: neutron - pr << 103 thePrimary.GetDefinition()->GetPDGEncoding(), 132 <<pDef->GetParticleName()<<G4endl; << 104 thePrimary.Get4Momentum()); 133 #endif << 105 134 << 106 if (result.first.e() == 0.) 135 G4LorentzVector pNucleon=nucleons[index].Get << 107 { 136 G4double residualNucleusEnergy=std::sqrt(sqr << 108 G4cout << "Warning - G4QuasiElasticChannel::Scatter no scattering" << G4endl; 137 pNu << 109 return 0; //no scatter 138 pNucleon.setE(targetNucleusMass-residualNucl << 110 } 139 G4LorentzVector residualNucleus4Mom=targetNu << 140 << 141 std::pair<G4LorentzVector,G4LorentzVector> r << 142 << 143 result=theQuasiElastic->Scatter(pDef->GetPDG << 144 thePrimary.G << 145 thePrimary.G << 146 G4LorentzVector scatteredHadron4Mom; << 147 if (result.first.e() > 0.) << 148 scatteredHadron4Mom=result.second; << 149 else { //scatter failed << 150 //G4cout << "Warning - G4QuasiElasticChann << 151 //return 0; //no scatter << 152 scatteredHadron4Mom=thePrimary.Get4Momentu << 153 residualNucleus4Mom=G4LorentzVector(0.,0., << 154 resDef=G4ParticleTable::GetParticleTable() << 155 } << 156 111 157 #ifdef debug_scatter 112 #ifdef debug_scatter 158 G4LorentzVector EpConservation=pNucleon+theP << 113 G4LorentzVector EpConservation=pNucleon+thePrimary.Get4Momentum() 159 - result.firs << 114 - result.first - result.second; 160 if ( (EpConservation.vect().mag2() > .01*M << 115 if ( (EpConservation.vect().mag2() > 0.01*MeV*MeV ) 161 || (std::abs(EpConservation.e()) > 0.1 * << 116 || (std::abs(EpConservation.e()) > 0.1 * MeV ) ) 162 { << 117 { 163 G4cout << "Warning - G4QuasiElasticChannel << 118 G4cout << "Warning - G4QuasiElasticChannel::Scatter E-p non conservation : " 164 << EpConservation << G4endl; << 119 << EpConservation << G4endl; 165 } << 120 } 166 #endif 121 #endif 167 122 168 G4KineticTrackVector * ktv = new G4KineticTr << 123 G4KineticTrackVector * ktv; 169 G4KineticTrack * sPrim=new G4KineticTrack(th << 124 ktv=new G4KineticTrackVector(); 170 0. << 125 G4KineticTrack * sPrim=new G4KineticTrack(thePrimary.GetDefinition(), 171 sPrim->SetCreatorModelID( secID ); << 126 0.,G4ThreeVector(0), result.second); 172 ktv->push_back(sPrim); << 127 ktv->push_back(sPrim); 173 if (result.first.e() > 0.) << 128 G4KineticTrack * sNuc=new G4KineticTrack(pDef, 174 { << 129 0.,G4ThreeVector(0), result.first); 175 G4KineticTrack * sNuc=new G4KineticTrack(p << 130 ktv->push_back(sNuc); 176 sNuc->SetCreatorModelID( secID ); << 131 177 ktv->push_back(sNuc); << 178 } << 179 if(resZ || resA==1) // For the only neutron << 180 { << 181 G4KineticTrack * rNuc=new G4KineticTrack(r << 182 0.,G4ThreeVector(0) << 183 rNuc->SetCreatorModelID( secID ); << 184 ktv->push_back(rNuc); << 185 } << 186 else // The residual nucleus consists of onl << 187 { << 188 residualNucleus4Mom/=resA; // Split 4- << 189 for(G4int in=0; in<resA; in++) // Loop ove << 190 { << 191 G4KineticTrack* rNuc=new G4KineticTrack( << 192 0.,G4ThreeVector(0) << 193 rNuc->SetCreatorModelID( secID ); << 194 ktv->push_back(rNuc); << 195 } << 196 } << 197 #ifdef debug_scatter 132 #ifdef debug_scatter 198 G4cout<<"G4QElC::Scat: Nucleon: "<<result.fi << 133 G4cout << " scattered Nucleon : " << result.first << " mass " <<result.first.mag() << G4endl; 199 G4cout<<"G4QElC::Scat: Project: "<<result.se << 134 G4cout << " scattered Project : " << result.second << " mass " <<result.second.mag() << G4endl; 200 #endif 135 #endif 201 return ktv; << 136 return ktv; 202 } 137 } 203 138