- Timestamp:
- May 28, 2009, 4:26:57 PM (15 years ago)
- Location:
- trunk/source/processes/hadronic/models/theo_high_energy/src
- Files:
-
- 2 edited
-
G4QuasiElasticChannel.cc (modified) (4 diffs)
-
G4TheoFSGenerator.cc (modified) (3 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/source/processes/hadronic/models/theo_high_energy/src/G4QuasiElasticChannel.cc
r1007 r1055 25 25 // 26 26 // 27 // $Id: G4QuasiElasticChannel.cc,v 1. 4 2008/04/24 13:26:19 gunterExp $28 // GEANT4 tag $Name: geant4-09-0 2$27 // $Id: G4QuasiElasticChannel.cc,v 1.7 2009/04/09 08:28:42 mkossov Exp $ 28 // GEANT4 tag $Name: geant4-09-03-beta-cand-01 $ 29 29 // 30 30 31 31 // Author : Gunter Folger March 2007 32 // Modified by Mikhail Kossov. Apr2009, E/M conservation: ResidualNucleus is added (ResNuc) 32 33 // Class Description 33 34 // Final state production model for theoretical models of hadron inelastic … … 42 43 43 44 44 #include "G4HadTmpUtil.hh" 45 #include "G4HadTmpUtil.hh" //lrint 45 46 46 47 //#define debug_scatter … … 48 49 G4QuasiElasticChannel::G4QuasiElasticChannel() 49 50 { 50 51 theQuasiElastic=G4QuasiFreeRatios::GetPointer(); 51 52 } 52 53 … … 55 56 56 57 G4double G4QuasiElasticChannel::GetFraction(G4Nucleus &theNucleus, 57 58 const G4DynamicParticle & thePrimary) 58 59 { 59 std::pair<G4double,G4double> ratios; 60 ratios=theQuasiElastic->GetRatios(thePrimary.GetTotalMomentum(), 61 thePrimary.GetDefinition()->GetPDGEncoding(), 62 G4lrint(theNucleus.GetZ()), 63 G4lrint(theNucleus.GetN()-theNucleus.GetZ())); 64 #ifdef debug_getFraction 65 G4cout << "G4QuasiElasticChannel::ratios " << ratios.first << " x " <<ratios.second 66 << " = " << ratios.first*ratios.second << G4endl; 60 std::pair<G4double,G4double> ratios; 61 #ifdef debug_scatter 62 G4cout << "G4QuasiElasticChannel:: P=" << thePrimary.GetTotalMomentum() 63 << ", pPDG=" << thePrimary.GetDefinition()->GetPDGEncoding() 64 << ", Z = " << G4lrint(theNucleus.GetZ()) 65 << ", N = " << G4lrint(theNucleus.GetN()-theNucleus.GetZ()) << G4endl; 67 66 #endif 68 69 return ratios.first*ratios.second; 67 ratios=theQuasiElastic->GetRatios(thePrimary.GetTotalMomentum(), 68 thePrimary.GetDefinition()->GetPDGEncoding(), 69 G4lrint(theNucleus.GetZ()), 70 G4lrint(theNucleus.GetN()-theNucleus.GetZ())); 71 #ifdef debug_scatter 72 G4cout << "G4QuasiElasticChannel::ratios " << ratios.first << " x " <<ratios.second 73 << " = " << ratios.first*ratios.second << G4endl; 74 #endif 75 76 return ratios.first*ratios.second; 77 //return 0.; // Switch off quasielastic (temporary) M.K. 78 //return 1.; // Only quasielastic (temporary) M.K. (DANGEROSE! Crashes at A=1!) 70 79 } 71 80 72 81 G4KineticTrackVector * G4QuasiElasticChannel::Scatter(G4Nucleus &theNucleus, 73 82 const G4DynamicParticle & thePrimary) 74 83 { 75 76 77 G4int A=G4lrint(theNucleus.GetN()); 78 // G4int Z=G4lrint(theNucleus.GetZ()); 79 // build Nucleus and choose random nucleon to scatter with 80 the3DNucleus.Init(theNucleus.GetN(),theNucleus.GetZ()); 81 const std::vector<G4Nucleon *> nucleons=the3DNucleus.GetNucleons(); 82 83 G4int index; 84 do { 85 index=G4lrint((A-1)*G4UniformRand()); 86 87 } while (index < 0 || index >= static_cast<G4int>(nucleons.size())); 88 89 // G4double an=G4UniformRand()*A; 90 G4ParticleDefinition * pDef= nucleons[index]->GetDefinition(); 91 84 G4int A=G4lrint(theNucleus.GetN()); 85 G4int Z=G4lrint(theNucleus.GetZ()); // M.K. ResNuc 86 // build Nucleus and choose random nucleon to scatter with 87 the3DNucleus.Init(theNucleus.GetN(),theNucleus.GetZ()); 92 88 #ifdef debug_scatter 93 G4cout << " neutron - proton? A, Z, an, pdg" <<" " 94 << A <<" "<<G4lrint(theNucleus.GetZ()) 95 << " "<<an <<" " << pDef->GetParticleName()<< G4endl; 89 G4cout<<"G4QElC::Scat: Before GetNucleons " << G4endl; 96 90 #endif 97 // G4LorentzVector pNucleon(G4ThreeVector(0,0,0),pDef->GetPDGMass()); 98 G4LorentzVector pNucleon=nucleons[index]->Get4Momentum(); 99 100 std::pair<G4LorentzVector,G4LorentzVector> result; 101 102 result=theQuasiElastic->Scatter(pDef->GetPDGEncoding(),pNucleon, 103 thePrimary.GetDefinition()->GetPDGEncoding(), 104 thePrimary.Get4Momentum()); 105 106 if (result.first.e() == 0.) 107 { 108 G4cout << "Warning - G4QuasiElasticChannel::Scatter no scattering" << G4endl; 109 return 0; //no scatter 110 } 91 const std::vector<G4Nucleon *> nucleons=the3DNucleus.GetNucleons(); 92 G4double targetNucleusMass=the3DNucleus.GetMass(); // M.K. ResNuc 93 #ifdef debug_scatter 94 G4cout<<"G4QElC::Scat: targetMass = " << targetNucleusMass << G4endl; 95 #endif 96 G4LorentzVector targetNucleus4Mom(0.,0.,0.,targetNucleusMass); // M.K. ResNuc 97 G4int index; 98 do 99 { 100 index=G4lrint((A-1)*G4UniformRand()); 101 } while (index < 0 || index >= static_cast<G4int>(nucleons.size())); 102 #ifdef debug_scatter 103 G4cout<<"G4QElC::Scat: index = " << index << G4endl; 104 #endif 105 G4ParticleDefinition * pDef= nucleons[index]->GetDefinition(); 106 G4int resA=A-1; // M.K. ResNuc 107 G4int resZ=Z-static_cast<int>(pDef->GetPDGCharge()); // M.K. ResNuc 108 G4ParticleDefinition* resDef=G4Neutron::Neutron(); // Resolve t-p=nn problem M.K. ResNuc 109 G4double residualNucleusMass=resDef->GetPDGMass(); // M.K. ResNuc 110 if(resZ) 111 { 112 resDef=G4ParticleTable::GetParticleTable()->FindIon(resZ,resA,0,resZ);// M.K. ResNuc 113 residualNucleusMass=resDef->GetPDGMass(); // M.K. ResNuc 114 } 115 else residualNucleusMass*=resA; // resA=resN M.K. ResNuc 116 #ifdef debug_scatter 117 G4cout<<"G4QElChan::Scatter: neutron - proton? A ="<<A<<", Z="<<Z<<", projName=" 118 <<pDef->GetParticleName()<<G4endl; 119 #endif 120 // G4LorentzVector pNucleon(G4ThreeVector(0,0,0),pDef->GetPDGMass()); 121 G4LorentzVector pNucleon=nucleons[index]->Get4Momentum(); 122 G4double residualNucleusEnergy=std::sqrt(residualNucleusMass*residualNucleusMass+ 123 pNucleon.vect().mag2()); // M.K. ResNuc 124 pNucleon.setE(targetNucleusMass-residualNucleusEnergy); // M.K. ResNuc 125 G4LorentzVector residualNucleus4Mom=targetNucleus4Mom-pNucleon; // M.K. ResNuc 126 127 std::pair<G4LorentzVector,G4LorentzVector> result; 128 129 result=theQuasiElastic->Scatter(pDef->GetPDGEncoding(),pNucleon, 130 thePrimary.GetDefinition()->GetPDGEncoding(), 131 thePrimary.Get4Momentum()); 132 G4LorentzVector scatteredHadron4Mom=result.second; // M.K. ResNuc 133 if (result.first.e() <= 0.) 134 { 135 //G4cout << "Warning - G4QuasiElasticChannel::Scatter no scattering" << G4endl; 136 //return 0; //no scatter 137 G4LorentzVector scatteredHadron4Mom=thePrimary.Get4Momentum(); // M.K. ResNuc 138 residualNucleus4Mom=G4LorentzVector(0.,0.,0.,targetNucleusMass); // M.K. ResNuc 139 resDef=G4ParticleTable::GetParticleTable()->FindIon(Z,A,0,Z); // M.K. ResNuc 140 } 111 141 112 142 #ifdef debug_scatter 113 114 115 if ( (EpConservation.vect().mag2() > 0.01*MeV*MeV )116 ||(std::abs(EpConservation.e()) > 0.1 * MeV ) )117 118 119 120 143 G4LorentzVector EpConservation=pNucleon+thePrimary.Get4Momentum() 144 - result.first - result.second; 145 if ( (EpConservation.vect().mag2() > .01*MeV*MeV ) 146 || (std::abs(EpConservation.e()) > 0.1 * MeV ) ) 147 { 148 G4cout << "Warning - G4QuasiElasticChannel::Scatter E-p non conservation : " 149 << EpConservation << G4endl; 150 } 121 151 #endif 122 152 123 G4KineticTrackVector * ktv; 124 ktv=new G4KineticTrackVector(); 125 G4KineticTrack * sPrim=new G4KineticTrack(thePrimary.GetDefinition(), 126 0.,G4ThreeVector(0), result.second); 127 ktv->push_back(sPrim); 128 G4KineticTrack * sNuc=new G4KineticTrack(pDef, 129 0.,G4ThreeVector(0), result.first); 130 ktv->push_back(sNuc); 131 153 G4KineticTrackVector * ktv; 154 ktv=new G4KineticTrackVector(); 155 G4KineticTrack * sPrim=new G4KineticTrack(thePrimary.GetDefinition(), 156 0.,G4ThreeVector(0), scatteredHadron4Mom); 157 ktv->push_back(sPrim); 158 if (result.first.e() > 0.) 159 { 160 G4KineticTrack * sNuc=new G4KineticTrack(pDef, 0.,G4ThreeVector(0), result.first); 161 ktv->push_back(sNuc); 162 } 163 if(resZ || resA==1) // For the only neutron or for tnuclei with Z>0 M.K. ResNuc 164 { 165 G4KineticTrack * rNuc=new G4KineticTrack(resDef, 166 0.,G4ThreeVector(0), residualNucleus4Mom); // M.K. ResNuc 167 ktv->push_back(rNuc); // M.K. ResNuc 168 } 169 else // The residual nucleus consists of only neutrons M.K. ResNuc 170 { 171 residualNucleus4Mom/=resA; // Split 4-mom of A*n system equally M.K. ResNuc 172 for(G4int in=0; in<resA; in++) // Loop over neutrons in A*n system. M.K. ResNuc 173 { 174 G4KineticTrack* rNuc=new G4KineticTrack(resDef, 175 0.,G4ThreeVector(0), residualNucleus4Mom); // M.K. ResNuc 176 ktv->push_back(rNuc); // M.K. ResNuc 177 } 178 } 132 179 #ifdef debug_scatter 133 G4cout << " scattered Nucleon : " << result.first << " mass "<<result.first.mag() << G4endl;134 G4cout << " scattered Project : " << result.second << " mass " <<result.second.mag()<< G4endl;180 G4cout<<"G4QElC::Scat: Nucleon: "<<result.first <<" mass "<<result.first.mag() << G4endl; 181 G4cout<<"G4QElC::Scat: Project: "<<result.second<<" mass "<<result.second.mag()<< G4endl; 135 182 #endif 136 return ktv; 183 return ktv; 137 184 } -
trunk/source/processes/hadronic/models/theo_high_energy/src/G4TheoFSGenerator.cc
r1007 r1055 25 25 // 26 26 // 27 // $Id: G4TheoFSGenerator.cc,v 1. 9 2007/11/13 16:01:36 gunterExp $28 // GEANT4 tag $Name: geant4-09-0 2$27 // $Id: G4TheoFSGenerator.cc,v 1.11 2009/04/09 08:28:42 mkossov Exp $ 28 // GEANT4 tag $Name: geant4-09-03-beta-cand-01 $ 29 29 // 30 30 // G4TheoFSGenerator … … 91 91 if ( theQuasielastic->GetFraction(theNucleus, *aPart) > G4UniformRand() ) 92 92 { 93 G4KineticTrackVector *result= theQuasielastic->Scatter(theNucleus, *aPart); 93 //G4cout<<"___G4TheoFSGenerator: before Scatter (1) QE=" << theQuasielastic<<G4endl; 94 G4KineticTrackVector *result= theQuasielastic->Scatter(theNucleus, *aPart); 95 //G4cout << "^^G4TheoFSGenerator: after Scatter (1) " << G4endl; 96 if (result) 97 { 98 for(unsigned int i=0; i<result->size(); i++) 99 { 100 G4DynamicParticle * aNew = 101 new G4DynamicParticle(result->operator[](i)->GetDefinition(), 102 result->operator[](i)->Get4Momentum().e(), 103 result->operator[](i)->Get4Momentum().vect()); 104 theParticleChange->AddSecondary(aNew); 105 delete result->operator[](i); 106 } 107 delete result; 108 109 } else 110 { 111 theParticleChange->SetStatusChange(isAlive); 112 theParticleChange->SetEnergyChange(thePrimary.GetKineticEnergy()); 113 theParticleChange->SetMomentumChange(thePrimary.Get4Momentum().vect().unit()); 114 115 } 116 return theParticleChange; 117 } 118 } 119 if ( theProjectileDiffraction) { 120 121 if ( theProjectileDiffraction->GetFraction(theNucleus, *aPart) > G4UniformRand() ) 122 // strictly, this returns fraction on inelastic, so quasielastic should 123 // already be substracted, ie. quasielastic should always be used 124 // before diffractive 125 { 126 //G4cout << "____G4TheoFSGenerator: before Scatter (2) " << G4endl; 127 G4KineticTrackVector *result= theProjectileDiffraction->Scatter(theNucleus, *aPart); 128 //G4cout << "^^^^G4TheoFSGenerator: after Scatter (2) " << G4endl; 94 129 if (result) 95 130 { … … 115 150 } 116 151 } 117 if ( theProjectileDiffraction) { 118 119 if ( theProjectileDiffraction->GetFraction(theNucleus, *aPart) > G4UniformRand() ) 120 // strictly, this returns fraction on inelastic, so quasielastic should 121 // already be substracted, ie. quasielastic should always be used 122 // before diffractive 123 { 124 G4KineticTrackVector *result= theProjectileDiffraction->Scatter(theNucleus, *aPart); 125 if (result) 126 { 127 for(unsigned int i=0; i<result->size(); i++) 128 { 129 G4DynamicParticle * aNew = 130 new G4DynamicParticle(result->operator[](i)->GetDefinition(), 131 result->operator[](i)->Get4Momentum().e(), 132 result->operator[](i)->Get4Momentum().vect()); 133 theParticleChange->AddSecondary(aNew); 134 delete result->operator[](i); 135 } 136 delete result; 137 138 } else 139 { 140 theParticleChange->SetStatusChange(isAlive); 141 theParticleChange->SetEnergyChange(thePrimary.GetKineticEnergy()); 142 theParticleChange->SetMomentumChange(thePrimary.Get4Momentum().vect().unit()); 143 144 } 145 return theParticleChange; 146 } 147 } 148 152 //G4cout << "____G4TheoFSGenerator: before Scatter (3) " << G4endl; 149 153 G4KineticTrackVector * theInitialResult = 150 154 theHighEnergyGenerator->Scatter(theNucleus, *aPart); 151 155 //G4cout << "^^^^G4TheoFSGenerator: after Scatter (3) " << G4endl; 152 156 153 157 G4ReactionProductVector * theTransportResult = NULL;
Note: See TracChangeset
for help on using the changeset viewer.
