Geant4 Cross Reference |
1 // 1 2 // ******************************************* 3 // * License and Disclaimer 4 // * 5 // * The Geant4 software is copyright of th 6 // * the Geant4 Collaboration. It is provided 7 // * conditions of the Geant4 Software License 8 // * LICENSE and available at http://cern.ch/ 9 // * include a list of copyright holders. 10 // * 11 // * Neither the authors of this software syst 12 // * institutes,nor the agencies providing fin 13 // * work make any representation or warran 14 // * regarding this software system or assum 15 // * use. Please see the license in the file 16 // * for the full disclaimer and the limitatio 17 // * 18 // * This code implementation is the result 19 // * technical work of the GEANT4 collaboratio 20 // * By using, copying, modifying or distri 21 // * any work based on the software) you ag 22 // * use in resulting scientific publicati 23 // * acceptance of all terms of the Geant4 Sof 24 // ******************************************* 25 // 26 // 27 //-------------------------------------------- 28 // 29 // ClassName: G4CRCoalescence ("CR" stand 30 // 31 // Author: 2020 Alberto Ribon , based on 32 // Diego Mauricio Gomez Coral fo 33 // 34 // Description: This class can be optionally 35 // 36 // G4TheoFSGenerator::ApplyYou 37 // 38 // to coalesce pairs of proton-n 39 // into deuterons and antideuter 40 // of secondaries produced by a 41 // This class can be useful in p 42 // applications. 43 // By default, this class is not 44 // However, it can be enabled vi 45 // 46 // /process/had/enableCRCoales 47 // 48 // It is assumed that the candid 49 // antiproton-antideuteron pairs 50 // spatial position, so the cond 51 // into account only their close 52 // 53 // This class is based entirely 54 // Diego Mauricio Gomez Coral fo 55 // The main application of this 56 // 57 // Notes: 58 // - In its current version, co 59 // proton projectile (because 60 // for deuteron and antideute 61 // only for the case of proto 62 // - This class is not meant be 63 // by intranuclear cascade mo 64 // INCL - which should have a 65 // 66 // Modified: 67 // 68 //-------------------------------------------- 69 // 70 #include "G4DynamicParticle.hh" 71 #include "G4Proton.hh" 72 #include "G4Neutron.hh" 73 #include "G4Deuteron.hh" 74 #include "G4AntiProton.hh" 75 #include "G4AntiNeutron.hh" 76 #include "G4CRCoalescence.hh" 77 #include "G4ReactionProduct.hh" 78 #include "G4IonTable.hh" 79 #include "G4PhysicsModelCatalog.hh" 80 81 82 G4CRCoalescence::G4CRCoalescence() : G4Hadroni 83 fP0_d( 0.0 ), fP0_dbar( 0.0 ), secID( -1 ) { 84 secID = G4PhysicsModelCatalog::GetModelID( " 85 } 86 87 88 G4CRCoalescence::~G4CRCoalescence() {} 89 90 91 void G4CRCoalescence::SetP0Coalescence( const 92 // Note by A.R. : in the present version, th 93 // proton projectile. If we w 94 // for other applications, be 95 // coalescence parameters als 96 // (Note that the method "Gen 97 // as it is for all projecti 98 fP0_dbar = 0.0; 99 fP0_d = 0.0; 100 if ( thePrimary.GetDefinition()->GetPDGEncod 101 G4double mproj = thePrimary.GetDefinition( 102 G4double pz = thePrimary.Get4Momentum().z( 103 G4double ekin = std::sqrt( pz*pz + mproj*m 104 if ( ekin > 10.0 ) { 105 fP0_dbar = 130.0 / ( 1.0 + std::exp( 21. 106 fP0_d = 118.1 * ( 1.0 + std::exp( 5.5 107 } 108 } 109 //G4cout << "Coalescence parameter p0 deuter 110 } 111 112 113 void G4CRCoalescence::GenerateDeuterons( G4Rea 114 // Deuteron clusters are made with the first 115 // the coalescence conditions, starting with 116 // A deuteron is a pair (i,j) where i is the 117 // with the relative momentum less than p0 ( 118 // The same applies for antideuteron cluster 119 // instead of protons and neutrons, respecti 120 121 // Vectors of index-position and 3-momentum 122 // protons, neutrons, antiprotons and antine 123 std::vector< std::pair< G4int, G4ThreeVector 124 std::vector< std::pair< G4int, G4ThreeVector 125 std::vector< std::pair< G4int, G4ThreeVector 126 std::vector< std::pair< G4int, G4ThreeVector 127 for ( unsigned int i = 0; i < result->size() 128 G4int pdgid = result->operator[](i)->GetDe 129 if ( pdgid == 2212 ) { // proton 130 proton.push_back( std::make_pair( i, res 131 result->erase( result->begin() + i ); 132 } 133 } 134 for ( unsigned int i = 0; i < result->size() 135 G4int pdgid = result->operator[](i)->GetDe 136 if ( pdgid == 2112 ) { // neutron 137 neutron.push_back( std::make_pair( i, re 138 result->erase( result->begin() + i ); 139 } 140 } 141 for ( unsigned int i = 0; i < result->size() 142 G4int pdgid = result->operator[](i)->GetDe 143 if ( pdgid == -2212 ) { // antiproton 144 antiproton.push_back( std::make_pair( i, 145 result->erase( result->begin() + i ); 146 } 147 } 148 for ( unsigned int i = 0; i < result->size() 149 G4int pdgid = result->operator[](i)->GetDe 150 if ( pdgid == -2112 ) { // antineutron 151 antineutron.push_back( std::make_pair( i 152 result->erase( result->begin() + i ); 153 } 154 } 155 156 for ( unsigned int i = 0; i < proton.size(); 157 if ( proton.at(i).first == -1 ) continue; 158 G4ThreeVector p1 = proton.at(i).second; 159 int partner1 = FindPartner( p1, G4Proton:: 160 G4Neutron::Neutron()->GetPDGMass(), 1 161 if ( partner1 == -1 ) { // if no partner 162 G4ParticleDefinition* prt = G4ParticleTa 163 G4ReactionProduct* finalp = new G4Reacti 164 finalp->SetDefinition( prt ); 165 G4double massp = prt->GetPDGMass(); 166 G4double totalEnergy = std::sqrt( p1.mag 167 finalp->SetMomentum( p1 ); 168 finalp->SetTotalEnergy( totalEnergy ); 169 finalp->SetMass( massp ); 170 result->push_back( finalp ); 171 continue; 172 } 173 G4ThreeVector p2 = neutron.at(partner1).se 174 PushDeuteron( p1, p2, 1, result ); 175 neutron.at(partner1).first = -1; // tag t 176 } 177 178 for ( unsigned int i = 0; i < neutron.size() 179 if ( neutron.at(i).first == -1 ) continue; 180 G4ParticleDefinition* nrt = G4ParticleTabl 181 G4ReactionProduct* finaln = new G4Reaction 182 finaln->SetDefinition( nrt ); 183 G4ThreeVector p2 = neutron.at(i).second; 184 G4double massn = nrt->GetPDGMass(); 185 G4double totalEnergy = std::sqrt( p2.mag() 186 finaln->SetMomentum( p2 ); 187 finaln->SetTotalEnergy( totalEnergy ); 188 finaln->SetMass( massn ); 189 result->push_back( finaln ); 190 } 191 192 for ( unsigned int i = 0; i < antiproton.siz 193 if ( antiproton.at(i).first == -1 ) contin 194 G4ThreeVector p1 = antiproton.at(i).second 195 int partner1 = FindPartner( p1, G4Proton:: 196 G4Neutron::Neutron()->GetPDGMass(), -1 197 if ( partner1 == -1 ) { // if no partner 198 G4ParticleDefinition* pbar = G4ParticleT 199 G4ReactionProduct* finalpbar = new G4Rea 200 finalpbar->SetDefinition( pbar ); 201 G4double massp = pbar->GetPDGMass(); 202 G4double totalEnergy = std::sqrt( p1.mag 203 finalpbar->SetMomentum( p1 ); 204 finalpbar->SetTotalEnergy( totalEnergy ) 205 finalpbar->SetMass( massp ); 206 result->push_back( finalpbar ); 207 continue; 208 } 209 G4ThreeVector p2 = antineutron.at(partner1 210 PushDeuteron( p1, p2, -1, result ); 211 antineutron.at(partner1).first = -1; // t 212 } 213 214 for ( unsigned int i = 0; i < antineutron.si 215 if ( antineutron.at(i).first == -1 ) conti 216 G4ParticleDefinition* nbar = G4ParticleTab 217 G4ReactionProduct* finalnbar = new G4React 218 finalnbar->SetDefinition( nbar ); 219 G4ThreeVector p2 = antineutron.at(i).secon 220 G4double massn = nbar->GetPDGMass(); 221 G4double totalEnergy = std::sqrt( p2.mag() 222 finalnbar->SetMomentum( p2 ); 223 finalnbar->SetTotalEnergy( totalEnergy ); 224 finalnbar->SetMass( massn ); 225 result->push_back( finalnbar ); 226 } 227 } 228 229 230 void G4CRCoalescence::PushDeuteron( const G4Th 231 G4ReactionProductVector* result ) 232 // Create a deuteron or antideuteron (depend 233 // from the two input momenta "p1" and "p2", 234 if ( charge > 0 ) { 235 G4ParticleDefinition* deuteron = G4Particl 236 G4ReactionProduct* finaldeut = new G4React 237 finaldeut->SetDefinition( deuteron ); 238 G4ThreeVector psum = p1 + p2; 239 G4double massd = deuteron->GetPDGMass(); 240 G4double totalEnergy = std::sqrt( psum.mag 241 finaldeut->SetMomentum( psum ); 242 finaldeut->SetTotalEnergy( totalEnergy ); 243 finaldeut->SetMass( massd ); 244 finaldeut->SetCreatorModelID( secID ); 245 result->push_back( finaldeut ); 246 } else { 247 G4ParticleDefinition* antideuteron = G4Par 248 G4ReactionProduct* finalantideut = new G4R 249 finalantideut->SetDefinition( antideuteron 250 G4ThreeVector psum = p1 + p2; 251 G4double massd = antideuteron->GetPDGMass( 252 G4double totalEnergy = std::sqrt( psum.mag 253 finalantideut->SetMomentum( psum ); 254 finalantideut->SetTotalEnergy( totalEnergy 255 finalantideut->SetMass( massd ); 256 finalantideut->SetCreatorModelID( secID ); 257 result->push_back( finalantideut ); 258 } 259 } 260 261 262 G4int G4CRCoalescence::FindPartner( const G4Th 263 std::vector< std::pair< G4int, G4T 264 G4double m2, G4int charge ) { 265 // Find a nucleon/antinucleon (depending on 266 // (which is a vector of either neutron or a 267 // within a sphere of radius p0 centered at 268 // particles (neutrons or antineutrons depen 269 for ( unsigned int j = 0; j < neutron.size() 270 if ( neutron.at(j).first == -1 ) continue; 271 G4ThreeVector p2 = neutron.at(j).second; 272 if ( ! Coalescence( p1, m1, p2, m2, charge 273 return j; 274 } 275 return -1; // no partner found 276 } 277 278 279 G4bool G4CRCoalescence::Coalescence( const G4T 280 const G4ThreeVector &p2, G4double 281 // Returns true if the momenta of the two nu 282 // inside of an sphere of radius p0 (assumin 283 return Coalescence( p1.x(), p1.y(), p1.z(), 284 } 285 286 287 G4bool G4CRCoalescence::Coalescence( G4double 288 G4double 289 G4int charge ) { 290 // Returns true if the momenta of the two nu 291 // inside of a sphere of radius p0 (assuming 292 G4double deltaP = GetPcm( p1x, p1y, p1z, m1, 293 if ( charge > 0 ) return ( deltaP < fP0_d ); 294 else return ( deltaP < fP0_dbar 295 } 296 297 298 G4double G4CRCoalescence::GetPcm( const G4Thre 299 const G4ThreeVector& p2, G4double m2 300 // Momentum in the center-of-mass frame of t 301 return GetPcm( p1.x(), p1.y(), p1.z(), m1, p 302 } 303 304 305 G4double G4CRCoalescence::GetPcm( G4double p1x 306 G4double p2x, G4double p2y, G4double 307 // Momentum in the center-of-mass frame of t 308 G4double scm = GetS( p1x, p1y, p1z, m1, p2x, 309 return std::sqrt( (scm - (m1-m2)*(m1-m2))*(s 310 } 311 312 313 G4double G4CRCoalescence::GetS( G4double p1x, 314 G4double p2x, G4double p2y, G4double p 315 // Square of center-of-mass energy of two pa 316 G4double E1 = std::sqrt( p1x*p1x + p1y*p1y + 317 G4double E2 = std::sqrt( p2x*p2x + p2y*p2y + 318 return (E1+E2)*(E1+E2) - (p1x+p2x)*(p1x+p2x) 319 } 320