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