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 // 27 // 28 //// ----------------------------------------- 28 //// ------------------------------------------------------------ 29 // GEANT 4 class implementation file 29 // GEANT 4 class implementation file 30 // 30 // 31 // ---------------- G4VPartonStringModel 31 // ---------------- G4VPartonStringModel ---------------- 32 // by Gunter Folger, May 1998. 32 // by Gunter Folger, May 1998. 33 // abstract class for all Parton String M 33 // abstract class for all Parton String Models 34 // ------------------------------------------- 34 // ------------------------------------------------------------ 35 // debug switch 35 // debug switch 36 //#define debug_PartonStringModel 36 //#define debug_PartonStringModel 37 //#define debug_heavyHadrons 37 //#define debug_heavyHadrons 38 // ------------------------------------------- 38 // ------------------------------------------------------------ 39 39 40 #include "G4VPartonStringModel.hh" 40 #include "G4VPartonStringModel.hh" 41 #include "G4ios.hh" 41 #include "G4ios.hh" 42 #include "G4ShortLivedConstructor.hh" 42 #include "G4ShortLivedConstructor.hh" 43 43 44 #include "G4ParticleTable.hh" 44 #include "G4ParticleTable.hh" 45 #include "G4IonTable.hh" 45 #include "G4IonTable.hh" 46 46 47 G4VPartonStringModel::G4VPartonStringModel(con 47 G4VPartonStringModel::G4VPartonStringModel(const G4String& modelName) 48 : G4VHighEnergyGenerator(modelName), 48 : G4VHighEnergyGenerator(modelName), 49 stringFragmentationModel(nullptr) 49 stringFragmentationModel(nullptr) 50 { 50 { 51 // Make shure Shotrylived particles are con 51 // Make shure Shotrylived particles are constructed. 52 // VI: should not instantiate particles by 52 // VI: should not instantiate particles by any model 53 /* 53 /* 54 G4ShortLivedConstructor ShortLived; 54 G4ShortLivedConstructor ShortLived; 55 ShortLived.ConstructParticle(); 55 ShortLived.ConstructParticle(); 56 */ 56 */ 57 } 57 } 58 58 59 G4VPartonStringModel::~G4VPartonStringModel() 59 G4VPartonStringModel::~G4VPartonStringModel() 60 { 60 { 61 } 61 } 62 62 63 G4KineticTrackVector * G4VPartonStringModel::S 63 G4KineticTrackVector * G4VPartonStringModel::Scatter(const G4Nucleus &theNucleus, 64 64 const G4DynamicParticle &aPrimary) 65 { 65 { 66 G4ExcitedStringVector * strings = nullptr; 66 G4ExcitedStringVector * strings = nullptr; 67 G4DynamicParticle thePrimary=aPrimary; 67 G4DynamicParticle thePrimary=aPrimary; 68 G4LorentzVector SumStringMom(0.,0.,0.,0.); 68 G4LorentzVector SumStringMom(0.,0.,0.,0.); 69 G4KineticTrackVector * theResult = 0; 69 G4KineticTrackVector * theResult = 0; 70 G4Nucleon * theNuclNucleon(nullptr); 70 G4Nucleon * theNuclNucleon(nullptr); 71 71 72 #ifdef debug_PartonStringModel 72 #ifdef debug_PartonStringModel 73 G4cout<<G4endl; 73 G4cout<<G4endl; 74 G4cout<<"-----------------------Parton-Strin 74 G4cout<<"-----------------------Parton-String model is runnung ------------"<<G4endl; 75 G4cout<<"Projectile Name Mass "<<thePrimary. 75 G4cout<<"Projectile Name Mass "<<thePrimary.GetDefinition()->GetParticleName()<<" " 76 <<thePrimary. 76 <<thePrimary.GetMass()<<G4endl; 77 G4cout<<" Momentum "<<thePrimary. 77 G4cout<<" Momentum "<<thePrimary.Get4Momentum()<<G4endl; 78 G4cout<<"Target nucleus A Z "<<theNucleus. 78 G4cout<<"Target nucleus A Z "<<theNucleus.GetA_asInt()<<" " 79 <<theNucleus. 79 <<theNucleus.GetZ_asInt()<<G4endl<<G4endl; 80 G4int Bsum=thePrimary.GetDefinition()->GetBa 80 G4int Bsum=thePrimary.GetDefinition()->GetBaryonNumber() + theNucleus.GetA_asInt(); 81 G4int Qsum=thePrimary.GetDefinition()->GetPD 81 G4int Qsum=thePrimary.GetDefinition()->GetPDGCharge() + theNucleus.GetZ_asInt(); 82 G4cout<<"Initial baryon number "<<Bsum<<G4en 82 G4cout<<"Initial baryon number "<<Bsum<<G4endl; 83 G4cout<<"Initial charge "<<Qsum<<G4en 83 G4cout<<"Initial charge "<<Qsum<<G4endl; 84 G4cout<<"-------------- Parton-String model: 84 G4cout<<"-------------- Parton-String model: Generation of strings -------"<<G4endl<<G4endl; 85 Bsum -= theNucleus.GetA_asInt(); Qsum -= th 85 Bsum -= theNucleus.GetA_asInt(); Qsum -= theNucleus.GetZ_asInt(); 86 if(GetProjectileNucleus()) { 86 if(GetProjectileNucleus()) { 87 Bsum -= thePrimary.GetDefinition()->GetBar 87 Bsum -= thePrimary.GetDefinition()->GetBaryonNumber(); 88 Qsum -= thePrimary.GetDefinition()->GetPDG 88 Qsum -= thePrimary.GetDefinition()->GetPDGCharge(); 89 } 89 } 90 G4int QsumSec(0), BsumSec(0); 90 G4int QsumSec(0), BsumSec(0); 91 G4LorentzVector SumPsecondr(0.,0.,0.,0.); 91 G4LorentzVector SumPsecondr(0.,0.,0.,0.); 92 #endif 92 #endif 93 93 94 G4LorentzRotation toZ; 94 G4LorentzRotation toZ; 95 G4LorentzVector Ptmp=thePrimary.Get4Momentum 95 G4LorentzVector Ptmp=thePrimary.Get4Momentum(); 96 toZ.rotateZ(-1*Ptmp.phi()); 96 toZ.rotateZ(-1*Ptmp.phi()); 97 toZ.rotateY(-1*Ptmp.theta()); 97 toZ.rotateY(-1*Ptmp.theta()); 98 thePrimary.Set4Momentum(toZ*Ptmp); 98 thePrimary.Set4Momentum(toZ*Ptmp); 99 G4LorentzRotation toLab(toZ.inverse()); 99 G4LorentzRotation toLab(toZ.inverse()); 100 100 101 G4bool Success=true; 101 G4bool Success=true; 102 G4int attempts = 0, maxAttempts=1000; 102 G4int attempts = 0, maxAttempts=1000; 103 do 103 do 104 { 104 { 105 if (attempts++ > maxAttempts ) 105 if (attempts++ > maxAttempts ) 106 { 106 { 107 Init(theNucleus,thePrimary); // To put 107 Init(theNucleus,thePrimary); // To put a nucleus into ground state 108 / 108 // But marks of hitted nucleons are left. They must be erased. 109 G4V3DNucleus * ResNucleus = GetWoundedNu 109 G4V3DNucleus * ResNucleus = GetWoundedNucleus(); 110 theNuclNucleon = ResNucleus->StartLoop() 110 theNuclNucleon = ResNucleus->StartLoop() ? ResNucleus->GetNextNucleon() : nullptr; 111 while( theNuclNucleon ) 111 while( theNuclNucleon ) 112 { 112 { 113 if(theNuclNucleon->AreYouHit()) theNuc 113 if(theNuclNucleon->AreYouHit()) theNuclNucleon->Hit(nullptr); 114 theNuclNucleon = ResNucleus->GetNextNu 114 theNuclNucleon = ResNucleus->GetNextNucleon(); 115 } 115 } 116 116 117 G4V3DNucleus * ProjResNucleus = GetProje 117 G4V3DNucleus * ProjResNucleus = GetProjectileNucleus(); 118 if(ProjResNucleus != 0) 118 if(ProjResNucleus != 0) 119 { 119 { 120 theNuclNucleon = ProjResNucleus->Start 120 theNuclNucleon = ProjResNucleus->StartLoop() ? ProjResNucleus->GetNextNucleon() : nullptr; 121 while( theNuclNucleon ) 121 while( theNuclNucleon ) 122 { 122 { 123 if(theNuclNucleon->AreYouHit()) theN 123 if(theNuclNucleon->AreYouHit()) theNuclNucleon->Hit(nullptr); 124 theNuclNucleon = ProjResNucleus->Get 124 theNuclNucleon = ProjResNucleus->GetNextNucleon(); 125 } 125 } 126 } 126 } 127 127 128 G4ExceptionDescription ed; 128 G4ExceptionDescription ed; 129 ed << "Projectile Name Mass " <<thePrima 129 ed << "Projectile Name Mass " <<thePrimary.GetDefinition()->GetParticleName() 130 << " " << thePrimary.GetMass()<< G4en 130 << " " << thePrimary.GetMass()<< G4endl; 131 ed << "          Momentum " 131 ed << "          Momentum " << thePrimary.Get4Momentum() <<G4endl; 132 ed << "Target nucleus  A Z " << theNu 132 ed << "Target nucleus  A Z " << theNucleus.GetA_asInt() << " " 133 << theNu 133 << theNucleus.GetZ_asInt() <<G4endl; 134 ed << "Initial states of projectile and 134 ed << "Initial states of projectile and target nucleus will be returned!"<<G4endl; 135 G4Exception( "G4VPartonStringModel::Scat 135 G4Exception( "G4VPartonStringModel::Scatter(): fails to generate or fragment strings ", 136 "HAD_PARTON_STRING_001", Ju 136 "HAD_PARTON_STRING_001", JustWarning, ed ); 137 137 138 G4ThreeVector Position(0.,0.,2*ResNucleu 138 G4ThreeVector Position(0.,0.,2*ResNucleus->GetOuterRadius()); 139 G4KineticTrack* Hadron = new G4KineticTr 139 G4KineticTrack* Hadron = new G4KineticTrack(aPrimary.GetParticleDefinition(), 0., 140 140 Position, aPrimary.Get4Momentum()); 141 if(theResult == nullptr) theResult = new 141 if(theResult == nullptr) theResult = new G4KineticTrackVector(); 142 theResult->push_back(Hadron); 142 theResult->push_back(Hadron); 143 return theResult; 143 return theResult; 144 } 144 } 145 145 146 Success=true; 146 Success=true; 147 147 148 Init(theNucleus,thePrimary); 148 Init(theNucleus,thePrimary); 149 149 150 strings = GetStrings(); 150 strings = GetStrings(); 151 151 152 if (strings->empty()) { Success=false; con 152 if (strings->empty()) { Success=false; continue; } 153 153 154 // G4double stringEnergy(0); 154 // G4double stringEnergy(0); 155 SumStringMom=G4LorentzVector(0.,0.,0.,0.); 155 SumStringMom=G4LorentzVector(0.,0.,0.,0.); 156 156 157 #ifdef debug_PartonStringModel 157 #ifdef debug_PartonStringModel 158 G4cout<<"------------ Parton-String model: 158 G4cout<<"------------ Parton-String model: Number of produced strings ---- "<<strings->size()<<G4endl; 159 #endif 159 #endif 160 160 161 #ifdef debug_heavyHadrons 161 #ifdef debug_heavyHadrons 162 // Check charm and bottom numbers of the p 162 // Check charm and bottom numbers of the projectile: 163 G4int count_charm_projectile = thePrimary 163 G4int count_charm_projectile = thePrimary.GetDefinition()->GetQuarkContent( 4 ) - 164 thePrimary 164 thePrimary.GetDefinition()->GetAntiQuarkContent( 4 ); 165 G4int count_bottom_projectile = thePrimary 165 G4int count_bottom_projectile = thePrimary.GetDefinition()->GetQuarkContent( 5 ) - 166 thePrimary 166 thePrimary.GetDefinition()->GetAntiQuarkContent( 5 ); 167 G4int count_charm_strings = 0, count_botto 167 G4int count_charm_strings = 0, count_bottom_strings = 0; 168 G4int count_charm_hadrons = 0, count_botto 168 G4int count_charm_hadrons = 0, count_bottom_hadrons = 0; 169 #endif 169 #endif 170 170 171 for ( unsigned int astring=0; astring < st 171 for ( unsigned int astring=0; astring < strings->size(); astring++) 172 { 172 { 173 // rotate string to lab frame, models 173 // rotate string to lab frame, models have it aligned to z 174 if((*strings)[astring]->IsExcited()) 174 if((*strings)[astring]->IsExcited()) 175 { 175 { 176 // stringEnergy += (*strings)[astring] 176 // stringEnergy += (*strings)[astring]->GetLeftParton()->Get4Momentum().t(); 177 // stringEnergy += (*strings)[astring] 177 // stringEnergy += (*strings)[astring]->GetRightParton()->Get4Momentum().t(); 178 (*strings)[astring]->LorentzRotate(toL 178 (*strings)[astring]->LorentzRotate(toLab); 179 SumStringMom+=(*strings)[astring]->Get 179 SumStringMom+=(*strings)[astring]->Get4Momentum(); 180 #ifdef debug_PartonStringModel 180 #ifdef debug_PartonStringModel 181 G4cout<<"String No "<<astring+1<<" "<< 181 G4cout<<"String No "<<astring+1<<" "<<(*strings)[astring]->Get4Momentum()<<" " 182 <<(*strings)[astri 182 <<(*strings)[astring]->Get4Momentum().mag() 183 <<" Partons "<<(*strings)[astr 183 <<" Partons "<<(*strings)[astring]->GetLeftParton()->GetDefinition()->GetPDGEncoding() 184 <<" "<<(*strings)[astri 184 <<" "<<(*strings)[astring]->GetRightParton()->GetDefinition()->GetPDGEncoding()<<G4endl; 185 #endif 185 #endif 186 186 187 #ifdef debug_heavyHadrons 187 #ifdef debug_heavyHadrons 188 G4int left_charm = (*strings)[astring]- 188 G4int left_charm = (*strings)[astring]->GetLeftParton()->GetDefinition()->GetQuarkContent( 4 ); 189 G4int left_anticharm = (*strings)[astring]- 189 G4int left_anticharm = (*strings)[astring]->GetLeftParton()->GetDefinition()->GetAntiQuarkContent( 4 ); 190 G4int right_charm = (*strings)[astring]- 190 G4int right_charm = (*strings)[astring]->GetRightParton()->GetDefinition()->GetQuarkContent( 4 ); 191 G4int right_anticharm = (*strings)[astring]- 191 G4int right_anticharm = (*strings)[astring]->GetRightParton()->GetDefinition()->GetAntiQuarkContent( 4 ); 192 G4int left_bottom = (*strings)[astring] 192 G4int left_bottom = (*strings)[astring]->GetLeftParton()->GetDefinition()->GetQuarkContent( 5 ); 193 G4int left_antibottom = (*strings)[astring] 193 G4int left_antibottom = (*strings)[astring]->GetLeftParton()->GetDefinition()->GetAntiQuarkContent( 5 ); 194 G4int right_bottom = (*strings)[astring] 194 G4int right_bottom = (*strings)[astring]->GetRightParton()->GetDefinition()->GetQuarkContent( 5 ); 195 G4int right_antibottom = (*strings)[astring] 195 G4int right_antibottom = (*strings)[astring]->GetRightParton()->GetDefinition()->GetAntiQuarkContent( 5 ); 196 if ( left_charm != 0 || left_anticharm ! 196 if ( left_charm != 0 || left_anticharm != 0 || right_charm != 0 || right_anticharm != 0 || 197 left_bottom != 0 || left_antibottom ! 197 left_bottom != 0 || left_antibottom != 0 || right_bottom != 0 || right_antibottom != 0 ) { 198 count_charm_strings += left_charm - left 198 count_charm_strings += left_charm - left_anticharm + right_charm - right_anticharm; 199 count_bottom_strings += left_bottom - left 199 count_bottom_strings += left_bottom - left_antibottom + right_bottom - right_antibottom; 200 G4cout << "G4VPartonStringModel::Scatter : 200 G4cout << "G4VPartonStringModel::Scatter : string #" << astring << " (" 201 << (*strings)[astring]->GetLeftParton()-> 201 << (*strings)[astring]->GetLeftParton()->GetDefinition()->GetParticleName() << " , " 202 << (*strings)[astring]->GetRightParton()- 202 << (*strings)[astring]->GetRightParton()->GetDefinition()->GetParticleName() << ")" << G4endl; 203 } 203 } 204 #endif 204 #endif 205 } 205 } 206 else 206 else 207 { 207 { 208 // stringEnergy += (*strings)[astring] 208 // stringEnergy += (*strings)[astring]->GetKineticTrack()->Get4Momentum().t(); 209 (*strings)[astring]->LorentzRotate(toL 209 (*strings)[astring]->LorentzRotate(toLab); 210 SumStringMom+=(*strings)[astring]->Get 210 SumStringMom+=(*strings)[astring]->GetKineticTrack()->Get4Momentum(); 211 #ifdef debug_PartonStringModel 211 #ifdef debug_PartonStringModel 212 G4cout<<"A track No "<<astring+1<<" " 212 G4cout<<"A track No "<<astring+1<<" " 213 <<(*strings)[astring]->GetKineti 213 <<(*strings)[astring]->GetKineticTrack()->Get4Momentum()<<" " 214 <<(*strings)[astring]->GetKineti 214 <<(*strings)[astring]->GetKineticTrack()->Get4Momentum().mag()<<" " 215 <<(*strings)[astring]->GetKineti 215 <<(*strings)[astring]->GetKineticTrack()->GetDefinition()->GetParticleName()<<G4endl; 216 #endif 216 #endif 217 217 218 #ifdef debug_heavyHadrons 218 #ifdef debug_heavyHadrons 219 G4int charm = (*strings)[astring]->GetK 219 G4int charm = (*strings)[astring]->GetKineticTrack()->GetDefinition()->GetQuarkContent( 4 ); 220 G4int anticharm = (*strings)[astring]->GetK 220 G4int anticharm = (*strings)[astring]->GetKineticTrack()->GetDefinition()->GetAntiQuarkContent( 4 ); 221 G4int bottom = (*strings)[astring]->GetK 221 G4int bottom = (*strings)[astring]->GetKineticTrack()->GetDefinition()->GetQuarkContent( 5 ); 222 G4int antibottom = (*strings)[astring]->GetK 222 G4int antibottom = (*strings)[astring]->GetKineticTrack()->GetDefinition()->GetAntiQuarkContent( 5 ); 223 if ( charm != 0 || anticharm != 0 | 223 if ( charm != 0 || anticharm != 0 || bottom != 0 || antibottom != 0 ) { 224 count_charm_strings += charm - anticharm 224 count_charm_strings += charm - anticharm; 225 count_bottom_strings += bottom - antibotto 225 count_bottom_strings += bottom - antibottom; 226 G4cout << "G4VPartonStringModel::Scatter : 226 G4cout << "G4VPartonStringModel::Scatter : track #" << astring << "\t" 227 << (*strings)[astring]->GetKi 227 << (*strings)[astring]->GetKineticTrack()->GetDefinition()->GetParticleName() << G4endl; 228 } 228 } 229 #endif 229 #endif 230 } 230 } 231 } 231 } 232 232 233 #ifdef debug_heavyHadrons 233 #ifdef debug_heavyHadrons 234 if ( count_charm_projectile != count_charm 234 if ( count_charm_projectile != count_charm_strings ) { 235 G4cout << "G4VPartonStringModel::Scatter 235 G4cout << "G4VPartonStringModel::Scatter : CHARM VIOLATION in String formation ! #projectile=" 236 << count_charm_projectile << " ; #strin 236 << count_charm_projectile << " ; #strings=" << count_charm_strings << G4endl; 237 } 237 } 238 if ( count_bottom_projectile != count_bott 238 if ( count_bottom_projectile != count_bottom_strings ) { 239 G4cout << "G4VPartonStringModel::Scatter 239 G4cout << "G4VPartonStringModel::Scatter : BOTTOM VIOLATION in String formation ! #projectile=" 240 << count_bottom_projectile << " ; #stri 240 << count_bottom_projectile << " ; #strings=" << count_bottom_strings << G4endl; 241 } 241 } 242 #endif 242 #endif 243 243 244 #ifdef debug_PartonStringModel 244 #ifdef debug_PartonStringModel 245 G4cout<<G4endl<<"SumString4Mom "<<SumStrin 245 G4cout<<G4endl<<"SumString4Mom "<<SumStringMom<<G4endl; 246 G4LorentzVector TargetResidual4Momentum(0. 246 G4LorentzVector TargetResidual4Momentum(0.,0.,0.,0.); 247 G4LorentzVector ProjectileResidual4Momentu 247 G4LorentzVector ProjectileResidual4Momentum(0.,0.,0.,0.); 248 G4int hitsT(0), charged_hitsT(0); 248 G4int hitsT(0), charged_hitsT(0); 249 G4int hitsP(0), charged_hitsP(0); 249 G4int hitsP(0), charged_hitsP(0); 250 G4double ExcitationEt(0.), ExcitationEp(0. 250 G4double ExcitationEt(0.), ExcitationEp(0.); 251 #endif 251 #endif 252 252 253 // We assume that the target nucleus is ne 253 // We assume that the target nucleus is never a hypernucleus, whereas 254 // the projectile nucleus can be a light h 254 // the projectile nucleus can be a light hypernucleus or anti-hypernucleus. 255 255 256 G4V3DNucleus * ProjResNucleus = GetProject 256 G4V3DNucleus * ProjResNucleus = GetProjectileNucleus(); 257 257 258 G4int numberProtonProjectileResidual( 0 ), 258 G4int numberProtonProjectileResidual( 0 ), numberNeutronProjectileResidual( 0 ); 259 G4int numberLambdaProjectileResidual( 0 ); 259 G4int numberLambdaProjectileResidual( 0 ); 260 if(ProjResNucleus != 0) 260 if(ProjResNucleus != 0) 261 { 261 { 262 theNuclNucleon = ProjResNucleus->StartLo 262 theNuclNucleon = ProjResNucleus->StartLoop() ? ProjResNucleus->GetNextNucleon() : nullptr; 263 G4int numberProtonProjectileHits( 0 ), n 263 G4int numberProtonProjectileHits( 0 ), numberNeutronProjectileHits( 0 ); 264 G4int numberLambdaProjectileHits( 0 ); 264 G4int numberLambdaProjectileHits( 0 ); 265 while( theNuclNucleon ) 265 while( theNuclNucleon ) 266 { 266 { 267 if(theNuclNucleon->AreYouHit()) 267 if(theNuclNucleon->AreYouHit()) 268 { 268 { 269 G4LorentzVector tmp=toLab*theNuclNuc 269 G4LorentzVector tmp=toLab*theNuclNucleon->Get4Momentum(); 270 const G4ParticleDefinition* def = theNuclN 270 const G4ParticleDefinition* def = theNuclNucleon->GetDefinition(); 271 #ifdef debug_PartonStringModel 271 #ifdef debug_PartonStringModel 272 ProjectileResidual4Momentum += tmp; 272 ProjectileResidual4Momentum += tmp; 273 hitsP++; 273 hitsP++; 274 if ( def == G4Proton::Definition() | 274 if ( def == G4Proton::Definition() || def == G4AntiProton::Definition() ) ++charged_hitsP; 275 ExcitationEp +=theNuclNucleon->GetBi 275 ExcitationEp +=theNuclNucleon->GetBindingEnergy(); 276 #endif 276 #endif 277 theNuclNucleon->SetMomentum(tmp); 277 theNuclNucleon->SetMomentum(tmp); 278 if ( def == G4Proton::Definition() 278 if ( def == G4Proton::Definition() || def == G4AntiProton::Definition() ) ++numberProtonProjectileHits; 279 if ( def == G4Neutron::Definition() 279 if ( def == G4Neutron::Definition() || def == G4AntiNeutron::Definition() ) ++numberNeutronProjectileHits; 280 if ( def == G4Lambda::Definition() 280 if ( def == G4Lambda::Definition() || def == G4AntiLambda::Definition() ) ++numberLambdaProjectileHits; 281 } 281 } 282 theNuclNucleon = ProjResNucleus->GetNe 282 theNuclNucleon = ProjResNucleus->GetNextNucleon(); 283 } 283 } 284 G4int numberLambdaProjectile = 0; 284 G4int numberLambdaProjectile = 0; 285 if ( thePrimary.GetDefinition()->IsHyper 285 if ( thePrimary.GetDefinition()->IsHypernucleus() ) { 286 numberLambdaProjectile = thePrimary.GetDefin 286 numberLambdaProjectile = thePrimary.GetDefinition()->GetNumberOfLambdasInHypernucleus(); 287 } else if ( thePrimary.GetDefinition()-> 287 } else if ( thePrimary.GetDefinition()->IsAntiHypernucleus() ) { 288 numberLambdaProjectile = thePrimary.GetDefin 288 numberLambdaProjectile = thePrimary.GetDefinition()->GetNumberOfAntiLambdasInAntiHypernucleus(); 289 } 289 } 290 #ifdef debug_PartonStringModel 290 #ifdef debug_PartonStringModel 291 G4cout<<"Projectile residual A, Z (numbe 291 G4cout<<"Projectile residual A, Z (numberOfLambdasOrAntiLambdas) and E* " 292 <<thePrimary.GetDefinition()->GetB 292 <<thePrimary.GetDefinition()->GetBaryonNumber() - hitsP<<" " 293 <<thePrimary.GetDefinition()->GetP 293 <<thePrimary.GetDefinition()->GetPDGCharge() - charged_hitsP<<" (" 294 << numberLambdaProjectile - numberLambda 294 << numberLambdaProjectile - numberLambdaProjectileHits << ") " 295 <<ExcitationEp<<G4endl; 295 <<ExcitationEp<<G4endl; 296 G4cout<<"Projectile residual 4 momentum 296 G4cout<<"Projectile residual 4 momentum "<<ProjectileResidual4Momentum<<G4endl; 297 #endif 297 #endif 298 numberProtonProjectileResidual = std::ma 298 numberProtonProjectileResidual = std::max( std::abs( G4int( thePrimary.GetDefinition()->GetPDGCharge() ) ) - 299 numberProtonProjectileHits, 0 ); 299 numberProtonProjectileHits, 0 ); 300 numberLambdaProjectileResidual = std::ma 300 numberLambdaProjectileResidual = std::max( numberLambdaProjectile - numberLambdaProjectileHits, 0 ); 301 numberNeutronProjectileResidual = std::m 301 numberNeutronProjectileResidual = std::max( std::abs( thePrimary.GetDefinition()->GetBaryonNumber() ) - 302 302 std::abs( G4int( thePrimary.GetDefinition()->GetPDGCharge() ) ) - 303 numberLambdaProjectile - numberN 303 numberLambdaProjectile - numberNeutronProjectileHits, 0 ); 304 } 304 } 305 305 306 G4V3DNucleus * ResNucleus = GetWoundedNucl 306 G4V3DNucleus * ResNucleus = GetWoundedNucleus(); 307 307 308 // loop over wounded nucleus 308 // loop over wounded nucleus 309 theNuclNucleon = ResNucleus->StartLoop() ? 309 theNuclNucleon = ResNucleus->StartLoop() ? ResNucleus->GetNextNucleon() : nullptr; 310 G4int numberProtonTargetHits( 0 ), numberN 310 G4int numberProtonTargetHits( 0 ), numberNeutronTargetHits( 0 ); 311 while( theNuclNucleon ) 311 while( theNuclNucleon ) 312 { 312 { 313 if(theNuclNucleon->AreYouHit()) 313 if(theNuclNucleon->AreYouHit()) 314 { 314 { 315 G4LorentzVector tmp=toLab*theNuclNucle 315 G4LorentzVector tmp=toLab*theNuclNucleon->Get4Momentum(); 316 #ifdef debug_PartonStringModel 316 #ifdef debug_PartonStringModel 317 TargetResidual4Momentum += tmp; 317 TargetResidual4Momentum += tmp; 318 hitsT++; 318 hitsT++; 319 if ( theNuclNucleon->GetDefinition() = 319 if ( theNuclNucleon->GetDefinition() == G4Proton::Proton() ) ++charged_hitsT; 320 ExcitationEt +=theNuclNucleon->GetBind 320 ExcitationEt +=theNuclNucleon->GetBindingEnergy(); 321 #endif 321 #endif 322 theNuclNucleon->SetMomentum(tmp); 322 theNuclNucleon->SetMomentum(tmp); 323 if ( theNuclNucleon->GetDefinition() = 323 if ( theNuclNucleon->GetDefinition() == G4Proton::Proton() ) ++numberProtonTargetHits; 324 if ( theNuclNucleon->GetDefinition() = 324 if ( theNuclNucleon->GetDefinition() == G4Neutron::Neutron() ) ++numberNeutronTargetHits; 325 } 325 } 326 theNuclNucleon = ResNucleus->GetNextNucl 326 theNuclNucleon = ResNucleus->GetNextNucleon(); 327 } 327 } 328 328 329 #ifdef debug_PartonStringModel 329 #ifdef debug_PartonStringModel 330 G4cout<<"Target residual A, Z and E* " 330 G4cout<<"Target residual A, Z and E* " 331 <<theNucleus.GetA_asInt() - hitsT<<" 331 <<theNucleus.GetA_asInt() - hitsT<<" " 332 <<theNucleus.GetZ_asInt() - charged_ 332 <<theNucleus.GetZ_asInt() - charged_hitsT<<" " 333 <<ExcitationEt<<G4endl; 333 <<ExcitationEt<<G4endl; 334 G4cout<<"Target residual 4 momentum "<<Tar 334 G4cout<<"Target residual 4 momentum "<<TargetResidual4Momentum<<G4endl; 335 Bsum+=( hitsT + hitsP); 335 Bsum+=( hitsT + hitsP); 336 Qsum+=(charged_hitsT + charged_hitsP); 336 Qsum+=(charged_hitsT + charged_hitsP); 337 G4cout<<"Hitted # of nucleons of projectil 337 G4cout<<"Hitted # of nucleons of projectile and target "<<hitsP<<" "<<hitsT<<G4endl; 338 G4cout<<"Hitted # of protons of projectile 338 G4cout<<"Hitted # of protons of projectile and target " 339 <<charged_hitsP<<" "<<charged_hitsT< 339 <<charged_hitsP<<" "<<charged_hitsT<<G4endl<<G4endl; 340 G4cout<<"Bsum Qsum "<<Bsum<<" "<<Qsum<<G4e 340 G4cout<<"Bsum Qsum "<<Bsum<<" "<<Qsum<<G4endl<<G4endl; 341 #endif 341 #endif 342 342 343 // Re-sample in the case of unphysical nuc 343 // Re-sample in the case of unphysical nuclear residual: 344 // 1 (H), 2 (2He), and 3 (3Li) protons alo 344 // 1 (H), 2 (2He), and 3 (3Li) protons alone without neutrons can exist, but not more; 345 // no bound states of 2 or more neutrons w 345 // no bound states of 2 or more neutrons without protons can exist. 346 G4int numberProtonTargetResidual = theNucl 346 G4int numberProtonTargetResidual = theNucleus.GetZ_asInt() - numberProtonTargetHits; 347 G4int numberNeutronTargetResidual = theNuc 347 G4int numberNeutronTargetResidual = theNucleus.GetA_asInt() - theNucleus.GetZ_asInt() - numberNeutronTargetHits; 348 G4bool unphysicalResidual = false; 348 G4bool unphysicalResidual = false; 349 if ( ( numberProtonTargetResidual > 3 && 349 if ( ( numberProtonTargetResidual > 3 && numberNeutronTargetResidual == 0 ) || 350 ( numberProtonTargetResidual == 0 && 350 ( numberProtonTargetResidual == 0 && numberNeutronTargetResidual > 1 ) ) { 351 unphysicalResidual = true; 351 unphysicalResidual = true; 352 //G4cout << "***UNPHYSICAL TARGET RESIDU 352 //G4cout << "***UNPHYSICAL TARGET RESIDUAL*** Z=" << numberProtonTargetResidual 353 // << " ; N=" << numberNeutronTarg 353 // << " ; N=" << numberNeutronTargetResidual; 354 } 354 } 355 // The projectile residual can be a hypern 355 // The projectile residual can be a hypernucleus or anti-hypernucleus: 356 // only the following combinations are cur 356 // only the following combinations are currently allowed in Geant4: 357 // p-n-lambda (hypertriton), p-n-n-lambda 357 // p-n-lambda (hypertriton), p-n-n-lambda (hyperH4), p-p-n-lambda (hyperAlpha), 358 // p-p-n-n-lambda (hyperHe5), n-n-lambda-l 358 // p-p-n-n-lambda (hyperHe5), n-n-lambda-lambda (doublehyperdoubleneutron), 359 // p-n-lambda-lambda (doubleHyperH4) 359 // p-n-lambda-lambda (doubleHyperH4) 360 if ( ( numberProtonProjectileResidual > 3 360 if ( ( numberProtonProjectileResidual > 3 && numberNeutronProjectileResidual == 0 ) || 361 ( numberProtonProjectileResidual == 0 361 ( numberProtonProjectileResidual == 0 && numberNeutronProjectileResidual > 1 && 362 numberLambdaProjectileResidual == 0 ) || 362 numberLambdaProjectileResidual == 0 ) || 363 ( numberProtonProjectileResidual == 0 && 363 ( numberProtonProjectileResidual == 0 && numberNeutronProjectileResidual <= 1 && 364 numberLambdaProjectileResidual > 0 ) || 364 numberLambdaProjectileResidual > 0 ) || 365 ( numberProtonProjectileResidual == 0 && 365 ( numberProtonProjectileResidual == 0 && numberNeutronProjectileResidual > 2 && 366 numberLambdaProjectileResidual > 0 ) || 366 numberLambdaProjectileResidual > 0 ) || 367 ( numberLambdaProjectileResidual > 2 ) || 367 ( numberLambdaProjectileResidual > 2 ) || 368 ( numberProtonProjectileResidual > 0 && n 368 ( numberProtonProjectileResidual > 0 && numberNeutronProjectileResidual == 0 && 369 numberLambdaProjectileResidual > 0 ) || 369 numberLambdaProjectileResidual > 0 ) || 370 ( numberProtonProjectileResidual > 1 && n 370 ( numberProtonProjectileResidual > 1 && numberNeutronProjectileResidual > 1 && 371 numberLambdaProjectileResidual > 1 ) 371 numberLambdaProjectileResidual > 1 ) 372 ) { 372 ) { 373 unphysicalResidual = true; 373 unphysicalResidual = true; 374 //G4cout << "***UNPHYSICAL PROJECTILE RE 374 //G4cout << "***UNPHYSICAL PROJECTILE RESIDUAL*** Z=" << numberProtonProjectileResidual 375 // << " ; N=" << numberNeutronProj 375 // << " ; N=" << numberNeutronProjectileResidual 376 // << " ; L=" << numberLambdaProje 376 // << " ; L=" << numberLambdaProjectileResidual; 377 } 377 } 378 if ( unphysicalResidual ) { 378 if ( unphysicalResidual ) { 379 //G4cout << " -> REJECTING COLLISION bec 379 //G4cout << " -> REJECTING COLLISION because of unphysical residual !" << G4endl; 380 Success = false; 380 Success = false; 381 continue; 381 continue; 382 } 382 } 383 383 384 //======================================== 384 //========================================================================================= 385 // Fragment s 385 // Fragment strings 386 #ifdef debug_PartonStringModel 386 #ifdef debug_PartonStringModel 387 G4cout<<"---------------- Attempt to fragm 387 G4cout<<"---------------- Attempt to fragment strings ------------- "<<attempts<<G4endl; 388 #endif 388 #endif 389 389 390 G4double InvMass=SumStringMom.mag(); 390 G4double InvMass=SumStringMom.mag(); 391 G4double SumMass(0.); 391 G4double SumMass(0.); 392 392 393 #ifdef debug_PartonStringModel 393 #ifdef debug_PartonStringModel 394 QsumSec=0; BsumSec=0; 394 QsumSec=0; BsumSec=0; 395 SumPsecondr=G4LorentzVector(0.,0.,0.,0.); 395 SumPsecondr=G4LorentzVector(0.,0.,0.,0.); 396 #endif 396 #endif 397 397 398 if(theResult != nullptr) 398 if(theResult != nullptr) 399 { 399 { 400 std::for_each(theResult->begin(), theRes 400 std::for_each(theResult->begin(), theResult->end(), DeleteKineticTrack()); 401 delete theResult; 401 delete theResult; 402 } 402 } 403 403 404 theResult = stringFragmentationModel->Frag 404 theResult = stringFragmentationModel->FragmentStrings(strings); 405 405 406 #ifdef debug_PartonStringModel 406 #ifdef debug_PartonStringModel 407 G4cout<<"String fragmentation success (OK 407 G4cout<<"String fragmentation success (OK - #0, Not OK - 0) : "<<theResult<<G4endl; 408 #endif 408 #endif 409 409 410 if(theResult == 0) {Success=false; continu 410 if(theResult == 0) {Success=false; continue;} 411 411 412 #ifdef debug_PartonStringModel 412 #ifdef debug_PartonStringModel 413 G4cout<<"Parton-String model: Number of pr 413 G4cout<<"Parton-String model: Number of produced particles "<<theResult->size()<<G4endl; 414 SumPsecondr = G4LorentzVector(0.,0.,0.,0.) 414 SumPsecondr = G4LorentzVector(0.,0.,0.,0.); 415 QsumSec = 0; BsumSec = 0; 415 QsumSec = 0; BsumSec = 0; 416 #endif 416 #endif 417 417 418 SumMass=0.; 418 SumMass=0.; 419 for ( unsigned int i=0; i < theResult->siz 419 for ( unsigned int i=0; i < theResult->size(); i++) 420 { 420 { 421 SumMass+=(*theResult)[i]->Get4Momentum() 421 SumMass+=(*theResult)[i]->Get4Momentum().mag(); 422 #ifdef debug_PartonStringModel 422 #ifdef debug_PartonStringModel 423 G4cout<<i<<" : "<<(*theResult)[i]->GetDe 423 G4cout<<i<<" : "<<(*theResult)[i]->GetDefinition()->GetParticleName()<<" " 424 <<(*theResult)[i]->Get4M 424 <<(*theResult)[i]->Get4Momentum()<<" " 425 <<(*theResult)[i]->Get4M 425 <<(*theResult)[i]->Get4Momentum().mag()<<" " 426 <<(*theResult)[i]->GetDe 426 <<(*theResult)[i]->GetDefinition()->GetPDGMass()<<G4endl; 427 SumPsecondr+=(*theResult)[i]->Get4Moment 427 SumPsecondr+=(*theResult)[i]->Get4Momentum(); 428 BsumSec += (*theResult)[i]->GetDefinitio 428 BsumSec += (*theResult)[i]->GetDefinition()->GetBaryonNumber(); 429 QsumSec += (*theResult)[i]->GetDefinitio 429 QsumSec += (*theResult)[i]->GetDefinition()->GetPDGCharge(); 430 #endif 430 #endif 431 431 432 #ifdef debug_heavyHadrons 432 #ifdef debug_heavyHadrons 433 G4int charm = (*theResult)[i]->GetD 433 G4int charm = (*theResult)[i]->GetDefinition()->GetQuarkContent( 4 ); 434 G4int anticharm = (*theResult)[i]->GetD 434 G4int anticharm = (*theResult)[i]->GetDefinition()->GetAntiQuarkContent( 4 ); 435 G4int bottom = (*theResult)[i]->GetD 435 G4int bottom = (*theResult)[i]->GetDefinition()->GetQuarkContent( 5 ); 436 G4int antibottom = (*theResult)[i]->GetD 436 G4int antibottom = (*theResult)[i]->GetDefinition()->GetAntiQuarkContent( 5 ); 437 if ( charm != 0 || anticharm != 0 || 437 if ( charm != 0 || anticharm != 0 || bottom != 0 || antibottom != 0 ) { 438 count_charm_hadrons += charm - antich 438 count_charm_hadrons += charm - anticharm; 439 count_bottom_hadrons += bottom - antib 439 count_bottom_hadrons += bottom - antibottom; 440 G4cout << "G4VPartonStringModel::Scatter : h 440 G4cout << "G4VPartonStringModel::Scatter : hadron #" << i << "\t" 441 << (*theResult)[i]->GetDefiniti 441 << (*theResult)[i]->GetDefinition()->GetParticleName() << G4endl; 442 } 442 } 443 #endif 443 #endif 444 } 444 } 445 445 446 #ifdef debug_heavyHadrons 446 #ifdef debug_heavyHadrons 447 if ( count_charm_projectile != count_charm 447 if ( count_charm_projectile != count_charm_hadrons ) { 448 G4cout << "G4VPartonStringModel::Scatter 448 G4cout << "G4VPartonStringModel::Scatter : CHARM VIOLATION in String hadronization ! #projectile=" 449 << count_charm_projectile << " ; #hadro 449 << count_charm_projectile << " ; #hadrons=" << count_charm_hadrons << G4endl; 450 } 450 } 451 if ( count_bottom_projectile != count_bott 451 if ( count_bottom_projectile != count_bottom_hadrons ) { 452 G4cout << "G4VPartonStringModel::Scatter 452 G4cout << "G4VPartonStringModel::Scatter : BOTTOM VIOLATION in String hadronization ! #projectile=" 453 << count_bottom_projectile << " ; #hadr 453 << count_bottom_projectile << " ; #hadrons=" << count_bottom_hadrons << G4endl; 454 } 454 } 455 #endif 455 #endif 456 456 457 #ifdef debug_PartonStringModel 457 #ifdef debug_PartonStringModel 458 G4cout<<G4endl<<"-----------------------Pa 458 G4cout<<G4endl<<"-----------------------Parton-String model: balances -------------"<<G4endl; 459 if(Qsum != QsumSec) { 459 if(Qsum != QsumSec) { 460 G4cout<<"Charge is not conserved!!! ---- 460 G4cout<<"Charge is not conserved!!! ----"<<G4endl; 461 G4cout<<" Qsum != QsumSec "<<Qsum<<" "<< 461 G4cout<<" Qsum != QsumSec "<<Qsum<<" "<<QsumSec<<G4endl; 462 } 462 } 463 if(Bsum != BsumSec) { 463 if(Bsum != BsumSec) { 464 G4cout<<"Baryon number is not conserved! 464 G4cout<<"Baryon number is not conserved!!!"<<G4endl; 465 G4cout<<" Bsum != BsumSec "<<Bsum<<" "<< 465 G4cout<<" Bsum != BsumSec "<<Bsum<<" "<<BsumSec<<G4endl; 466 } 466 } 467 #endif 467 #endif 468 468 469 if((SumMass > InvMass)||(SumMass == 0.)) { 469 if((SumMass > InvMass)||(SumMass == 0.)) {Success=false;} 470 std::for_each(strings->begin(), strings->e 470 std::for_each(strings->begin(), strings->end(), DeleteString() ); 471 delete strings; 471 delete strings; 472 472 473 } while(!Success); 473 } while(!Success); 474 474 475 #ifdef debug_PartonStringModel 475 #ifdef debug_PartonStringModel 476 G4cout <<"Baryon number balance "<<Bs 476 G4cout <<"Baryon number balance "<<Bsum-BsumSec<<G4endl; 477 G4cout <<"Charge balance "<<Qs 477 G4cout <<"Charge balance "<<Qsum-QsumSec<<G4endl; 478 G4cout <<"4 momentum balance "<<Su 478 G4cout <<"4 momentum balance "<<SumStringMom-SumPsecondr<<G4endl; 479 G4cout<<"---------------------End of Parton- 479 G4cout<<"---------------------End of Parton-String model work -------------"<<G4endl<<G4endl; 480 #endif 480 #endif 481 481 482 return theResult; 482 return theResult; 483 } 483 } 484 484 485 void G4VPartonStringModel::ModelDescription(st 485 void G4VPartonStringModel::ModelDescription(std::ostream& outFile) const 486 { 486 { 487 outFile << GetModelName() << " has no descri 487 outFile << GetModelName() << " has no description yet.\n"; 488 } 488 } 489 489 490 G4V3DNucleus * G4VPartonStringModel::GetProjec 490 G4V3DNucleus * G4VPartonStringModel::GetProjectileNucleus() const 491 { return nullptr; } 491 { return nullptr; } 492 492 493 493