Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/hadronic/models/theo_high_energy/src/G4QuasiElasticChannel.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /processes/hadronic/models/theo_high_energy/src/G4QuasiElasticChannel.cc (Version 11.3.0) and /processes/hadronic/models/theo_high_energy/src/G4QuasiElasticChannel.cc (Version 9.6.p4)


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