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: G4DiffractiveSplitableHadron.cc 107317 2017-11-08 16:25:57Z gcosmo $ >> 28 // GEANT4 tag $Name: $ 27 // 29 // 28 30 29 // ------------------------------------------- 31 // ------------------------------------------------------------ 30 // GEANT 4 class implementation file 32 // GEANT 4 class implementation file 31 // 33 // 32 // ---------------- G4DiffractiveSplitabl 34 // ---------------- G4DiffractiveSplitableHadron---------------- 33 // by Gunter Folger, August 1998. 35 // by Gunter Folger, August 1998. 34 // class splitting an interacting partic 36 // class splitting an interacting particle. Used by FTF String Model. 35 // ------------------------------------------- 37 // ------------------------------------------------------------ 36 38 37 #include "G4DiffractiveSplitableHadron.hh" 39 #include "G4DiffractiveSplitableHadron.hh" 38 40 39 #include "G4ParticleDefinition.hh" 41 #include "G4ParticleDefinition.hh" 40 #include "Randomize.hh" 42 #include "Randomize.hh" 41 43 42 44 43 //============================================ 45 //============================================================================ 44 46 45 G4DiffractiveSplitableHadron::G4DiffractiveSpl 47 G4DiffractiveSplitableHadron::G4DiffractiveSplitableHadron() 46 { 48 { 47 PartonIndex = -1; 49 PartonIndex = -1; 48 G4LorentzVector tmp=G4LorentzVector(0.,0.,0. << 49 Parton[0] = new G4Parton( 1 ); 50 Parton[0] = new G4Parton( 1 ); 50 Parton[1] = new G4Parton(-1 ); 51 Parton[1] = new G4Parton(-1 ); 51 << 52 Parton[0]->Set4Momentum(tmp); Parton[1]->Set << 53 } 52 } 54 53 55 54 56 //============================================ 55 //============================================================================ 57 56 58 G4DiffractiveSplitableHadron::G4DiffractiveSpl 57 G4DiffractiveSplitableHadron::G4DiffractiveSplitableHadron( const G4ReactionProduct& aPrimary ) : 59 G4VSplitableHadron( aPrimary ) 58 G4VSplitableHadron( aPrimary ) 60 { 59 { 61 PartonIndex = -1; << 60 PartonIndex = -2; 62 Parton[0] = nullptr; << 61 Parton[0] = NULL; 63 Parton[1] = nullptr; << 64 } 62 } 65 63 66 64 67 //============================================ 65 //============================================================================ 68 66 69 G4DiffractiveSplitableHadron::G4DiffractiveSpl 67 G4DiffractiveSplitableHadron::G4DiffractiveSplitableHadron( const G4Nucleon& aNucleon ) : 70 G4VSplitableHadron( aNucleon ) 68 G4VSplitableHadron( aNucleon ) 71 { 69 { 72 PartonIndex = -1; << 70 PartonIndex = -2; 73 Parton[0] = nullptr; << 71 Parton[0] = NULL; 74 Parton[1] = nullptr; << 75 } 72 } 76 73 77 74 78 //============================================ 75 //============================================================================ 79 76 80 G4DiffractiveSplitableHadron::G4DiffractiveSpl 77 G4DiffractiveSplitableHadron::G4DiffractiveSplitableHadron( const G4VKineticNucleon* aNucleon ) : 81 G4VSplitableHadron( aNucleon ) 78 G4VSplitableHadron( aNucleon ) 82 { 79 { 83 PartonIndex = -1; << 80 PartonIndex = -2; 84 Parton[0] = nullptr; << 81 Parton[0] = NULL; 85 Parton[1] = nullptr; << 86 } 82 } 87 83 88 84 89 //============================================ 85 //============================================================================ 90 86 91 G4DiffractiveSplitableHadron::~G4DiffractiveSp << 87 G4DiffractiveSplitableHadron::~G4DiffractiveSplitableHadron() { >> 88 //G4cout << "Destruct G4DiffractiveSplitableHadron" << Parton[0] << " " << Parton[1] << G4endl; >> 89 //if ( Parton[0] != NULL ) { delete Parton[0]; delete Parton[1]; } >> 90 } 92 91 93 92 94 //============================================ 93 //============================================================================ 95 94 96 void G4DiffractiveSplitableHadron::SplitUp() { 95 void G4DiffractiveSplitableHadron::SplitUp() { 97 << 96 //G4cout << "SplitUp() IsSplit() Parton[0] " << IsSplit() << " " << Parton[0] << G4endl; 98 if ( IsSplit() ) return; 97 if ( IsSplit() ) return; 99 Splitting(); 98 Splitting(); 100 // Split once only... 99 // Split once only... 101 if ( Parton[0] != nullptr ) return; << 100 if ( Parton[0] != NULL ) return; 102 101 103 // flavours of quark ends 102 // flavours of quark ends 104 G4int PDGcode = GetDefinition()->GetPDGEncod 103 G4int PDGcode = GetDefinition()->GetPDGEncoding(); 105 G4int stringStart, stringEnd; 104 G4int stringStart, stringEnd; 106 ChooseStringEnds( PDGcode, &stringStart, &st 105 ChooseStringEnds( PDGcode, &stringStart, &stringEnd ); 107 106 108 Parton[0] = new G4Parton( stringStart ); 107 Parton[0] = new G4Parton( stringStart ); 109 Parton[1] = new G4Parton( stringEnd ); 108 Parton[1] = new G4Parton( stringEnd ); 110 109 111 G4LorentzVector tmp=G4LorentzVector(0.,0.,0. << 112 Parton[0]->Set4Momentum(tmp); Parton[1]->Set << 113 << 114 /* // 110 /* // Inversion of a string 115 if ( G4UniformRand() < 1.75 ) { //0.75 111 if ( G4UniformRand() < 1.75 ) { //0.75 116 Parton[0] = new G4Parton( stringStart ); 112 Parton[0] = new G4Parton( stringStart ); 117 Parton[1] = new G4Parton( stringEnd ); 113 Parton[1] = new G4Parton( stringEnd ); 118 } else { 114 } else { 119 Parton[0] = new G4Parton( stringEnd ); 115 Parton[0] = new G4Parton( stringEnd ); 120 Parton[1] = new G4Parton( stringStart ); 116 Parton[1] = new G4Parton( stringStart ); 121 } 117 } 122 */ 118 */ 123 119 124 PartonIndex = -1; 120 PartonIndex = -1; 125 } 121 } 126 122 127 123 128 //============================================ 124 //============================================================================ 129 125 130 G4Parton* G4DiffractiveSplitableHadron::GetNex 126 G4Parton* G4DiffractiveSplitableHadron::GetNextParton() { 131 ++PartonIndex; 127 ++PartonIndex; 132 if ( PartonIndex > 1 || PartonIndex < 0 ) << 128 if ( PartonIndex > 1 || PartonIndex < 0 ) return NULL; 133 G4int PartonInd( PartonIndex ); 129 G4int PartonInd( PartonIndex ); 134 if ( PartonIndex == 1 ) PartonIndex = -1; 130 if ( PartonIndex == 1 ) PartonIndex = -1; 135 return Parton[ PartonInd ]; 131 return Parton[ PartonInd ]; 136 } 132 } 137 133 138 134 139 //============================================ 135 //============================================================================ 140 136 141 G4Parton* G4DiffractiveSplitableHadron::GetNex 137 G4Parton* G4DiffractiveSplitableHadron::GetNextAntiParton() { 142 ++PartonIndex; 138 ++PartonIndex; 143 if ( PartonIndex > 1 || PartonIndex < 0 ) << 139 if ( PartonIndex > 1 || PartonIndex < 0 ) return NULL; 144 G4int PartonInd( PartonIndex ); 140 G4int PartonInd( PartonIndex ); 145 if ( PartonIndex == 1 ) PartonIndex = -1; 141 if ( PartonIndex == 1 ) PartonIndex = -1; 146 return Parton[ PartonInd ]; 142 return Parton[ PartonInd ]; 147 } 143 } 148 144 149 145 150 //============================================ 146 //============================================================================ 151 147 152 void G4DiffractiveSplitableHadron::SetFirstPar 148 void G4DiffractiveSplitableHadron::SetFirstParton( G4int PDGcode ) { 153 delete Parton[0]; 149 delete Parton[0]; 154 Parton[0] = new G4Parton( PDGcode ); 150 Parton[0] = new G4Parton( PDGcode ); 155 G4LorentzVector tmp=G4LorentzVector(0.,0.,0. << 156 Parton[0]->Set4Momentum(tmp); << 157 } 151 } 158 152 159 153 160 //============================================ 154 //============================================================================ 161 155 162 void G4DiffractiveSplitableHadron::SetSecondPa 156 void G4DiffractiveSplitableHadron::SetSecondParton( G4int PDGcode ) { 163 delete Parton[1]; 157 delete Parton[1]; 164 Parton[1] = new G4Parton( PDGcode ); 158 Parton[1] = new G4Parton( PDGcode ); 165 G4LorentzVector tmp=G4LorentzVector(0.,0.,0. << 166 Parton[1]->Set4Momentum(tmp); << 167 } 159 } 168 160 169 161 170 //============================================ 162 //============================================================================ 171 163 172 void G4DiffractiveSplitableHadron::ChooseStrin 164 void G4DiffractiveSplitableHadron::ChooseStringEnds( G4int PDGcode, G4int* aEnd, 173 165 G4int* bEnd ) const { 174 G4int absPDGcode = std::abs( PDGcode ); 166 G4int absPDGcode = std::abs( PDGcode ); 175 167 176 if ( absPDGcode < 1000 ) { //-------------- 168 if ( absPDGcode < 1000 ) { //-------------------- Meson ------------- 177 G4int heavy(0), light(0); 169 G4int heavy(0), light(0); 178 if (!((absPDGcode == 111)||(absPDGcode == << 170 if(!((absPDGcode == 111)||(absPDGcode == 221)||(absPDGcode == 331))) 179 { // Ordinary mes 171 { // Ordinary mesons ======================= 180 heavy = absPDGcode/100; 172 heavy = absPDGcode/100; 181 light = (absPDGcode % 100)/10; 173 light = (absPDGcode % 100)/10; 182 //G4int anti = std::pow( -1 , std::max( h 174 //G4int anti = std::pow( -1 , std::max( heavy, light ) ); 183 G4int anti = 1 - 2*( std::max( heavy, lig 175 G4int anti = 1 - 2*( std::max( heavy, light ) % 2 ); 184 if (PDGcode < 0 ) anti *= -1; 176 if (PDGcode < 0 ) anti *= -1; 185 heavy *= anti; 177 heavy *= anti; 186 light *= -1 * anti; 178 light *= -1 * anti; 187 } 179 } 188 else 180 else 189 { // Pi0, Eta, Eta 181 { // Pi0, Eta, Eta' ======================= 190 if ( G4UniformRand() < 0.5 ) {heavy = 1; << 182 if( G4UniformRand() < 0.5 ) {heavy = 1; light = -1;} 191 else {heavy = 2; << 183 else {heavy = 2; light = -2;} 192 } 184 } 193 if ( G4UniformRand() < 0.5 ) { 185 if ( G4UniformRand() < 0.5 ) { 194 *aEnd = heavy; 186 *aEnd = heavy; 195 *bEnd = light; 187 *bEnd = light; 196 } else { 188 } else { 197 *aEnd = light; 189 *aEnd = light; 198 *bEnd = heavy; 190 *bEnd = heavy; 199 } 191 } 200 } else { //-------------- 192 } else { //-------------------- Baryon -------------- 201 G4int j1000 = PDGcode/1000; 193 G4int j1000 = PDGcode/1000; 202 G4int j100 = (PDGcode % 1000)/100; 194 G4int j100 = (PDGcode % 1000)/100; 203 G4int j10 = (PDGcode % 100)/10; 195 G4int j10 = (PDGcode % 100)/10; 204 196 205 if ( absPDGcode > 4000 ) { << 206 *aEnd = j10; << 207 if ( G4UniformRand() > 0.25 ) { << 208 *bEnd = Diquark( j1000, j100, 0 ); << 209 } else { << 210 *bEnd = Diquark( j1000, j100, 1 ); << 211 } << 212 return; << 213 } << 214 << 215 G4double SuppresUUDDSS=1.0/2.0; 197 G4double SuppresUUDDSS=1.0/2.0; 216 if ((j1000 == j100) && (j1000 == j10)) Sup << 198 if((j1000 == j100) && (j1000 == j10)) SuppresUUDDSS=1.; 217 199 218 const G4int maxNumberOfLoops = 1000; 200 const G4int maxNumberOfLoops = 1000; 219 G4int loopCounter = 0; 201 G4int loopCounter = 0; 220 do 202 do 221 { 203 { 222 G4double random = G4UniformRand(); 204 G4double random = G4UniformRand(); 223 205 224 if (random < 0.33333) << 206 if(random < 0.33333) 225 { 207 { 226 if (( j100 == j10 ) && ( G4UniformRand << 208 if(( j100 == j10 ) && ( G4UniformRand() > SuppresUUDDSS )) continue; 227 *aEnd = j1000; 209 *aEnd = j1000; 228 if ( j100 == j10 ) {*bEnd << 210 if( j100 == j10 ) {*bEnd = Diquark( j100, j10, 1 );} 229 else 211 else 230 if ( G4UniformRand() > 0.25) {*bEnd << 212 if( G4UniformRand() > 0.25) {*bEnd = Diquark( j100, j10, 0 );} 231 else {*bEnd = 213 else {*bEnd = Diquark( j100, j10, 1 );} 232 break; 214 break; 233 } 215 } 234 else if (random < 0.66667) << 216 else if(random < 0.66667) 235 { 217 { 236 if (( j1000 == j10 ) && ( G4UniformRan << 218 if(( j1000 == j10 ) && ( G4UniformRand() > SuppresUUDDSS )) continue; 237 *aEnd = j100; 219 *aEnd = j100; 238 if ( j1000 == j10 ) {*bEnd << 220 if( j1000 == j10 ) {*bEnd = Diquark( j1000, j10, 1 );} 239 else 221 else 240 if ( G4UniformRand() > 0.25) {*bEnd << 222 if( G4UniformRand() > 0.25) {*bEnd = Diquark( j1000, j10, 0 );} 241 else {*bEnd = 223 else {*bEnd = Diquark( j1000, j10, 1 );} 242 break; 224 break; 243 } 225 } 244 else 226 else 245 { 227 { 246 if (( j1000 == j100 ) && ( G4UniformRa << 228 if(( j1000 == j100 ) && ( G4UniformRand() > SuppresUUDDSS )) continue; 247 *aEnd = j10; 229 *aEnd = j10; 248 if ( j1000 == j100 ) {*bEnd << 230 if( j1000 == j100 ) {*bEnd = Diquark( j1000, j100, 1 );} 249 else 231 else 250 if ( G4UniformRand() > 0.25) {*bEnd << 232 if( G4UniformRand() > 0.25) {*bEnd = Diquark( j1000, j100, 0 );} 251 else {*bEnd = 233 else {*bEnd = Diquark( j1000, j100, 1 );} 252 break; 234 break; 253 } 235 } 254 } while ( (true) && 236 } while ( (true) && 255 ++loopCounter < maxNumberOfLoops 237 ++loopCounter < maxNumberOfLoops ); /* Loop checking, 10.08.2015, A.Ribon */ 256 if ( loopCounter >= maxNumberOfLoops ) { 238 if ( loopCounter >= maxNumberOfLoops ) { 257 *aEnd = j10; *bEnd = Diquark( j1000, j10 239 *aEnd = j10; *bEnd = Diquark( j1000, j100, 1 ); // Just something acceptable, without any physics consideration. 258 } 240 } >> 241 >> 242 /* >> 243 if ( std::abs( j100 ) >= std::abs( j10 ) ) { >> 244 if ( random < udspin1 ) { >> 245 *aEnd = j1000; >> 246 *bEnd = Diquark( j100, j10, 1 ); >> 247 } else if ( random < udspin1 + uuspin1 ) { >> 248 *aEnd = j10; >> 249 *bEnd = Diquark( j1000, j100, 1 ); >> 250 } else { >> 251 *aEnd = j100; >> 252 // Careful, there is no diquark q1=q2, (q1 q2)0, possible for Omega- >> 253 *bEnd = Diquark( j1000, j10, j1000 != j100 ? 0 : 1 ); >> 254 } >> 255 } else { >> 256 // Lambda-like hadrons have two lightest quarks in spin 0 >> 257 if ( random < udspin1 ) { >> 258 *aEnd = j1000; >> 259 // as above, but with charmed baryons >> 260 *bEnd = Diquark( j100, j10, j100 != j10 ? 0 : 10 ); >> 261 } else if ( random < udspin1 + uuspin1 ) { >> 262 *aEnd = j10; >> 263 *bEnd = Diquark( j1000, j100, 1 ); >> 264 } else { >> 265 *aEnd = j100; >> 266 *bEnd = Diquark( j1000, j10, 1 ); >> 267 } >> 268 } >> 269 */ 259 270 260 } 271 } 261 } 272 } 262 273 263 274 264 //============================================ 275 //============================================================================ 265 276 266 G4int G4DiffractiveSplitableHadron::Diquark( G 277 G4int G4DiffractiveSplitableHadron::Diquark( G4int aquark, G4int bquark, G4int Spin) const { 267 G4int diquarkPDG = std::max( std::abs( aquar 278 G4int diquarkPDG = std::max( std::abs( aquark ), std::abs( bquark ) ) * 1000 + 268 std::min( std::abs( aquar 279 std::min( std::abs( aquark ), std::abs( bquark ) ) * 100 + 269 2*Spin + 1; 280 2*Spin + 1; 270 return ( aquark > 0 && bquark > 0 ) ? diqu 281 return ( aquark > 0 && bquark > 0 ) ? diquarkPDG : -1*diquarkPDG; 271 } 282 } 272 283 273 284