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: G4ElasticHNScattering.cc 100828 2016-11-02 15:25:59Z gcosmo $ 27 // 28 // 28 29 29 // ------------------------------------------- 30 // ------------------------------------------------------------ 30 // GEANT 4 class implemetation file 31 // GEANT 4 class implemetation file 31 // 32 // 32 // ---------------- G4ElasticHNScattering 33 // ---------------- G4ElasticHNScattering -------------- 33 // by V. Uzhinsky, March 200 34 // by V. Uzhinsky, March 2008. 34 // elastic scattering used by Frit 35 // elastic scattering used by Fritiof model 35 // Take a projectile and a tar 36 // Take a projectile and a target 36 // scatter the projectile and 37 // scatter the projectile and target 37 // ------------------------------------------- 38 // --------------------------------------------------------------------- 38 39 39 #include "globals.hh" 40 #include "globals.hh" 40 #include "Randomize.hh" 41 #include "Randomize.hh" 41 #include "G4PhysicalConstants.hh" 42 #include "G4PhysicalConstants.hh" 42 #include "G4SystemOfUnits.hh" << 43 43 44 #include "G4ElasticHNScattering.hh" 44 #include "G4ElasticHNScattering.hh" 45 #include "G4LorentzRotation.hh" 45 #include "G4LorentzRotation.hh" 46 #include "G4ThreeVector.hh" 46 #include "G4ThreeVector.hh" 47 #include "G4ParticleDefinition.hh" 47 #include "G4ParticleDefinition.hh" 48 #include "G4VSplitableHadron.hh" 48 #include "G4VSplitableHadron.hh" 49 #include "G4ExcitedString.hh" 49 #include "G4ExcitedString.hh" 50 #include "G4FTFParameters.hh" 50 #include "G4FTFParameters.hh" 51 51 52 #include "G4SampleResonance.hh" 52 #include "G4SampleResonance.hh" 53 53 54 #include "G4Exp.hh" 54 #include "G4Exp.hh" 55 #include "G4Log.hh" 55 #include "G4Log.hh" 56 56 57 //============================================ 57 //============================================================================ 58 58 59 G4ElasticHNScattering::G4ElasticHNScattering() 59 G4ElasticHNScattering::G4ElasticHNScattering() {} 60 60 61 61 62 //============================================ 62 //============================================================================ 63 63 64 G4bool G4ElasticHNScattering::ElasticScatterin 64 G4bool G4ElasticHNScattering::ElasticScattering( G4VSplitableHadron* projectile, 65 65 G4VSplitableHadron* target, 66 66 G4FTFParameters* theParameters ) const { 67 projectile->IncrementCollisionCount( 1 ); 67 projectile->IncrementCollisionCount( 1 ); 68 target->IncrementCollisionCount( 1 ); 68 target->IncrementCollisionCount( 1 ); 69 69 70 if ( projectile->Get4Momentum().z() < 0.0 ) << 70 G4SampleResonance BrW; 71 71 72 // Projectile parameters 72 // Projectile parameters 73 G4LorentzVector Pprojectile = projectile->Ge 73 G4LorentzVector Pprojectile = projectile->Get4Momentum(); 74 G4double M0projectile = Pprojectile.mag(); << 74 if ( Pprojectile.z() < 0.0 ) return false; >> 75 G4bool PutOnMassShell( false ); >> 76 G4double M0projectile = Pprojectile.mag(); >> 77 //if ( M0projectile < projectile->GetDefinition()->GetPDGMass() ) { >> 78 >> 79 G4double MminProjectile=BrW.GetMinimumMass(projectile->GetDefinition()); >> 80 >> 81 if ( M0projectile < MminProjectile ) { >> 82 PutOnMassShell = true; >> 83 M0projectile = projectile->GetDefinition()->GetPDGMass(); >> 84 } 75 G4double M0projectile2 = M0projectile * M0pr 85 G4double M0projectile2 = M0projectile * M0projectile; >> 86 G4double AveragePt2 = theParameters->GetAvaragePt2ofElasticScattering(); 76 87 77 // Target parameters 88 // Target parameters 78 G4LorentzVector Ptarget = target->Get4Moment 89 G4LorentzVector Ptarget = target->Get4Momentum(); 79 G4double M0target = Ptarget.mag(); 90 G4double M0target = Ptarget.mag(); 80 G4double M0target2 = M0target * M0target; << 91 //if ( M0target < target->GetDefinition()->GetPDGMass() ) { 81 92 82 G4double AveragePt2 = theParameters->GetAvar << 93 G4double MminTarget=BrW.GetMinimumMass(target->GetDefinition()); >> 94 >> 95 if ( M0target < MminTarget ) { >> 96 PutOnMassShell = true; >> 97 M0target = target->GetDefinition()->GetPDGMass(); >> 98 } >> 99 G4double M0target2 = M0target * M0target; 83 100 84 // Transform momenta to cms and then rotate 101 // Transform momenta to cms and then rotate parallel to z axis; 85 G4LorentzVector Psum; 102 G4LorentzVector Psum; 86 Psum = Pprojectile + Ptarget; 103 Psum = Pprojectile + Ptarget; 87 G4LorentzRotation toCms( -1*Psum.boostVector 104 G4LorentzRotation toCms( -1*Psum.boostVector() ); 88 G4LorentzVector Ptmp = toCms*Pprojectile; 105 G4LorentzVector Ptmp = toCms*Pprojectile; 89 if ( Ptmp.pz() <= 0.0 ) return false; 106 if ( Ptmp.pz() <= 0.0 ) return false; 90 //"String" moving backwards in CMS, abort c 107 //"String" moving backwards in CMS, abort collision ! 91 //G4cout << " abort Collision! " << G4endl; 108 //G4cout << " abort Collision! " << G4endl; 92 toCms.rotateZ( -1*Ptmp.phi() ); 109 toCms.rotateZ( -1*Ptmp.phi() ); 93 toCms.rotateY( -1*Ptmp.theta() ); 110 toCms.rotateY( -1*Ptmp.theta() ); 94 G4LorentzRotation toLab( toCms.inverse() ); 111 G4LorentzRotation toLab( toCms.inverse() ); 95 Pprojectile.transform( toCms ); 112 Pprojectile.transform( toCms ); 96 Ptarget.transform( toCms ); 113 Ptarget.transform( toCms ); 97 114 >> 115 // Putting on mass-on-shell, if needed 98 G4double PZcms2, PZcms; 116 G4double PZcms2, PZcms; 99 G4double S = Psum.mag2(); 117 G4double S = Psum.mag2(); 100 G4double SqrtS = std::sqrt( S ); 118 G4double SqrtS = std::sqrt( S ); 101 if ( SqrtS < M0projectile + M0target ) retur 119 if ( SqrtS < M0projectile + M0target ) return false; 102 120 103 PZcms2 = ( S*S + sqr( M0projectile2 ) + sqr( 121 PZcms2 = ( S*S + sqr( M0projectile2 ) + sqr( M0target2 ) 104 - 2*S*M0projectile2 - 2*S*M0targe 122 - 2*S*M0projectile2 - 2*S*M0target2 - 2*M0projectile2*M0target2 ) / 4.0 / S; 105 123 106 PZcms = ( PZcms2 > 0.0 ? std::sqrt( PZcms2 ) << 124 if ( PZcms2 < 0.0 ) { // It can be in an interaction with off-shell nuclear nucleon >> 125 if ( M0projectile > projectile->GetDefinition()->GetPDGMass() ) { >> 126 // An attempt to de-excite the projectile >> 127 // It is assumed that the target is in the ground state >> 128 M0projectile = projectile->GetDefinition()->GetPDGMass(); >> 129 M0projectile2 = M0projectile * M0projectile; >> 130 PZcms2= ( S*S + sqr( M0projectile2 ) + sqr( M0target2 ) >> 131 - 2*S*M0projectile2 - 2*S*M0target2 - 2*M0projectile2*M0target2 ) / 4.0 / S; >> 132 if ( PZcms2 < 0.0 ) { return false; } // Nonsuccesful attempt to de-excitate the projectile >> 133 } else { >> 134 return false; // The projectile was not excited, but the energy was too low to put >> 135 // the target nucleon on mass-shell >> 136 } >> 137 } >> 138 >> 139 PZcms = std::sqrt( PZcms2 ); >> 140 >> 141 if ( PutOnMassShell ) { >> 142 if ( Pprojectile.z() > 0.0 ) { >> 143 Pprojectile.setPz( PZcms ); >> 144 Ptarget.setPz( -PZcms ); >> 145 } else { >> 146 Pprojectile.setPz( -PZcms ); >> 147 Ptarget.setPz( PZcms ); >> 148 }; >> 149 Pprojectile.setE( std::sqrt( M0projectile2 + Pprojectile.x() * Pprojectile.x() + >> 150 Pprojectile.y() * Pprojectile.y() + PZcms2 ) ); >> 151 Ptarget.setE( std::sqrt( M0target2 + Ptarget.x() * Ptarget.x() + Ptarget.y() * Ptarget.y() + >> 152 PZcms2 ) ); >> 153 } 107 154 108 G4double maxPtSquare = PZcms2; 155 G4double maxPtSquare = PZcms2; 109 156 110 // Now we can calculate the transferred Pt 157 // Now we can calculate the transferred Pt 111 G4double Pt2; 158 G4double Pt2; 112 G4double ProjMassT2, ProjMassT; 159 G4double ProjMassT2, ProjMassT; 113 G4double TargMassT2, TargMassT; 160 G4double TargMassT2, TargMassT; 114 G4LorentzVector Qmomentum; 161 G4LorentzVector Qmomentum; 115 162 116 const G4int maxNumberOfLoops = 1000; 163 const G4int maxNumberOfLoops = 1000; 117 G4int loopCounter = 0; 164 G4int loopCounter = 0; 118 do { 165 do { 119 Qmomentum = G4LorentzVector( GaussianPt( A 166 Qmomentum = G4LorentzVector( GaussianPt( AveragePt2, maxPtSquare ), 0.0 ); 120 Pt2 = G4ThreeVector( Qmomentum.vect() ).ma 167 Pt2 = G4ThreeVector( Qmomentum.vect() ).mag2(); 121 ProjMassT2 = M0projectile2 + Pt2; 168 ProjMassT2 = M0projectile2 + Pt2; 122 ProjMassT = std::sqrt( ProjMassT2 ); 169 ProjMassT = std::sqrt( ProjMassT2 ); 123 TargMassT2 = M0target2 + Pt2; 170 TargMassT2 = M0target2 + Pt2; 124 TargMassT = std::sqrt( TargMassT2 ); 171 TargMassT = std::sqrt( TargMassT2 ); 125 } while ( ( SqrtS < ProjMassT + TargMassT ) 172 } while ( ( SqrtS < ProjMassT + TargMassT ) && 126 ++loopCounter < maxNumberOfLoops ) 173 ++loopCounter < maxNumberOfLoops ); /* Loop checking, 10.08.2015, A.Ribon */ 127 if ( loopCounter >= maxNumberOfLoops ) { 174 if ( loopCounter >= maxNumberOfLoops ) { 128 return false; 175 return false; 129 } 176 } 130 177 131 PZcms2 = ( S*S + sqr( ProjMassT2 ) + sqr( Ta 178 PZcms2 = ( S*S + sqr( ProjMassT2 ) + sqr( TargMassT2 ) 132 - 2.0*S*ProjMassT2 - 2.0*S*TargMa 179 - 2.0*S*ProjMassT2 - 2.0*S*TargMassT2 - 2.0*ProjMassT2*TargMassT2 ) / 4.0 / S; 133 180 134 if ( PZcms2 < 0.0 ) { PZcms2 = 0.0; }; // t 181 if ( PZcms2 < 0.0 ) { PZcms2 = 0.0; }; // to avoid the exactness problem 135 PZcms = std::sqrt( PZcms2 ); 182 PZcms = std::sqrt( PZcms2 ); 136 Pprojectile.setPz( PZcms ); 183 Pprojectile.setPz( PZcms ); 137 Ptarget.setPz( -PZcms ); 184 Ptarget.setPz( -PZcms ); 138 Pprojectile += Qmomentum; 185 Pprojectile += Qmomentum; 139 Ptarget -= Qmomentum; 186 Ptarget -= Qmomentum; 140 187 141 // Transform back and update SplitableHadron 188 // Transform back and update SplitableHadron Participant. 142 Pprojectile.transform( toLab ); 189 Pprojectile.transform( toLab ); 143 Ptarget.transform( toLab ); 190 Ptarget.transform( toLab ); 144 191 145 // Calculation of the creation time 192 // Calculation of the creation time 146 projectile->SetTimeOfCreation( target->GetTi 193 projectile->SetTimeOfCreation( target->GetTimeOfCreation() ); 147 projectile->SetPosition( target->GetPosition 194 projectile->SetPosition( target->GetPosition() ); 148 195 149 // Creation time and position of target nucl 196 // Creation time and position of target nucleon were determined at 150 // ReggeonCascade() of G4FTFModel 197 // ReggeonCascade() of G4FTFModel 151 198 152 projectile->Set4Momentum( Pprojectile ); 199 projectile->Set4Momentum( Pprojectile ); 153 target->Set4Momentum( Ptarget ); 200 target->Set4Momentum( Ptarget ); 154 201 155 //projectile->IncrementCollisionCount( 1 ); 202 //projectile->IncrementCollisionCount( 1 ); 156 //target->IncrementCollisionCount( 1 ); 203 //target->IncrementCollisionCount( 1 ); 157 204 158 return true; 205 return true; 159 } 206 } 160 207 161 208 162 //============================================ 209 //============================================================================ 163 210 164 G4ThreeVector G4ElasticHNScattering::GaussianP 211 G4ThreeVector G4ElasticHNScattering::GaussianPt( G4double AveragePt2, 165 212 G4double maxPtSquare ) const { 166 // @@ this method is used in FTFModel as wel 213 // @@ this method is used in FTFModel as well. Should go somewhere common! 167 G4double Pt2( 0.0 ); 214 G4double Pt2( 0.0 ); 168 if ( AveragePt2 <= 0.0 ) { 215 if ( AveragePt2 <= 0.0 ) { 169 Pt2 = 0.0; 216 Pt2 = 0.0; 170 } else { 217 } else { 171 Pt2 = -AveragePt2 * G4Log( 1.0 + G4Uniform << 218 Pt2 = -AveragePt2 * G4Log( 1.0 + G4UniformRand() * >> 219 ( G4Exp( -maxPtSquare/AveragePt2 ) -1.0 ) ); 172 } 220 } 173 G4double Pt = ( Pt2 > 0.0 ? std::sqrt( Pt2 ) << 221 G4double Pt = std::sqrt( Pt2 ); 174 G4double phi = G4UniformRand() * twopi; 222 G4double phi = G4UniformRand() * twopi; 175 return G4ThreeVector( Pt * std::cos( phi ), 223 return G4ThreeVector( Pt * std::cos( phi ), Pt * std::sin( phi ), 0.0 ); 176 } 224 } 177 225 178 226 179 //============================================ 227 //============================================================================ 180 228 181 G4ElasticHNScattering::G4ElasticHNScattering( 229 G4ElasticHNScattering::G4ElasticHNScattering( const G4ElasticHNScattering& ) { 182 throw G4HadronicException( __FILE__, __LINE_ 230 throw G4HadronicException( __FILE__, __LINE__, 183 "G4ElasticHNScatt << 231 "G4ElasticHNScattering copy contructor not meant to be called" ); 184 } 232 } 185 233 186 234 187 //============================================ 235 //============================================================================ 188 236 189 G4ElasticHNScattering::~G4ElasticHNScattering( 237 G4ElasticHNScattering::~G4ElasticHNScattering() {} 190 238 191 239 192 //============================================ 240 //============================================================================ 193 241 194 const G4ElasticHNScattering & G4ElasticHNScatt 242 const G4ElasticHNScattering & G4ElasticHNScattering::operator=( const G4ElasticHNScattering& ) { 195 throw G4HadronicException( __FILE__, __LINE_ 243 throw G4HadronicException( __FILE__, __LINE__, 196 "G4ElasticHNScatt 244 "G4ElasticHNScattering = operator not meant to be called" ); 197 } 245 } 198 246 199 247 200 //============================================ 248 //============================================================================ 201 249 202 G4bool G4ElasticHNScattering::operator==( cons << 250 int G4ElasticHNScattering::operator==( const G4ElasticHNScattering& ) const { 203 throw G4HadronicException( __FILE__, __LINE__ 251 throw G4HadronicException( __FILE__, __LINE__, 204 "G4ElasticHNScatte 252 "G4ElasticHNScattering == operator not meant to be called" ); 205 } 253 } 206 254 207 255 208 //============================================ 256 //============================================================================ 209 257 210 G4bool G4ElasticHNScattering::operator!=( cons << 258 int G4ElasticHNScattering::operator!=( const G4ElasticHNScattering& ) const { 211 throw G4HadronicException( __FILE__, __LINE_ 259 throw G4HadronicException( __FILE__, __LINE__, 212 "G4ElasticHNScatte 260 "G4ElasticHNScattering != operator not meant to be called" ); 213 } 261 } 214 262 215 263