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 #include "globals.hh" 27 #include "globals.hh" 28 #include "G4PhysicalConstants.hh" 28 #include "G4PhysicalConstants.hh" 29 #include "G4SystemOfUnits.hh" 29 #include "G4SystemOfUnits.hh" 30 #include "G4Proton.hh" 30 #include "G4Proton.hh" 31 #include "G4Neutron.hh" 31 #include "G4Neutron.hh" 32 #include "G4LorentzRotation.hh" 32 #include "G4LorentzRotation.hh" 33 #include "G4BinaryCascade.hh" 33 #include "G4BinaryCascade.hh" 34 #include "G4KineticTrackVector.hh" 34 #include "G4KineticTrackVector.hh" 35 #include "G4DecayKineticTracks.hh" 35 #include "G4DecayKineticTracks.hh" 36 #include "G4ReactionProductVector.hh" 36 #include "G4ReactionProductVector.hh" 37 #include "G4Track.hh" 37 #include "G4Track.hh" 38 #include "G4V3DNucleus.hh" 38 #include "G4V3DNucleus.hh" 39 #include "G4Fancy3DNucleus.hh" 39 #include "G4Fancy3DNucleus.hh" 40 #include "G4Scatterer.hh" 40 #include "G4Scatterer.hh" 41 #include "G4MesonAbsorption.hh" 41 #include "G4MesonAbsorption.hh" 42 #include "G4ping.hh" 42 #include "G4ping.hh" 43 #include "G4Delete.hh" 43 #include "G4Delete.hh" 44 44 45 #include "G4CollisionManager.hh" 45 #include "G4CollisionManager.hh" 46 #include "G4Absorber.hh" 46 #include "G4Absorber.hh" 47 47 48 #include "G4CollisionInitialState.hh" 48 #include "G4CollisionInitialState.hh" 49 #include "G4ListOfCollisions.hh" 49 #include "G4ListOfCollisions.hh" 50 #include "G4Fragment.hh" 50 #include "G4Fragment.hh" 51 #include "G4RKPropagation.hh" 51 #include "G4RKPropagation.hh" 52 52 53 #include "G4NuclearShellModelDensity.hh" 53 #include "G4NuclearShellModelDensity.hh" 54 #include "G4NuclearFermiDensity.hh" 54 #include "G4NuclearFermiDensity.hh" 55 #include "G4FermiMomentum.hh" 55 #include "G4FermiMomentum.hh" 56 56 57 #include "G4PreCompoundModel.hh" 57 #include "G4PreCompoundModel.hh" 58 #include "G4ExcitationHandler.hh" 58 #include "G4ExcitationHandler.hh" 59 #include "G4HadronicInteractionRegistry.hh" 59 #include "G4HadronicInteractionRegistry.hh" 60 60 61 #include "G4FermiPhaseSpaceDecay.hh" 61 #include "G4FermiPhaseSpaceDecay.hh" 62 62 63 #include "G4PreCompoundModel.hh" 63 #include "G4PreCompoundModel.hh" 64 #include "G4HadronicParameters.hh" << 65 64 66 #include <algorithm> 65 #include <algorithm> 67 #include "G4ShortLivedConstructor.hh" 66 #include "G4ShortLivedConstructor.hh" 68 #include <typeinfo> 67 #include <typeinfo> 69 68 70 #include "G4PhysicsModelCatalog.hh" << 71 << 72 // turn on general debugging info, and consi 69 // turn on general debugging info, and consistency checks 73 << 74 //#define debug_G4BinaryCascade 1 70 //#define debug_G4BinaryCascade 1 75 71 76 // more detailed debugging -- deprecated 72 // more detailed debugging -- deprecated 77 //#define debug_H1_BinaryCascade 1 73 //#define debug_H1_BinaryCascade 1 78 74 79 // specific debugging info per method or func 75 // specific debugging info per method or functionality 80 //#define debug_BIC_ApplyCollision 1 76 //#define debug_BIC_ApplyCollision 1 81 //#define debug_BIC_CheckPauli 1 77 //#define debug_BIC_CheckPauli 1 82 //#define debug_BIC_CorrectFinalPandE 1 78 //#define debug_BIC_CorrectFinalPandE 1 83 //#define debug_BIC_Propagate 1 79 //#define debug_BIC_Propagate 1 84 //#define debug_BIC_Propagate_Excitation 1 80 //#define debug_BIC_Propagate_Excitation 1 85 //#define debug_BIC_Propagate_Collisions 1 81 //#define debug_BIC_Propagate_Collisions 1 86 //#define debug_BIC_Propagate_finals 1 82 //#define debug_BIC_Propagate_finals 1 87 //#define debug_BIC_DoTimeStep 1 83 //#define debug_BIC_DoTimeStep 1 88 //#define debug_BIC_CorrectBarionsOnBoundary 1 84 //#define debug_BIC_CorrectBarionsOnBoundary 1 89 //#define debug_BIC_GetExcitationEnergy 1 85 //#define debug_BIC_GetExcitationEnergy 1 90 //#define debug_BIC_DeexcitationProducts 1 << 91 //#define debug_BIC_FinalNucleusMomentum 1 86 //#define debug_BIC_FinalNucleusMomentum 1 92 //#define debug_BIC_Final4Momentum 1 87 //#define debug_BIC_Final4Momentum 1 93 //#define debug_BIC_FillVoidnucleus 1 << 94 //#define debug_BIC_FindFragments 1 88 //#define debug_BIC_FindFragments 1 95 //#define debug_BIC_BuildTargetList 1 89 //#define debug_BIC_BuildTargetList 1 96 //#define debug_BIC_FindCollision 1 90 //#define debug_BIC_FindCollision 1 97 //#define debug_BIC_return 1 91 //#define debug_BIC_return 1 98 92 99 //------- << 93 #if defined(debug_G4BinaryCascade) 100 //#if defined(debug_G4BinaryCascade) << 101 #if 0 << 102 #define _CheckChargeAndBaryonNumber_(val) Ch 94 #define _CheckChargeAndBaryonNumber_(val) CheckChargeAndBaryonNumber(val) 103 //#define debugCheckChargeAndBaryonNumberver << 104 #else 95 #else 105 #define _CheckChargeAndBaryonNumber_(val) 96 #define _CheckChargeAndBaryonNumber_(val) 106 #endif 97 #endif 107 //#if defined(debug_G4BinaryCascade) << 108 #if 0 << 109 #define _DebugEpConservation(val) DebugEpCo << 110 //#define debugCheckChargeAndBaryonNumberver << 111 #else << 112 #define _DebugEpConservation(val) << 113 #endif << 114 << 115 G4int G4BinaryCascade::theBIC_ID = -1; << 116 << 117 // 98 // 118 // C O N S T R U C T O R S A N D D E S T 99 // C O N S T R U C T O R S A N D D E S T R U C T O R S 119 // 100 // >> 101 120 G4BinaryCascade::G4BinaryCascade(G4VPreCompoun 102 G4BinaryCascade::G4BinaryCascade(G4VPreCompoundModel* ptr) : 121 G4VIntraNuclearTransportModel("Binary Cascade" 103 G4VIntraNuclearTransportModel("Binary Cascade", ptr) 122 { 104 { 123 // initialise the resonance sector 105 // initialise the resonance sector 124 G4ShortLivedConstructor ShortLived; 106 G4ShortLivedConstructor ShortLived; 125 ShortLived.ConstructParticle(); 107 ShortLived.ConstructParticle(); 126 108 127 theCollisionMgr = new G4CollisionManager; 109 theCollisionMgr = new G4CollisionManager; 128 theDecay=new G4BCDecay; 110 theDecay=new G4BCDecay; 129 theImR.push_back(theDecay); 111 theImR.push_back(theDecay); 130 theLateParticle= new G4BCLateParticle; 112 theLateParticle= new G4BCLateParticle; 131 G4MesonAbsorption * aAb=new G4MesonAbsorpt 113 G4MesonAbsorption * aAb=new G4MesonAbsorption; 132 theImR.push_back(aAb); 114 theImR.push_back(aAb); 133 G4Scatterer * aSc=new G4Scatterer; 115 G4Scatterer * aSc=new G4Scatterer; 134 theH1Scatterer = new G4Scatterer; 116 theH1Scatterer = new G4Scatterer; 135 theImR.push_back(aSc); 117 theImR.push_back(aSc); 136 118 137 thePropagator = new G4RKPropagation; 119 thePropagator = new G4RKPropagation; 138 theCurrentTime = 0.; 120 theCurrentTime = 0.; 139 theBCminP = 45*MeV; 121 theBCminP = 45*MeV; 140 theCutOnP = 90*MeV; 122 theCutOnP = 90*MeV; 141 theCutOnPAbsorb= 0*MeV; // No Absorption 123 theCutOnPAbsorb= 0*MeV; // No Absorption of slow Mesons, other than above G4MesonAbsorption 142 124 143 // reuse existing pre-compound model 125 // reuse existing pre-compound model 144 if(!ptr) { 126 if(!ptr) { 145 G4HadronicInteraction* p = 127 G4HadronicInteraction* p = 146 G4HadronicInteractionRegistry::Instance()->F 128 G4HadronicInteractionRegistry::Instance()->FindModel("PRECO"); 147 G4VPreCompoundModel* pre = static_cast<G 129 G4VPreCompoundModel* pre = static_cast<G4VPreCompoundModel*>(p); 148 if(!pre) { pre = new G4PreCompoundModel( 130 if(!pre) { pre = new G4PreCompoundModel(); } 149 SetDeExcitation(pre); 131 SetDeExcitation(pre); 150 } 132 } 151 theExcitationHandler = GetDeExcitation()-> 133 theExcitationHandler = GetDeExcitation()->GetExcitationHandler(); 152 SetMinEnergy(0.0*GeV); 134 SetMinEnergy(0.0*GeV); 153 SetMaxEnergy(10.1*GeV); 135 SetMaxEnergy(10.1*GeV); 154 //PrintWelcomeMessage(); 136 //PrintWelcomeMessage(); 155 thePrimaryEscape = true; 137 thePrimaryEscape = true; 156 thePrimaryType = 0; 138 thePrimaryType = 0; 157 139 158 SetEnergyMomentumCheckLevels(1.0*perCent, << 140 SetEnergyMomentumCheckLevels(1*perCent, 1*MeV); 159 141 160 // init data members 142 // init data members 161 currentA=currentZ=0; 143 currentA=currentZ=0; 162 lateA=lateZ=0; 144 lateA=lateZ=0; 163 initialA=initialZ=0; 145 initialA=initialZ=0; 164 projectileA=projectileZ=0; 146 projectileA=projectileZ=0; 165 currentInitialEnergy=initial_nuclear_mass= 147 currentInitialEnergy=initial_nuclear_mass=0.; 166 massInNucleus=0.; 148 massInNucleus=0.; 167 theOuterRadius=0.; 149 theOuterRadius=0.; 168 theBIC_ID = G4PhysicsModelCatalog::GetMode << 169 fBCDEBUG = G4HadronicParameters::Instance( << 170 } 150 } 171 151 >> 152 /* >> 153 G4BinaryCascade::G4BinaryCascade(const G4BinaryCascade& ) >> 154 : G4VIntraNuclearTransportModel("Binary Cascade") >> 155 { >> 156 } >> 157 */ >> 158 172 G4BinaryCascade::~G4BinaryCascade() 159 G4BinaryCascade::~G4BinaryCascade() 173 { 160 { 174 ClearAndDestroy(&theTargetList); 161 ClearAndDestroy(&theTargetList); 175 ClearAndDestroy(&theSecondaryList); 162 ClearAndDestroy(&theSecondaryList); 176 ClearAndDestroy(&theCapturedList); 163 ClearAndDestroy(&theCapturedList); 177 delete thePropagator; 164 delete thePropagator; 178 delete theCollisionMgr; 165 delete theCollisionMgr; 179 for(auto & ptr : theImR) { delete ptr; } << 166 std::for_each(theImR.begin(), theImR.end(), Delete<G4BCAction>()); 180 theImR.clear(); << 181 delete theLateParticle; 167 delete theLateParticle; >> 168 //delete theExcitationHandler; 182 delete theH1Scatterer; 169 delete theH1Scatterer; 183 } 170 } 184 171 185 void G4BinaryCascade::ModelDescription(std::os 172 void G4BinaryCascade::ModelDescription(std::ostream& outFile) const 186 { 173 { 187 outFile << "G4BinaryCascade is an intra-nu 174 outFile << "G4BinaryCascade is an intra-nuclear cascade model in which\n" 188 << "an incident hadron collides wi 175 << "an incident hadron collides with a nucleon, forming two\n" 189 << "final-state particles, one or 176 << "final-state particles, one or both of which may be resonances.\n" 190 << "The resonances then decay hadr 177 << "The resonances then decay hadronically and the decay products\n" 191 << "are then propagated through th 178 << "are then propagated through the nuclear potential along curved\n" 192 << "trajectories until they re-int 179 << "trajectories until they re-interact or leave the nucleus.\n" 193 << "This model is valid for incide 180 << "This model is valid for incident pions up to 1.5 GeV and\n" 194 << "nucleons up to 10 GeV.\n" << 181 << "nucleons up to 10 GeV.\n"; 195 << "The remaining excited nucleus << 196 if (theDeExcitation) << 197 { << 198 outFile << theDeExcitation->GetM << 199 theDeExcitation->DeExciteModelDe << 200 } << 201 else if (theExcitationHandler) << 202 { << 203 outFile << "G4ExcitationHandler << 204 theExcitationHandler->ModelDesc << 205 } << 206 else << 207 { << 208 outFile << "void.\n"; << 209 } << 210 outFile<< " \n"; << 211 } 182 } 212 void G4BinaryCascade::PropagateModelDescriptio 183 void G4BinaryCascade::PropagateModelDescription(std::ostream& outFile) const 213 { 184 { 214 outFile << "G4BinaryCascade propagtes seco 185 outFile << "G4BinaryCascade propagtes secondaries produced by a high\n" 215 << "energy model through the wound 186 << "energy model through the wounded nucleus.\n" 216 << "Secondaries are followed after 187 << "Secondaries are followed after the formation time and if\n" 217 << "within the nucleus are propaga 188 << "within the nucleus are propagated through the nuclear\n" 218 << "potential along curved traject 189 << "potential along curved trajectories until they interact\n" 219 << "with a nucleon, decay, or leav 190 << "with a nucleon, decay, or leave the nucleus.\n" 220 << "An interaction of a secondary 191 << "An interaction of a secondary with a nucleon produces two\n" 221 << "final-state particles, one or 192 << "final-state particles, one or both of which may be resonances.\n" 222 << "Resonances decay hadronically 193 << "Resonances decay hadronically and the decay products\n" 223 << "are in turn propagated through 194 << "are in turn propagated through the nuclear potential along curved\n" 224 << "trajectories until they re-int 195 << "trajectories until they re-interact or leave the nucleus.\n" 225 << "This model is valid for pions 196 << "This model is valid for pions up to 1.5 GeV and\n" 226 << "nucleons up to about 3.5 GeV.\ << 197 << "nucleons up to about 3.5 GeV.\n"; 227 << "The remaining excited nucleus << 228 if (theDeExcitation) // pre << 229 { << 230 outFile << theDeExcitation->GetModelName << 231 theDeExcitation->DeExciteModelDescriptio << 232 } << 233 else if (theExcitationHandler) // de-ex << 234 { << 235 outFile << "G4ExcitationHandler"; // << 236 theExcitationHandler->ModelDescription( << 237 } << 238 else << 239 { << 240 outFile << "void.\n"; << 241 } << 242 outFile<< " \n"; << 243 } 198 } 244 199 245 //-------------------------------------------- 200 //---------------------------------------------------------------------------- 246 201 247 // 202 // 248 // I M P L E M E N T A T I O N 203 // I M P L E M E N T A T I O N 249 // 204 // 250 205 251 206 252 //-------------------------------------------- 207 //---------------------------------------------------------------------------- 253 G4HadFinalState * G4BinaryCascade::ApplyYourse 208 G4HadFinalState * G4BinaryCascade::ApplyYourself(const G4HadProjectile & aTrack, 254 G4Nucleus & aNucleus) 209 G4Nucleus & aNucleus) 255 //-------------------------------------------- 210 //---------------------------------------------------------------------------- 256 { 211 { 257 if(fBCDEBUG) G4cerr << " ######### Binary << 212 static G4int eventcounter=0; >> 213 >> 214 // if ( eventcounter == 0 ) { >> 215 // SetEpReportLevel(3); // report non conservation with model etc. >> 216 // G4double relativeLevel = 1*perCent; >> 217 // G4double absoluteLevel = 2*MeV; >> 218 // SetEnergyMomentumCheckLevels(relativeLevel,absoluteLevel); >> 219 // } >> 220 >> 221 //if(eventcounter == 100*(eventcounter/100) ) >> 222 eventcounter++; >> 223 if(getenv("BCDEBUG") ) G4cerr << " ######### Binary Cascade Reaction number starts ######### "<<eventcounter<<G4endl; 258 224 259 G4LorentzVector initial4Momentum = aTrack. 225 G4LorentzVector initial4Momentum = aTrack.Get4Momentum(); 260 const G4ParticleDefinition * definition = << 226 G4ParticleDefinition * definition = const_cast<G4ParticleDefinition *>(aTrack.GetDefinition()); 261 227 262 if(initial4Momentum.e()-initial4Momentum.m 228 if(initial4Momentum.e()-initial4Momentum.m()<theBCminP && 263 ( definition==G4Neutron::NeutronDe 229 ( definition==G4Neutron::NeutronDefinition() || definition==G4Proton::ProtonDefinition() ) ) 264 { 230 { 265 return theDeExcitation->ApplyYourself( 231 return theDeExcitation->ApplyYourself(aTrack, aNucleus); 266 } 232 } 267 233 268 theParticleChange.Clear(); 234 theParticleChange.Clear(); 269 // initialize the G4V3DNucleus from G4Nucl 235 // initialize the G4V3DNucleus from G4Nucleus 270 the3DNucleus = new G4Fancy3DNucleus; 236 the3DNucleus = new G4Fancy3DNucleus; 271 237 272 // Build a KineticTrackVector with the G4T 238 // Build a KineticTrackVector with the G4Track 273 G4KineticTrackVector * secondaries;// = ne 239 G4KineticTrackVector * secondaries;// = new G4KineticTrackVector; 274 G4ThreeVector initialPosition(0., 0., 0.); 240 G4ThreeVector initialPosition(0., 0., 0.); // will be set later 275 241 276 if(!fBCDEBUG) << 242 if(!getenv("I_Am_G4BinaryCascade_Developer") ) 277 { 243 { 278 if(definition!=G4Neutron::NeutronDefin 244 if(definition!=G4Neutron::NeutronDefinition() && 279 definition!=G4Proton::ProtonDe 245 definition!=G4Proton::ProtonDefinition() && 280 definition!=G4PionPlus::PionPl 246 definition!=G4PionPlus::PionPlusDefinition() && 281 definition!=G4PionMinus::PionM 247 definition!=G4PionMinus::PionMinusDefinition() ) 282 { 248 { 283 G4cerr << "You are trying to use G << 249 G4cerr << "You are using G4BinaryCascade for projectiles other than nucleons or pions."<<G4endl; 284 G4cerr << "G4BinaryCascade should << 285 G4cerr << "If you want to continue 250 G4cerr << "If you want to continue, please switch on the developer environment: "<<G4endl; 286 G4cerr << "setenv I_Am_G4BinaryCas 251 G4cerr << "setenv I_Am_G4BinaryCascade_Developer 1 "<<G4endl<<G4endl; 287 throw G4HadronicException(__FILE__ 252 throw G4HadronicException(__FILE__, __LINE__, "G4BinaryCascade - used for unvalid particle type - Fatal"); 288 } 253 } 289 } 254 } 290 255 291 // keep primary 256 // keep primary 292 thePrimaryType = definition; 257 thePrimaryType = definition; 293 thePrimaryEscape = false; 258 thePrimaryEscape = false; 294 259 295 G4double timePrimary=aTrack.GetGlobalTime( << 296 << 297 // try until an interaction will happen 260 // try until an interaction will happen 298 G4ReactionProductVector * products = nullp << 261 G4ReactionProductVector * products = 0; 299 G4int interactionCounter = 0,collisionLoop << 262 G4int interactionCounter = 0; 300 do 263 do 301 { 264 { 302 // reset status that could be changed 265 // reset status that could be changed in previous loop event 303 theCollisionMgr->ClearAndDestroy(); 266 theCollisionMgr->ClearAndDestroy(); 304 267 305 if(products != nullptr) << 268 if(products != 0) 306 { // free memory from previous loop e 269 { // free memory from previous loop event 307 ClearAndDestroy(products); 270 ClearAndDestroy(products); 308 delete products; 271 delete products; >> 272 products=0; 309 } 273 } 310 274 311 G4int massNumber=aNucleus.GetA_asInt() 275 G4int massNumber=aNucleus.GetA_asInt(); 312 the3DNucleus->Init(massNumber, aNucleu 276 the3DNucleus->Init(massNumber, aNucleus.GetZ_asInt()); 313 thePropagator->Init(the3DNucleus); 277 thePropagator->Init(the3DNucleus); 314 G4KineticTrack * kt; << 278 // GF Leak on kt??? but where to delete? 315 collisionLoopMaxCount = 200; << 279 G4KineticTrack * kt;// = new G4KineticTrack(definition, 0., initialPosition, initial4Momentum); 316 do // sample impact p 280 do // sample impact parameter until collisions are found 317 { 281 { 318 theCurrentTime=0; 282 theCurrentTime=0; 319 G4double radius = the3DNucleus->Ge 283 G4double radius = the3DNucleus->GetOuterRadius()+3*fermi; 320 initialPosition=GetSpherePoint(1.1 284 initialPosition=GetSpherePoint(1.1*radius, initial4Momentum); // get random position 321 kt = new G4KineticTrack(definition 285 kt = new G4KineticTrack(definition, 0., initialPosition, initial4Momentum); 322 kt->SetState(G4KineticTrack::outsi 286 kt->SetState(G4KineticTrack::outside); 323 // secondaries has been cleared by 287 // secondaries has been cleared by Propagate() in the previous loop event 324 secondaries= new G4KineticTrackVec 288 secondaries= new G4KineticTrackVector; 325 secondaries->push_back(kt); 289 secondaries->push_back(kt); 326 if(massNumber > 1) // 1H1 is speci 290 if(massNumber > 1) // 1H1 is special case 327 { 291 { 328 products = Propagate(secondari 292 products = Propagate(secondaries, the3DNucleus); 329 } else { 293 } else { 330 products = Propagate1H1(second 294 products = Propagate1H1(secondaries,the3DNucleus); 331 } 295 } 332 // until we FIND a collision ... o << 296 } while(! products ); // until we FIND a collision... 333 } while(! products && --collisionLoopM << 334 297 335 if(++interactionCounter>99) break; 298 if(++interactionCounter>99) break; 336 // ...until we find an ALLOWED collisi << 299 } while(products->size() == 0); // ...until we find an ALLOWED collision 337 } while(products && products->size() == 0) << 338 300 339 if(products && products->size()>0) << 301 if(products->size()>0) 340 { 302 { 341 // G4cout << "BIC Applyyourself: numb 303 // G4cout << "BIC Applyyourself: number of products " << products->size() << G4endl; 342 304 343 // Fill the G4ParticleChange * with pr 305 // Fill the G4ParticleChange * with products 344 theParticleChange.SetStatusChange(stop 306 theParticleChange.SetStatusChange(stopAndKill); 345 G4ReactionProductVector::iterator iter 307 G4ReactionProductVector::iterator iter; 346 308 347 for(iter = products->begin(); iter != 309 for(iter = products->begin(); iter != products->end(); ++iter) 348 { 310 { 349 G4DynamicParticle * aNewDP = << 311 G4DynamicParticle * aNew = 350 new G4DynamicParticle((*it 312 new G4DynamicParticle((*iter)->GetDefinition(), 351 (*iter)->GetTotalE 313 (*iter)->GetTotalEnergy(), 352 (*iter)->GetMoment 314 (*iter)->GetMomentum()); 353 G4HadSecondary aNew = G4HadSecondary << 354 G4double time=(*iter)->GetFormatio << 355 if(time < 0.0) { time = 0.0; } << 356 aNew.SetTime(timePrimary + time); << 357 aNew.SetCreatorModelID((*iter)->Ge << 358 aNew.SetParentResonanceDef((*iter) << 359 aNew.SetParentResonanceID((*iter)- << 360 theParticleChange.AddSecondary(aNe 315 theParticleChange.AddSecondary(aNew); 361 } 316 } 362 317 363 //DebugFinalEpConservation(aTrack, pr << 318 // DebugEpConservation(track, products); 364 319 365 320 366 } else { // no interaction, return primar 321 } else { // no interaction, return primary 367 if(fBCDEBUG) G4cerr << " ######### Bin << 322 if(getenv("BCDEBUG") ) G4cerr << " ######### Binary Cascade Reaction number void ######### "<<eventcounter<<G4endl; 368 theParticleChange.SetStatusChange(isAl 323 theParticleChange.SetStatusChange(isAlive); 369 theParticleChange.SetEnergyChange(aTra 324 theParticleChange.SetEnergyChange(aTrack.GetKineticEnergy()); 370 theParticleChange.SetMomentumChange(aT 325 theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit()); 371 } 326 } 372 327 373 if ( products ) << 328 ClearAndDestroy(products); 374 { << 329 delete products; 375 ClearAndDestroy(products); << 330 376 delete products; << 377 } << 378 << 379 delete the3DNucleus; 331 delete the3DNucleus; 380 the3DNucleus = nullptr; << 332 the3DNucleus = NULL; 381 333 382 if(fBCDEBUG) G4cerr << " ######### Binary << 334 if(getenv("BCDEBUG") ) G4cerr << " ######### Binary Cascade Reaction number ends ######### "<<eventcounter<<G4endl; 383 335 384 return &theParticleChange; 336 return &theParticleChange; 385 } 337 } 386 //-------------------------------------------- 338 //---------------------------------------------------------------------------- 387 G4ReactionProductVector * G4BinaryCascade::Pro 339 G4ReactionProductVector * G4BinaryCascade::Propagate( 388 G4KineticTrackVector * secondaries, G4 340 G4KineticTrackVector * secondaries, G4V3DNucleus * aNucleus) 389 //-------------------------------------------- 341 //---------------------------------------------------------------------------- 390 { 342 { 391 G4ping debug("debug_G4BinaryCascade"); 343 G4ping debug("debug_G4BinaryCascade"); 392 #ifdef debug_BIC_Propagate 344 #ifdef debug_BIC_Propagate 393 G4cout << "G4BinaryCascade Propagate start 345 G4cout << "G4BinaryCascade Propagate starting -------------------------------------------------------" <<G4endl; 394 #endif 346 #endif 395 347 396 the3DNucleus=aNucleus; 348 the3DNucleus=aNucleus; 397 G4ReactionProductVector * products = new G 349 G4ReactionProductVector * products = new G4ReactionProductVector; 398 theOuterRadius = the3DNucleus->GetOuterRad 350 theOuterRadius = the3DNucleus->GetOuterRadius(); 399 theCurrentTime=0; 351 theCurrentTime=0; 400 theProjectile4Momentum=G4LorentzVector(0,0 352 theProjectile4Momentum=G4LorentzVector(0,0,0,0); 401 theMomentumTransfer=G4ThreeVector(0,0,0); << 402 // build theSecondaryList, theProjectileLi 353 // build theSecondaryList, theProjectileList and theCapturedList 403 ClearAndDestroy(&theCapturedList); 354 ClearAndDestroy(&theCapturedList); 404 ClearAndDestroy(&theSecondaryList); 355 ClearAndDestroy(&theSecondaryList); 405 theSecondaryList.clear(); 356 theSecondaryList.clear(); 406 ClearAndDestroy(&theFinalState); 357 ClearAndDestroy(&theFinalState); 407 std::vector<G4KineticTrack *>::iterator it 358 std::vector<G4KineticTrack *>::iterator iter; 408 theCollisionMgr->ClearAndDestroy(); 359 theCollisionMgr->ClearAndDestroy(); 409 360 410 theCutOnP=90*MeV; 361 theCutOnP=90*MeV; 411 if(the3DNucleus->GetMass()>30) theCutOnP = 362 if(the3DNucleus->GetMass()>30) theCutOnP = 70*MeV; 412 if(the3DNucleus->GetMass()>60) theCutOnP = 363 if(the3DNucleus->GetMass()>60) theCutOnP = 50*MeV; 413 if(the3DNucleus->GetMass()>120) theCutOnP 364 if(the3DNucleus->GetMass()>120) theCutOnP = 45*MeV; 414 365 415 366 416 BuildTargetList(); 367 BuildTargetList(); 417 368 418 #ifdef debug_BIC_GetExcitationEnergy 369 #ifdef debug_BIC_GetExcitationEnergy 419 G4cout << "ExcitationEnergy0 " << GetExcit 370 G4cout << "ExcitationEnergy0 " << GetExcitationEnergy() << G4endl; 420 #endif 371 #endif 421 372 422 thePropagator->Init(the3DNucleus); 373 thePropagator->Init(the3DNucleus); 423 374 424 G4bool success = BuildLateParticleCollisio 375 G4bool success = BuildLateParticleCollisions(secondaries); 425 if (! success ) // fails if no excitatio 376 if (! success ) // fails if no excitation energy left.... 426 { 377 { 427 products=HighEnergyModelFSProducts(prod 378 products=HighEnergyModelFSProducts(products, secondaries); 428 ClearAndDestroy(secondaries); 379 ClearAndDestroy(secondaries); 429 delete secondaries; 380 delete secondaries; 430 381 431 #ifdef debug_G4BinaryCascade 382 #ifdef debug_G4BinaryCascade 432 G4cout << "G4BinaryCascade::Propagate: 383 G4cout << "G4BinaryCascade::Propagate: warning - high energy model failed energy conservation, returning unchanged high energy final state" << G4endl; 433 #endif 384 #endif 434 385 435 return products; 386 return products; 436 } 387 } 437 // check baryon and charge ... 388 // check baryon and charge ... 438 389 439 _CheckChargeAndBaryonNumber_("lateparticle 390 _CheckChargeAndBaryonNumber_("lateparticles"); 440 _DebugEpConservation(" be4 findcollisions" << 441 391 442 // if called stand alone find first collis 392 // if called stand alone find first collisions 443 FindCollisions(&theSecondaryList); 393 FindCollisions(&theSecondaryList); 444 394 445 395 446 if(theCollisionMgr->Entries() == 0 ) 396 if(theCollisionMgr->Entries() == 0 ) //late particles ALWAYS create Entries 447 { 397 { 448 //G4cout << " no collsions -> return 0 398 //G4cout << " no collsions -> return 0" << G4endl; 449 delete products; 399 delete products; 450 #ifdef debug_BIC_return 400 #ifdef debug_BIC_return 451 G4cout << "return @ begin2, no collis 401 G4cout << "return @ begin2, no collisions "<< G4endl; 452 #endif 402 #endif 453 return nullptr; << 403 return 0; 454 } 404 } 455 405 456 // end of initialization: do the job now 406 // end of initialization: do the job now 457 // loop until there are no more collisions 407 // loop until there are no more collisions 458 408 459 409 460 G4bool haveProducts = false; 410 G4bool haveProducts = false; 461 #ifdef debug_BIC_Propagate_Collisions << 462 G4int collisionCount=0; 411 G4int collisionCount=0; 463 #endif << 412 theMomentumTransfer=G4ThreeVector(0,0,0); 464 G4int collisionLoopMaxCount=1000000; << 413 while(theCollisionMgr->Entries() > 0 && currentZ) 465 while(theCollisionMgr->Entries() > 0 && cu << 466 { 414 { 467 if(Absorb()) { // absorb secondaries, p << 415 if(Absorb()) { // absorb secondaries, pions only 468 haveProducts = true; 416 haveProducts = true; 469 } << 417 } 470 if(Capture()) { // capture secondaries, << 418 _CheckChargeAndBaryonNumber_("absorbed"); >> 419 if(Capture()) { // capture secondaries, nucleons only 471 haveProducts = true; 420 haveProducts = true; 472 } 421 } >> 422 _CheckChargeAndBaryonNumber_("captured"); 473 423 474 // propagate to the next collision if 424 // propagate to the next collision if any (collisions could have been deleted 475 // by previous absorption or capture) 425 // by previous absorption or capture) 476 if(theCollisionMgr->Entries() > 0) 426 if(theCollisionMgr->Entries() > 0) 477 { 427 { 478 G4CollisionInitialState * 428 G4CollisionInitialState * 479 nextCollision = theCollisionMgr->G 429 nextCollision = theCollisionMgr->GetNextCollision(); 480 #ifdef debug_BIC_Propagate_Collisions 430 #ifdef debug_BIC_Propagate_Collisions 481 G4cout << " NextCollision * , Tim 431 G4cout << " NextCollision * , Time, curtime = " << nextCollision << " " 482 <<nextCollision->GetCollis 432 <<nextCollision->GetCollisionTime()<< " " << 483 theCurrentTime<< G4endl; 433 theCurrentTime<< G4endl; 484 #endif 434 #endif 485 if (!DoTimeStep(nextCollision->Get 435 if (!DoTimeStep(nextCollision->GetCollisionTime()-theCurrentTime) ) 486 { 436 { 487 // Check if nextCollision is s 437 // Check if nextCollision is still valid, ie. particle did not leave nucleus 488 if (theCollisionMgr->GetNextCo 438 if (theCollisionMgr->GetNextCollision() != nextCollision ) 489 { 439 { 490 nextCollision = nullptr; << 440 nextCollision = 0; 491 } 441 } 492 } 442 } 493 //_DebugEpConservation("Stepped"); << 494 443 495 if( nextCollision ) 444 if( nextCollision ) 496 { 445 { 497 if (ApplyCollision(nextCollisi 446 if (ApplyCollision(nextCollision)) 498 { 447 { 499 //G4cerr << "ApplyCollisio 448 //G4cerr << "ApplyCollision success " << G4endl; 500 haveProducts = true; 449 haveProducts = true; 501 #ifdef debug_BIC_Propagate_Collisions << 502 collisionCount++; 450 collisionCount++; 503 #endif << 451 _CheckChargeAndBaryonNumber_("ApplyCollision"); 504 452 505 } else { 453 } else { 506 //G4cerr << "ApplyCollisio 454 //G4cerr << "ApplyCollision failure " << G4endl; 507 theCollisionMgr->RemoveCol 455 theCollisionMgr->RemoveCollision(nextCollision); 508 } 456 } 509 } 457 } 510 } 458 } 511 } 459 } 512 460 513 //--------- end of on Collisions << 461 //--------- end of while on Collisions 514 //G4cout << "currentZ @ end loop " << curr 462 //G4cout << "currentZ @ end loop " << currentZ << G4endl; 515 G4int nProtons(0); << 463 if ( ! currentZ ){ 516 for(iter = theTargetList.begin(); iter != << 517 { << 518 if ( (*iter)->GetDefinition() == G4Proto << 519 } << 520 if ( ! theTargetList.size() || ! nProtons << 521 // nucleus completely destroyed, fill 464 // nucleus completely destroyed, fill in ReactionProductVector 522 products = FillVoidNucleusProducts(prod << 465 products = FillVoidNucleusProducts(products); 523 #ifdef debug_BIC_return 466 #ifdef debug_BIC_return 524 G4cout << "return @ Z=0 after collisio 467 G4cout << "return @ Z=0 after collision loop "<< G4endl; 525 PrintKTVector(&theSecondaryList,std::s 468 PrintKTVector(&theSecondaryList,std::string(" theSecondaryList")); 526 G4cout << "theTargetList size: " << th 469 G4cout << "theTargetList size: " << theTargetList.size() << G4endl; 527 PrintKTVector(&theTargetList,std::stri 470 PrintKTVector(&theTargetList,std::string(" theTargetList")); 528 PrintKTVector(&theCapturedList,std::st 471 PrintKTVector(&theCapturedList,std::string(" theCapturedList")); 529 472 530 G4cout << " ExcitE be4 Correct : " <<G 473 G4cout << " ExcitE be4 Correct : " <<GetExcitationEnergy() << G4endl; 531 G4cout << " Mom Transfered to nucleus 474 G4cout << " Mom Transfered to nucleus : " << theMomentumTransfer << " " << theMomentumTransfer.mag() << G4endl; 532 PrintKTVector(&theFinalState,std::stri 475 PrintKTVector(&theFinalState,std::string(" FinalState uncorrected")); 533 G4cout << "returned products: " << pro 476 G4cout << "returned products: " << products->size() << G4endl; 534 _CheckChargeAndBaryonNumber_("destroye << 535 _DebugEpConservation("destroyed Nucleu << 536 #endif 477 #endif 537 << 478 // @@GF FixMe cannot return here. 538 return products; 479 return products; 539 } 480 } 540 481 541 // No more collisions: absorb, capture and 482 // No more collisions: absorb, capture and propagate the secondaries out of the nucleus 542 if(Absorb()) { 483 if(Absorb()) { 543 haveProducts = true; 484 haveProducts = true; 544 // G4cout << "Absorb sucess " << G4end 485 // G4cout << "Absorb sucess " << G4endl; 545 } 486 } 546 487 547 if(Capture()) { 488 if(Capture()) { 548 haveProducts = true; 489 haveProducts = true; 549 // G4cout << "Capture sucess " << G4en 490 // G4cout << "Capture sucess " << G4endl; 550 } 491 } 551 492 552 if(!haveProducts) // no collisions happen 493 if(!haveProducts) // no collisions happened. Return an empty vector. 553 { 494 { 554 #ifdef debug_BIC_return 495 #ifdef debug_BIC_return 555 G4cout << "return 3, no products "<< G 496 G4cout << "return 3, no products "<< G4endl; 556 #endif 497 #endif 557 return products; 498 return products; 558 } 499 } 559 500 560 501 561 #ifdef debug_BIC_Propagate 502 #ifdef debug_BIC_Propagate 562 G4cout << " Momentum transfer to Nucleus " 503 G4cout << " Momentum transfer to Nucleus " << theMomentumTransfer << " " << theMomentumTransfer.mag() << G4endl; 563 G4cout << " Stepping particles out...... 504 G4cout << " Stepping particles out...... " << G4endl; 564 #endif 505 #endif 565 506 566 StepParticlesOut(); 507 StepParticlesOut(); 567 _DebugEpConservation("stepped out"); << 568 508 569 509 570 if ( theSecondaryList.size() > 0 ) 510 if ( theSecondaryList.size() > 0 ) 571 { 511 { 572 #ifdef debug_G4BinaryCascade 512 #ifdef debug_G4BinaryCascade 573 G4cerr << "G4BinaryCascade: Warning, h 513 G4cerr << "G4BinaryCascade: Warning, have active particles at end" << G4endl; 574 PrintKTVector(&theSecondaryList, "acti 514 PrintKTVector(&theSecondaryList, "active particles @ end added to theFinalState"); 575 #endif 515 #endif 576 // add left secondaries to FinalSate 516 // add left secondaries to FinalSate 577 for ( iter =theSecondaryList.begin(); 517 for ( iter =theSecondaryList.begin(); iter != theSecondaryList.end(); ++iter) 578 { 518 { 579 theFinalState.push_back(*iter); 519 theFinalState.push_back(*iter); 580 } 520 } 581 theSecondaryList.clear(); 521 theSecondaryList.clear(); 582 522 583 } 523 } 584 while ( theCollisionMgr->Entries() > 0 ) << 524 while ( theCollisionMgr->Entries() > 0 ) 585 { 525 { 586 #ifdef debug_G4BinaryCascade 526 #ifdef debug_G4BinaryCascade 587 G4cerr << " Warning: remove left over 527 G4cerr << " Warning: remove left over collision(s) " << G4endl; 588 #endif 528 #endif 589 theCollisionMgr->RemoveCollision(theCo 529 theCollisionMgr->RemoveCollision(theCollisionMgr->GetNextCollision()); 590 } 530 } 591 531 592 #ifdef debug_BIC_Propagate_Excitation 532 #ifdef debug_BIC_Propagate_Excitation 593 533 594 PrintKTVector(&theSecondaryList,std::strin 534 PrintKTVector(&theSecondaryList,std::string(" theSecondaryList")); 595 G4cout << "theTargetList size: " << theTar 535 G4cout << "theTargetList size: " << theTargetList.size() << G4endl; 596 // PrintKTVector(&theTargetList,std::stri 536 // PrintKTVector(&theTargetList,std::string(" theTargetList")); 597 PrintKTVector(&theCapturedList,std::string 537 PrintKTVector(&theCapturedList,std::string(" theCapturedList")); 598 538 599 G4cout << " ExcitE be4 Correct : " <<GetEx 539 G4cout << " ExcitE be4 Correct : " <<GetExcitationEnergy() << G4endl; 600 G4cout << " Mom Transfered to nucleus : " 540 G4cout << " Mom Transfered to nucleus : " << theMomentumTransfer << " " << theMomentumTransfer.mag() << G4endl; 601 PrintKTVector(&theFinalState,std::string(" 541 PrintKTVector(&theFinalState,std::string(" FinalState uncorrected")); 602 #endif 542 #endif 603 543 604 // 544 // 605 545 606 546 607 G4double ExcitationEnergy=GetExcitationEne 547 G4double ExcitationEnergy=GetExcitationEnergy(); 608 548 609 #ifdef debug_BIC_Propagate_finals 549 #ifdef debug_BIC_Propagate_finals 610 PrintKTVector(&theFinalState,std::string(" 550 PrintKTVector(&theFinalState,std::string(" FinalState be4 corr")); 611 G4cout << " Excitation Energy prefinal, # 551 G4cout << " Excitation Energy prefinal, #collisions:, out, captured " 612 << ExcitationEnergy << " " 552 << ExcitationEnergy << " " 613 << collisionCount << " " 553 << collisionCount << " " 614 << theFinalState.size() << " " 554 << theFinalState.size() << " " 615 << theCapturedList.size()<<G4endl; 555 << theCapturedList.size()<<G4endl; 616 #endif 556 #endif 617 557 618 if (ExcitationEnergy < 0 ) 558 if (ExcitationEnergy < 0 ) 619 { 559 { 620 G4int maxtry=5, ntry=0; 560 G4int maxtry=5, ntry=0; 621 do { 561 do { 622 CorrectFinalPandE(); 562 CorrectFinalPandE(); 623 ExcitationEnergy=GetExcitationEner 563 ExcitationEnergy=GetExcitationEnergy(); 624 } while ( ++ntry < maxtry && Excitatio << 564 } while ( ++ntry < maxtry && ExcitationEnergy < 0 ); 625 } 565 } 626 _DebugEpConservation("corrected"); << 627 566 628 #ifdef debug_BIC_Propagate_finals 567 #ifdef debug_BIC_Propagate_finals 629 PrintKTVector(&theFinalState,std::string(" 568 PrintKTVector(&theFinalState,std::string(" FinalState corrected")); 630 G4cout << " Excitation Energy final, #col 569 G4cout << " Excitation Energy final, #collisions:, out, captured " 631 << ExcitationEnergy << " " 570 << ExcitationEnergy << " " 632 << collisionCount << " " 571 << collisionCount << " " 633 << theFinalState.size() << " " 572 << theFinalState.size() << " " 634 << theCapturedList.size()<<G4endl; 573 << theCapturedList.size()<<G4endl; 635 #endif 574 #endif 636 575 637 576 638 if ( ExcitationEnergy < 0. ) 577 if ( ExcitationEnergy < 0. ) 639 { 578 { 640 #ifdef debug_G4BinaryCascade << 579 // if ( ExcitationEnergy < 0. ) 641 G4cerr << "G4BinaryCascade-W << 580 { 642 G4cerr <<ExcitationEnergy<<G << 581 //#ifdef debug_G4BinaryCascade 643 PrintKTVector(&theFinalState, << 582 // G4cerr << "G4BinaryCascade-Warning: negative excitation energy "; 644 PrintKTVector(&theCapturedLis << 583 // G4cerr <<ExcitationEnergy<<G4endl; 645 G4cout << "negative ExE:Final << 584 // PrintKTVector(&theFinalState,std::string("FinalState")); 646 << " "<< GetFinal4Mome << 585 // PrintKTVector(&theCapturedList,std::string("captured")); 647 << "negative ExE:Final << 586 // G4cout << "negative ExE:Final 4Momentum .mag: " << GetFinal4Momentum() 648 << " "<< GetFinalNucleusMome << 587 // << " "<< GetFinal4Momentum().mag()<< G4endl 649 #endif << 588 // << "negative ExE:FinalNucleusMom .mag: " << GetFinalNucleusMomentum() 650 #ifdef debug_BIC_return << 589 // << " "<< GetFinalNucleusMomentum().mag()<< G4endl; 651 G4cout << " negative Exci << 590 //#endif 652 G4cout << "return 4, excit << 591 } 653 #endif << 654 << 655 ClearAndDestroy(products); 592 ClearAndDestroy(products); 656 return products; // return empty pro << 593 //G4cout << " negative Excitation E return empty products " << products << G4endl; >> 594 #ifdef debug_BIC_return >> 595 G4cout << "return 4, excit < 0 "<< G4endl; >> 596 #endif >> 597 return products; // return empty products- FIXME 657 } 598 } 658 599 659 G4ReactionProductVector * precompoundProdu 600 G4ReactionProductVector * precompoundProducts=DeExcite(); 660 601 661 602 662 G4DecayKineticTracks decay(&theFinalState) 603 G4DecayKineticTracks decay(&theFinalState); 663 _DebugEpConservation("decayed"); << 664 604 665 products= ProductsAddFinalState(products, 605 products= ProductsAddFinalState(products, theFinalState); 666 606 667 products= ProductsAddPrecompound(products, 607 products= ProductsAddPrecompound(products, precompoundProducts); 668 608 669 // products=ProductsAddFakeGamma(products); 609 // products=ProductsAddFakeGamma(products); 670 610 671 611 672 thePrimaryEscape = true; 612 thePrimaryEscape = true; 673 613 674 #ifdef debug_BIC_return 614 #ifdef debug_BIC_return 675 G4cout << "BIC: return @end, all ok "<< G4 << 615 G4cout << "return @end, all ok "<< G4endl; 676 //G4cout << " return products " << produc 616 //G4cout << " return products " << products << G4endl; 677 #endif 617 #endif 678 618 679 return products; 619 return products; 680 } 620 } 681 621 682 //-------------------------------------------- 622 //---------------------------------------------------------------------------- 683 G4double G4BinaryCascade::GetExcitationEnergy( 623 G4double G4BinaryCascade::GetExcitationEnergy() 684 //-------------------------------------------- 624 //---------------------------------------------------------------------------- 685 { 625 { 686 626 687 // get A and Z for the residual nucleus 627 // get A and Z for the residual nucleus 688 #if defined(debug_G4BinaryCascade) || defined( 628 #if defined(debug_G4BinaryCascade) || defined(debug_BIC_GetExcitationEnergy) 689 G4int finalA = theTargetList.size()+theCap 629 G4int finalA = theTargetList.size()+theCapturedList.size(); 690 G4int finalZ = GetTotalCharge(theTargetLis 630 G4int finalZ = GetTotalCharge(theTargetList)+GetTotalCharge(theCapturedList); 691 if ( (currentA - finalA) != 0 || (currentZ 631 if ( (currentA - finalA) != 0 || (currentZ - finalZ) != 0 ) 692 { 632 { 693 G4cerr << "G4BIC:GetExcitationEnergy() 633 G4cerr << "G4BIC:GetExcitationEnergy(): Nucleon counting error current/final{A,Z} " 694 << "("<< currentA << "," << fi << 634 << currentA << " " << finalA << " "<< currentZ << " " << finalZ << G4endl; 695 } 635 } 696 636 697 #endif 637 #endif 698 638 699 G4double excitationE(0); 639 G4double excitationE(0); 700 G4double nucleusMass(0); 640 G4double nucleusMass(0); 701 if(currentZ>.5) 641 if(currentZ>.5) 702 { 642 { 703 nucleusMass = GetIonMass(currentZ,curr << 643 nucleusMass = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(currentZ,currentA); 704 } << 705 else if (currentZ==0 ) << 706 { << 707 if(currentA == 1) {nucleusMass = G4Neu << 708 else {nucleusMass = GetFi << 709 } 644 } >> 645 else if (currentZ==0 ) // Uzhi && currentA==1 ) // Uzhi >> 646 { // Uzhi >> 647 if(currentA == 1) {nucleusMass = G4Neutron::Neutron()->GetPDGMass();}// Uzhi >> 648 else {nucleusMass = GetFinalNucleusMomentum().mag() // Uzhi >> 649 - 3.*MeV*currentA;} // Uzhi >> 650 } // Uzhi 710 else 651 else 711 { 652 { 712 #ifdef debug_G4BinaryCascade 653 #ifdef debug_G4BinaryCascade 713 G4cout << "G4BinaryCascade::GetExcitat 654 G4cout << "G4BinaryCascade::GetExcitationEnergy(): Warning - invalid nucleus (A,Z)=(" 714 << currentA << "," << currentZ 655 << currentA << "," << currentZ << ")" << G4endl; 715 #endif 656 #endif 716 return 0; 657 return 0; 717 } 658 } 718 659 719 #ifdef debug_BIC_GetExcitationEnergy 660 #ifdef debug_BIC_GetExcitationEnergy 720 G4ping debug("debug_ExcitationEnergy"); 661 G4ping debug("debug_ExcitationEnergy"); 721 debug.push_back("====> current A, Z"); 662 debug.push_back("====> current A, Z"); 722 debug.push_back(currentZ); 663 debug.push_back(currentZ); 723 debug.push_back(currentA); 664 debug.push_back(currentA); 724 debug.push_back("====> final A, Z"); 665 debug.push_back("====> final A, Z"); 725 debug.push_back(finalZ); 666 debug.push_back(finalZ); 726 debug.push_back(finalA); 667 debug.push_back(finalA); 727 debug.push_back(nucleusMass); 668 debug.push_back(nucleusMass); 728 debug.push_back(GetFinalNucleusMomentum(). 669 debug.push_back(GetFinalNucleusMomentum().mag()); 729 debug.dump(); 670 debug.dump(); 730 // PrintKTVector(&theTargetList, std::str 671 // PrintKTVector(&theTargetList, std::string(" current target list info")); 731 //PrintKTVector(&theCapturedList, std::str 672 //PrintKTVector(&theCapturedList, std::string(" current captured list info")); 732 #endif 673 #endif 733 674 734 excitationE = GetFinalNucleusMomentum().ma 675 excitationE = GetFinalNucleusMomentum().mag() - nucleusMass; 735 676 736 //G4double exE2 = GetFinal4Momentum().mag( 677 //G4double exE2 = GetFinal4Momentum().mag() - nucleusMass; 737 678 738 //G4cout << "old/new excitE " << excitatio 679 //G4cout << "old/new excitE " << excitationE << " / "<< exE2 << G4endl; 739 680 740 #ifdef debug_BIC_GetExcitationEnergy 681 #ifdef debug_BIC_GetExcitationEnergy 741 // ------ debug 682 // ------ debug 742 if ( excitationE < 0 ) 683 if ( excitationE < 0 ) 743 { 684 { 744 G4cout << "negative ExE final Ion mass 685 G4cout << "negative ExE final Ion mass " <<nucleusMass<< G4endl; 745 G4LorentzVector Nucl_mom=GetFinalNucle 686 G4LorentzVector Nucl_mom=GetFinalNucleusMomentum(); 746 if(finalZ>.5) G4cout << " Final nuclmo 687 if(finalZ>.5) G4cout << " Final nuclmom/mass " << Nucl_mom << " " << Nucl_mom.mag() 747 << " (A,Z)=("<< finalA 688 << " (A,Z)=("<< finalA <<","<<finalZ <<")" 748 << " mass " << nucleusM 689 << " mass " << nucleusMass << " " 749 << " excitE " << excita 690 << " excitE " << excitationE << G4endl; 750 691 751 692 752 G4int A = the3DNucleus->GetMassNumber( 693 G4int A = the3DNucleus->GetMassNumber(); 753 G4int Z = the3DNucleus->GetCharge(); 694 G4int Z = the3DNucleus->GetCharge(); 754 G4double initialExc(0); 695 G4double initialExc(0); 755 if(Z>.5) 696 if(Z>.5) 756 { 697 { 757 initialExc = theInitial4Mom.mag()- << 698 initialExc = theInitial4Mom.mag()- >> 699 G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(Z, A); 758 G4cout << "GetExcitationEnergy: In 700 G4cout << "GetExcitationEnergy: Initial nucleus A Z " << A << " " << Z << " " << initialExc << G4endl; 759 } 701 } 760 } 702 } 761 703 762 #endif 704 #endif 763 705 764 return excitationE; 706 return excitationE; 765 } 707 } 766 708 767 709 768 //-------------------------------------------- 710 //---------------------------------------------------------------------------- 769 // 711 // 770 // P R I V A T E M E M B E R F U N C 712 // P R I V A T E M E M B E R F U N C T I O N S 771 // 713 // 772 //-------------------------------------------- 714 //---------------------------------------------------------------------------- 773 715 774 //-------------------------------------------- 716 //---------------------------------------------------------------------------- 775 void G4BinaryCascade::BuildTargetList() 717 void G4BinaryCascade::BuildTargetList() 776 //-------------------------------------------- 718 //---------------------------------------------------------------------------- 777 { 719 { 778 720 779 if(!the3DNucleus->StartLoop()) 721 if(!the3DNucleus->StartLoop()) 780 { 722 { 781 // G4cerr << "G4BinaryCascade::Buil 723 // G4cerr << "G4BinaryCascade::BuildTargetList(): StartLoop() error!" 782 // << G4endl; 724 // << G4endl; 783 return; 725 return; 784 } 726 } 785 727 786 ClearAndDestroy(&theTargetList); // clear 728 ClearAndDestroy(&theTargetList); // clear theTargetList before rebuilding 787 729 788 G4Nucleon * nucleon; 730 G4Nucleon * nucleon; 789 const G4ParticleDefinition * definition; << 731 G4ParticleDefinition * definition; 790 G4ThreeVector pos; 732 G4ThreeVector pos; 791 G4LorentzVector mom; 733 G4LorentzVector mom; 792 // if there are nucleon hit by higher ener 734 // if there are nucleon hit by higher energy models, then SUM(momenta) != 0 793 initialZ=the3DNucleus->GetCharge(); 735 initialZ=the3DNucleus->GetCharge(); 794 initialA=the3DNucleus->GetMassNumber(); 736 initialA=the3DNucleus->GetMassNumber(); 795 initial_nuclear_mass=GetIonMass(initialZ,i << 737 initial_nuclear_mass=G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(initialZ,initialA); 796 theInitial4Mom = G4LorentzVector(0,0,0,ini 738 theInitial4Mom = G4LorentzVector(0,0,0,initial_nuclear_mass); 797 currentA=0; 739 currentA=0; 798 currentZ=0; 740 currentZ=0; 799 while((nucleon = the3DNucleus->GetNextNucl << 741 while((nucleon = the3DNucleus->GetNextNucleon()) != NULL) 800 { 742 { 801 // check if nucleon is hit by higher e 743 // check if nucleon is hit by higher energy model. 802 if ( ! nucleon->AreYouHit() ) 744 if ( ! nucleon->AreYouHit() ) 803 { 745 { 804 definition = nucleon->GetDefinitio 746 definition = nucleon->GetDefinition(); 805 pos = nucleon->GetPosition(); 747 pos = nucleon->GetPosition(); 806 mom = nucleon->GetMomentum(); 748 mom = nucleon->GetMomentum(); 807 // G4cout << "Nucleus " << pos. 749 // G4cout << "Nucleus " << pos.mag()/fermi << " " << mom.e() << G4endl; 808 //theInitial4Mom += mom; 750 //theInitial4Mom += mom; 809 // the potential inside the 751 // the potential inside the nucleus is taken into account, and nucleons are on mass shell. 810 mom.setE( std::sqrt( mom.vect().ma 752 mom.setE( std::sqrt( mom.vect().mag2() + sqr(definition->GetPDGMass()) ) ); 811 G4KineticTrack * kt = new G4Kineti 753 G4KineticTrack * kt = new G4KineticTrack(definition, 0., pos, mom); 812 kt->SetState(G4KineticTrack::insid 754 kt->SetState(G4KineticTrack::inside); 813 kt->SetNucleon(nucleon); 755 kt->SetNucleon(nucleon); 814 theTargetList.push_back(kt); 756 theTargetList.push_back(kt); 815 ++currentA; 757 ++currentA; 816 if (definition->GetPDGCharge() > . 758 if (definition->GetPDGCharge() > .5 ) ++currentZ; 817 } 759 } 818 760 819 #ifdef debug_BIC_BuildTargetList 761 #ifdef debug_BIC_BuildTargetList 820 else { G4cout << "nucleon is hit" << n 762 else { G4cout << "nucleon is hit" << nucleon << G4endl;} 821 #endif 763 #endif 822 } 764 } 823 massInNucleus = 0; 765 massInNucleus = 0; 824 if(currentZ>.5) 766 if(currentZ>.5) 825 { 767 { 826 massInNucleus = GetIonMass(currentZ,cu << 768 massInNucleus = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(currentZ,currentA); 827 } else if (currentZ==0 && currentA>=1 ) 769 } else if (currentZ==0 && currentA>=1 ) 828 { 770 { 829 massInNucleus = currentA * G4Neutron:: 771 massInNucleus = currentA * G4Neutron::Neutron()->GetPDGMass(); 830 } else 772 } else 831 { 773 { 832 G4cerr << "G4BinaryCascade::BuildTarge 774 G4cerr << "G4BinaryCascade::BuildTargetList(): Fatal Error - invalid nucleus (A,Z)=(" 833 << currentA << "," << currentZ 775 << currentA << "," << currentZ << ")" << G4endl; 834 throw G4HadronicException(__FILE__, __ 776 throw G4HadronicException(__FILE__, __LINE__, "G4BinaryCasacde::BuildTargetList()"); 835 } 777 } 836 currentInitialEnergy= theInitial4Mom.e() + 778 currentInitialEnergy= theInitial4Mom.e() + theProjectile4Momentum.e(); 837 779 838 #ifdef debug_BIC_BuildTargetList 780 #ifdef debug_BIC_BuildTargetList 839 G4cout << "G4BinaryCascade::BuildTargetLis 781 G4cout << "G4BinaryCascade::BuildTargetList(): nucleus (A,Z)=(" 840 << currentA << "," << currentZ << 782 << currentA << "," << currentZ << ") mass: " << massInNucleus << 841 ", theInitial4Mom " << theInitial4 783 ", theInitial4Mom " << theInitial4Mom << 842 ", currentInitialEnergy " << curre 784 ", currentInitialEnergy " << currentInitialEnergy << G4endl; 843 #endif 785 #endif 844 786 845 } 787 } 846 788 847 //-------------------------------------------- 789 //---------------------------------------------------------------------------- 848 G4bool G4BinaryCascade::BuildLateParticleColl 790 G4bool G4BinaryCascade::BuildLateParticleCollisions(G4KineticTrackVector * secondaries) 849 //-------------------------------------------- 791 //---------------------------------------------------------------------------- 850 { 792 { 851 G4bool success(false); 793 G4bool success(false); 852 std::vector<G4KineticTrack *>::iterator ite 794 std::vector<G4KineticTrack *>::iterator iter; 853 795 854 lateA=lateZ=0; 796 lateA=lateZ=0; 855 projectileA=projectileZ=0; 797 projectileA=projectileZ=0; 856 798 857 G4double StartingTime=DBL_MAX; // Se 799 G4double StartingTime=DBL_MAX; // Search for minimal formation time 858 for(iter = secondaries->begin(); iter != se 800 for(iter = secondaries->begin(); iter != secondaries->end(); ++iter) 859 { 801 { 860 if((*iter)->GetFormationTime() < Startin 802 if((*iter)->GetFormationTime() < StartingTime) 861 StartingTime = (*iter)->GetFormationT 803 StartingTime = (*iter)->GetFormationTime(); 862 } 804 } 863 805 864 //PrintKTVector(secondaries, "initial late 806 //PrintKTVector(secondaries, "initial late particles "); 865 G4LorentzVector lateParticles4Momentum(0,0, 807 G4LorentzVector lateParticles4Momentum(0,0,0,0); 866 for(iter = secondaries->begin(); iter != se 808 for(iter = secondaries->begin(); iter != secondaries->end(); ++iter) 867 { 809 { 868 // G4cout << " Formation time : " << (* 810 // G4cout << " Formation time : " << (*iter)->GetDefinition()->GetParticleName() << " " 869 // << (*iter)->GetFormationTime() << G 811 // << (*iter)->GetFormationTime() << G4endl; 870 G4double FormTime = (*iter)->GetFormatio 812 G4double FormTime = (*iter)->GetFormationTime() - StartingTime; 871 (*iter)->SetFormationTime(FormTime); 813 (*iter)->SetFormationTime(FormTime); 872 if( (*iter)->GetState() == G4KineticTrac 814 if( (*iter)->GetState() == G4KineticTrack::undefined ) // particles from high energy generator 873 { 815 { 874 FindLateParticleCollision(*iter); 816 FindLateParticleCollision(*iter); 875 lateParticles4Momentum += (*iter)->Ge 817 lateParticles4Momentum += (*iter)->GetTrackingMomentum(); 876 lateA += (*iter)->GetDefinition()->Ge 818 lateA += (*iter)->GetDefinition()->GetBaryonNumber(); 877 lateZ += G4lrint((*iter)->GetDefiniti 819 lateZ += G4lrint((*iter)->GetDefinition()->GetPDGCharge()/eplus); 878 //PrintKTVector(*iter, "late particle 820 //PrintKTVector(*iter, "late particle "); 879 } else 821 } else 880 { 822 { 881 theSecondaryList.push_back(*iter); 823 theSecondaryList.push_back(*iter); 882 //PrintKTVector(*iter, "incoming part 824 //PrintKTVector(*iter, "incoming particle "); 883 theProjectile4Momentum += (*iter)->Ge 825 theProjectile4Momentum += (*iter)->GetTrackingMomentum(); 884 projectileA += (*iter)->GetDefinition 826 projectileA += (*iter)->GetDefinition()->GetBaryonNumber(); 885 projectileZ += G4lrint((*iter)->GetDe 827 projectileZ += G4lrint((*iter)->GetDefinition()->GetPDGCharge()/eplus); 886 #ifdef debug_BIC_Propagate 828 #ifdef debug_BIC_Propagate 887 G4cout << " Adding initial secondary 829 G4cout << " Adding initial secondary " << *iter 888 << " time" << (*iter)->GetForma 830 << " time" << (*iter)->GetFormationTime() 889 << ", state " << (*iter)->GetSt 831 << ", state " << (*iter)->GetState() << G4endl; 890 #endif 832 #endif 891 } 833 } 892 } 834 } 893 //theCollisionMgr->Print(); 835 //theCollisionMgr->Print(); 894 const G4HadProjectile * primary = GetPrimar << 895 836 >> 837 const G4HadProjectile * primary = GetPrimaryProjectile(); // check for primary from TheoHE model 896 if (primary){ 838 if (primary){ 897 G4LorentzVector mom=primary->Get4Momentu 839 G4LorentzVector mom=primary->Get4Momentum(); 898 theProjectile4Momentum += mom; 840 theProjectile4Momentum += mom; 899 projectileA = primary->GetDefinition()-> 841 projectileA = primary->GetDefinition()->GetBaryonNumber(); 900 projectileZ = G4lrint(primary->GetDefini 842 projectileZ = G4lrint(primary->GetDefinition()->GetPDGCharge()/eplus); 901 // now check if "excitation" energy left 843 // now check if "excitation" energy left by TheoHE model 902 G4double excitation= theProjectile4Momen 844 G4double excitation= theProjectile4Momentum.e() + initial_nuclear_mass - lateParticles4Momentum.e() - massInNucleus; 903 #ifdef debug_BIC_GetExcitationEnergy 845 #ifdef debug_BIC_GetExcitationEnergy 904 G4cout << "BIC: Proj.e, nucl initial, nu 846 G4cout << "BIC: Proj.e, nucl initial, nucl final, lateParticles" 905 << theProjectile4Momentum << ", " 847 << theProjectile4Momentum << ", " 906 << initial_nuclear_mass<< ", " << 848 << initial_nuclear_mass<< ", " << massInNucleus << ", " 907 << lateParticles4Momentum << G4end 849 << lateParticles4Momentum << G4endl; 908 G4cout << "BIC: Proj.e / initial excitat 850 G4cout << "BIC: Proj.e / initial excitation: " << theProjectile4Momentum.e() << " / " << excitation << G4endl; 909 #endif 851 #endif 910 success = excitation > 0; 852 success = excitation > 0; 911 #ifdef debug_G4BinaryCascade 853 #ifdef debug_G4BinaryCascade 912 if ( ! success ) { 854 if ( ! success ) { 913 G4cout << "G4BinaryCascade::BuildLate 855 G4cout << "G4BinaryCascade::BuildLateParticleCollisions(): Proj.e / initial excitation: " << theProjectile4Momentum.e() << " / " << excitation << G4endl; 914 //PrintKTVector(secondaries); 856 //PrintKTVector(secondaries); 915 } 857 } 916 #endif 858 #endif 917 } else { 859 } else { 918 // no primary from HE model -> cascade 860 // no primary from HE model -> cascade 919 success=true; 861 success=true; 920 } 862 } 921 863 922 if (success) { 864 if (success) { 923 secondaries->clear(); // Don't leave "G4 865 secondaries->clear(); // Don't leave "G4KineticTrack *"s in two vectors 924 delete secondaries; 866 delete secondaries; 925 } 867 } 926 return success; 868 return success; 927 } 869 } 928 870 929 //-------------------------------------------- 871 //---------------------------------------------------------------------------- 930 G4ReactionProductVector * G4BinaryCascade::De 872 G4ReactionProductVector * G4BinaryCascade::DeExcite() 931 //-------------------------------------------- 873 //---------------------------------------------------------------------------- 932 { 874 { 933 // find a fragment and call the precompound 875 // find a fragment and call the precompound model. 934 G4Fragment * fragment = nullptr; << 876 G4Fragment * fragment = 0; 935 G4ReactionProductVector * precompoundProduc << 877 G4ReactionProductVector * precompoundProducts = 0; 936 878 937 G4LorentzVector pFragment(0); 879 G4LorentzVector pFragment(0); 938 // G4cout << " final4mon " << GetFinal4Mome 880 // G4cout << " final4mon " << GetFinal4Momentum() /MeV << G4endl; 939 881 >> 882 // if ( ExcitationEnergy >= 0 ) // closed by Uzhi >> 883 // { // closed by Uzhi 940 fragment = FindFragments(); 884 fragment = FindFragments(); 941 << 885 if(fragment) // Uzhi 942 if(fragment) << 886 { // Uzhi 943 { << 887 if(fragment->GetA() >1.5) // Uzhi 944 if(fragment->GetA_asInt() >1) << 945 { 888 { 946 pFragment=fragment->GetMomentum(); 889 pFragment=fragment->GetMomentum(); 947 // G4cout << " going to preco with fr 890 // G4cout << " going to preco with fragment 4 mom " << pFragment << G4endl; 948 if (theDeExcitation) / 891 if (theDeExcitation) // pre-compound 949 { 892 { 950 precompoundProducts= theDeExcitati 893 precompoundProducts= theDeExcitation->DeExcite(*fragment); 951 } 894 } 952 else if (theExcitationHandler) // 895 else if (theExcitationHandler) // de-excitation 953 { 896 { 954 precompoundProducts=theExcitationH 897 precompoundProducts=theExcitationHandler->BreakItUp(*fragment); 955 } 898 } 956 899 957 } else 900 } else 958 { // f << 901 { // fragment->GetA() < 1.5, so a single proton, as a fragment must have Z>0 959 if (theTargetList.size() + theCapture 902 if (theTargetList.size() + theCapturedList.size() > 1 ) { 960 throw G4HadronicException(__FILE__ 903 throw G4HadronicException(__FILE__, __LINE__, "G4BinaryCasacde:: Invalid Fragment"); 961 } 904 } 962 905 963 std::vector<G4KineticTrack *>::iterat 906 std::vector<G4KineticTrack *>::iterator i; 964 if ( theTargetList.size() == 1 ) {i= 907 if ( theTargetList.size() == 1 ) {i=theTargetList.begin();} 965 if ( theCapturedList.size() == 1 ) {i << 908 if ( theCapturedList.size() == 1 ) {i=theCapturedList.begin();} // Uzhi 966 G4ReactionProduct * aNew = new G4Reac 909 G4ReactionProduct * aNew = new G4ReactionProduct((*i)->GetDefinition()); 967 aNew->SetTotalEnergy((*i)->GetDefinit 910 aNew->SetTotalEnergy((*i)->GetDefinition()->GetPDGMass()); 968 aNew->SetCreatorModelID(theBIC_ID); << 969 aNew->SetParentResonanceDef((*i)->GetParent << 970 aNew->SetParentResonanceID((*i)->GetParentR << 971 aNew->SetMomentum(G4ThreeVector(0));/ 911 aNew->SetMomentum(G4ThreeVector(0));// see boost for preCompoundProducts below.. 972 precompoundProducts = new G4ReactionP 912 precompoundProducts = new G4ReactionProductVector(); 973 precompoundProducts->push_back(aNew); 913 precompoundProducts->push_back(aNew); 974 } // End of f 914 } // End of fragment->GetA() < 1.5 975 delete fragment; 915 delete fragment; 976 fragment=nullptr; << 916 fragment=0; 977 917 978 } else // End of 918 } else // End of if(fragment) 979 { // No fra << 919 { // No fragment, can be neutrons only // Uzhi 980 920 981 precompoundProducts = DecayVoidNucleus() 921 precompoundProducts = DecayVoidNucleus(); 982 } 922 } 983 #ifdef debug_BIC_DeexcitationProducts << 984 << 985 G4LorentzVector fragment_momentum=GetFi << 986 G4LorentzVector Preco_momentum; << 987 if ( precompoundProducts ) << 988 { << 989 std::vector<G4ReactionProduct *>::it << 990 for(j = precompoundProducts->begin() << 991 { << 992 G4LorentzVector pProduct((*j)->Ge << 993 Preco_momentum += pProduct; << 994 } << 995 } << 996 G4cout << "finalNuclMom / sum preco pro << 997 << " delta E "<< fragment_momentum. << 998 << 999 #endif << 1000 << 1001 return precompoundProducts; 923 return precompoundProducts; 1002 } 924 } 1003 925 1004 //------------------------------------------- 926 //---------------------------------------------------------------------------- 1005 G4ReactionProductVector * G4BinaryCascade::D 927 G4ReactionProductVector * G4BinaryCascade::DecayVoidNucleus() 1006 //------------------------------------------- 928 //---------------------------------------------------------------------------- 1007 { 929 { 1008 G4ReactionProductVector * result = nullptr << 930 G4ReactionProductVector * result=0; 1009 if ( (theTargetList.size()+theCapturedList 931 if ( (theTargetList.size()+theCapturedList.size()) > 0 ) 1010 { 932 { 1011 result = new G4ReactionProductVector; 933 result = new G4ReactionProductVector; 1012 std::vector<G4KineticTrack *>::iterator 934 std::vector<G4KineticTrack *>::iterator aNuc; 1013 G4LorentzVector aVec; 935 G4LorentzVector aVec; 1014 std::vector<G4double> masses; 936 std::vector<G4double> masses; 1015 G4double sumMass(0); 937 G4double sumMass(0); 1016 938 1017 if ( theTargetList.size() != 0) << 939 if ( theTargetList.size() != 0) // Uzhi 1018 { << 940 { 1019 for ( aNuc=theTargetList.begin(); aN 941 for ( aNuc=theTargetList.begin(); aNuc != theTargetList.end(); aNuc++) 1020 { 942 { 1021 G4double mass=(*aNuc)->GetDefinit 943 G4double mass=(*aNuc)->GetDefinition()->GetPDGMass(); 1022 masses.push_back(mass); 944 masses.push_back(mass); 1023 sumMass += mass; 945 sumMass += mass; 1024 } 946 } 1025 } << 947 } // Uzhi 1026 948 1027 if ( theCapturedList.size() != 0) << 949 if ( theCapturedList.size() != 0) // Uzhi 1028 { << 950 { // Uzhi 1029 for(aNuc = theCapturedList.begin(); << 951 for(aNuc = theCapturedList.begin(); // Uzhi 1030 { << 952 aNuc != theCapturedList.end(); aNuc++) // Uzhi 1031 G4double mass=(*aNuc)->GetDefinit << 953 { // Uzhi 1032 masses.push_back(mass); << 954 G4double mass=(*aNuc)->GetDefinition()->GetPDGMass(); // Uzhi 1033 sumMass += mass; << 955 masses.push_back(mass); // Uzhi >> 956 sumMass += mass; // Uzhi 1034 } 957 } 1035 } 958 } 1036 959 1037 G4LorentzVector finalP=GetFinal4Momentu 960 G4LorentzVector finalP=GetFinal4Momentum(); 1038 G4FermiPhaseSpaceDecay decay; 961 G4FermiPhaseSpaceDecay decay; 1039 // G4cout << " some neutrons? " << mass 962 // G4cout << " some neutrons? " << masses.size() <<" " ; 1040 // G4cout<< theTargetList.size()<<" "<< 963 // G4cout<< theTargetList.size()<<" "<<finalP <<" " << finalP.mag()<<G4endl; 1041 964 1042 G4double eCMS=finalP.mag(); 965 G4double eCMS=finalP.mag(); 1043 if ( eCMS < sumMass ) 966 if ( eCMS < sumMass ) // @@GF --- Cheat!! 1044 { 967 { 1045 eCMS=sumMass + 2*MeV*masses.size(); 968 eCMS=sumMass + 2*MeV*masses.size(); 1046 finalP.setE(std::sqrt(finalP.vect(). 969 finalP.setE(std::sqrt(finalP.vect().mag2() + sqr(eCMS))); 1047 } 970 } 1048 971 1049 precompoundLorentzboost.set(finalP.boos 972 precompoundLorentzboost.set(finalP.boostVector()); 1050 std::vector<G4LorentzVector*> * momenta 973 std::vector<G4LorentzVector*> * momenta=decay.Decay(eCMS,masses); 1051 std::vector<G4LorentzVector*>::iterator 974 std::vector<G4LorentzVector*>::iterator aMom=momenta->begin(); 1052 975 1053 976 1054 if ( theTargetList.size() != 0) 977 if ( theTargetList.size() != 0) 1055 { 978 { 1056 for ( aNuc=theTargetList.begin(); 979 for ( aNuc=theTargetList.begin(); 1057 (aNuc != theTargetList.end()) 980 (aNuc != theTargetList.end()) && (aMom!=momenta->end()); 1058 aNuc++, aMom++ ) 981 aNuc++, aMom++ ) 1059 { 982 { 1060 G4ReactionProduct * aNew = new G4 983 G4ReactionProduct * aNew = new G4ReactionProduct((*aNuc)->GetDefinition()); 1061 aNew->SetTotalEnergy((*aMom)->e() 984 aNew->SetTotalEnergy((*aMom)->e()); 1062 aNew->SetMomentum((*aMom)->vect() 985 aNew->SetMomentum((*aMom)->vect()); 1063 aNew->SetCreatorModelID(theBIC_ID << 1064 aNew->SetParentResonanceDef((*aNuc)->Ge << 1065 aNew->SetParentResonanceID((*aNuc)->Get << 1066 result->push_back(aNew); 986 result->push_back(aNew); 1067 delete *aMom; << 1068 } << 1069 } << 1070 987 1071 if ( theCapturedList.size() != 0) << 1072 { << 1073 for ( aNuc=theCapturedList.begin(); << 1074 (aNuc != theCapturedList.end() << 1075 aNuc++, aMom++ ) << 1076 { << 1077 G4ReactionProduct * aNew = new G4 << 1078 aNew->SetTotalEnergy((*aMom)->e() << 1079 aNew->SetMomentum((*aMom)->vect() << 1080 aNew->SetCreatorModelID(theBIC_ID << 1081 aNew->SetParentResonanceDef((*aNu << 1082 aNew->SetParentResonanceID((*aNuc << 1083 result->push_back(aNew); << 1084 delete *aMom; 988 delete *aMom; 1085 } 989 } 1086 } 990 } 1087 991 >> 992 if ( theCapturedList.size() != 0) // Uzhi >> 993 { // Uzhi >> 994 for ( aNuc=theCapturedList.begin(); // Uzhi >> 995 (aNuc != theCapturedList.end()) && (aMom!=momenta->end()); // Uzhi >> 996 aNuc++, aMom++ ) // Uzhi >> 997 { // Uzhi >> 998 G4ReactionProduct * aNew = new G4ReactionProduct( // Uzhi >> 999 (*aNuc)->GetDefinition()); // Uzhi >> 1000 aNew->SetTotalEnergy((*aMom)->e()); // Uzhi >> 1001 aNew->SetMomentum((*aMom)->vect()); // Uzhi >> 1002 result->push_back(aNew); // Uzhi >> 1003 delete *aMom; // Uzhi >> 1004 } // Uzhi >> 1005 } // Uzhi >> 1006 1088 delete momenta; 1007 delete momenta; 1089 } 1008 } 1090 return result; 1009 return result; 1091 } // End if(!fragment) 1010 } // End if(!fragment) 1092 1011 1093 //------------------------------------------- 1012 //---------------------------------------------------------------------------- 1094 G4ReactionProductVector * G4BinaryCascade::Pr 1013 G4ReactionProductVector * G4BinaryCascade::ProductsAddFinalState(G4ReactionProductVector * products, G4KineticTrackVector & fs) 1095 //------------------------------------------- 1014 //---------------------------------------------------------------------------- 1096 { 1015 { 1097 // fill in products the outgoing particles 1016 // fill in products the outgoing particles 1098 std::size_t i(0); << 1017 size_t i(0); 1099 #ifdef debug_BIC_Propagate_finals << 1100 G4LorentzVector mom_fs; << 1101 #endif << 1102 for(i = 0; i< fs.size(); i++) 1018 for(i = 0; i< fs.size(); i++) 1103 { 1019 { 1104 G4KineticTrack * kt = fs[i]; 1020 G4KineticTrack * kt = fs[i]; 1105 G4ReactionProduct * aNew = new G4Reac 1021 G4ReactionProduct * aNew = new G4ReactionProduct(kt->GetDefinition()); 1106 aNew->SetMomentum(kt->Get4Momentum(). 1022 aNew->SetMomentum(kt->Get4Momentum().vect()); 1107 aNew->SetTotalEnergy(kt->Get4Momentum 1023 aNew->SetTotalEnergy(kt->Get4Momentum().e()); 1108 aNew->SetNewlyAdded(kt->IsParticipant 1024 aNew->SetNewlyAdded(kt->IsParticipant()); 1109 aNew->SetCreatorModelID(theBIC_ID); << 1110 aNew->SetParentResonanceDef(kt->GetPa << 1111 aNew->SetParentResonanceID(kt->GetPar << 1112 products->push_back(aNew); 1025 products->push_back(aNew); 1113 1026 1114 #ifdef debug_BIC_Propagate_finals 1027 #ifdef debug_BIC_Propagate_finals 1115 mom_fs += kt->Get4Momentum(); << 1116 G4cout <<kt->GetDefinition()->GetPart 1028 G4cout <<kt->GetDefinition()->GetParticleName(); 1117 G4cout << " Particle Ekin " << aNew-> 1029 G4cout << " Particle Ekin " << aNew->GetKineticEnergy(); 1118 G4cout << ", is " << (kt->GetDefiniti << 1030 G4cout << "final is " << kt->GetDefinition()->GetPDGStable() ? "stable" : 1119 (kt->GetDefinition()->IsShort << 1031 ( kt->GetDefinition()->IsShortLived() ? "short lived " : "non stable") << G4endl;; 1120 G4cout << G4endl; << 1121 #endif 1032 #endif 1122 1033 1123 } 1034 } 1124 #ifdef debug_BIC_Propagate_finals << 1125 G4cout << " Final state momentum " << mom << 1126 #endif << 1127 << 1128 return products; 1035 return products; 1129 } 1036 } 1130 //------------------------------------------- 1037 //---------------------------------------------------------------------------- 1131 G4ReactionProductVector * G4BinaryCascade::Pr 1038 G4ReactionProductVector * G4BinaryCascade::ProductsAddPrecompound(G4ReactionProductVector * products, G4ReactionProductVector * precompoundProducts) 1132 //------------------------------------------- 1039 //---------------------------------------------------------------------------- 1133 { 1040 { 1134 G4LorentzVector pSumPreco(0), pPreco(0); 1041 G4LorentzVector pSumPreco(0), pPreco(0); 1135 1042 1136 if ( precompoundProducts ) 1043 if ( precompoundProducts ) 1137 { 1044 { 1138 for(auto j = precompoundProducts->cbegi << 1045 std::vector<G4ReactionProduct *>::iterator j; >> 1046 for(j = precompoundProducts->begin(); j != precompoundProducts->end(); ++j) 1139 { 1047 { 1140 // boost back to system of moving nu 1048 // boost back to system of moving nucleus 1141 G4LorentzVector pProduct((*j)->GetMo 1049 G4LorentzVector pProduct((*j)->GetMomentum(),(*j)->GetTotalEnergy()); 1142 pPreco+= pProduct; 1050 pPreco+= pProduct; 1143 #ifdef debug_BIC_Propagate_finals 1051 #ifdef debug_BIC_Propagate_finals 1144 G4cout << "BIC: pProduct be4 boost " << 1052 G4cout << " pProduct be4 boost " <<pProduct << G4endl; 1145 #endif 1053 #endif 1146 pProduct *= precompoundLorentzboost; 1054 pProduct *= precompoundLorentzboost; 1147 #ifdef debug_BIC_Propagate_finals 1055 #ifdef debug_BIC_Propagate_finals 1148 G4cout << "BIC: pProduct aft boost " << 1056 G4cout << " pProduct aft boost " <<pProduct << G4endl; 1149 #endif 1057 #endif 1150 pSumPreco += pProduct; 1058 pSumPreco += pProduct; 1151 (*j)->SetTotalEnergy(pProduct.e()); 1059 (*j)->SetTotalEnergy(pProduct.e()); 1152 (*j)->SetMomentum(pProduct.vect()); 1060 (*j)->SetMomentum(pProduct.vect()); 1153 (*j)->SetNewlyAdded(true); 1061 (*j)->SetNewlyAdded(true); 1154 products->push_back(*j); 1062 products->push_back(*j); 1155 } 1063 } 1156 // G4cout << " unboosted preco result m 1064 // G4cout << " unboosted preco result mom " << pPreco / MeV << " ..- fragmentMom " << (pPreco - pFragment)/MeV<< G4endl; 1157 // G4cout << " preco result mom " << pS 1065 // G4cout << " preco result mom " << pSumPreco / MeV << " ..-file4Mom " << (pSumPreco - GetFinal4Momentum())/MeV<< G4endl; 1158 precompoundProducts->clear(); 1066 precompoundProducts->clear(); 1159 delete precompoundProducts; 1067 delete precompoundProducts; 1160 } 1068 } 1161 return products; 1069 return products; 1162 } 1070 } 1163 //------------------------------------------- 1071 //---------------------------------------------------------------------------- 1164 void G4BinaryCascade::FindCollisions(G4Kinet 1072 void G4BinaryCascade::FindCollisions(G4KineticTrackVector * secondaries) 1165 //------------------------------------------- 1073 //---------------------------------------------------------------------------- 1166 { 1074 { 1167 for(auto i=secondaries->cbegin(); i!=seco << 1075 for(std::vector<G4KineticTrack *>::iterator i = secondaries->begin(); >> 1076 i != secondaries->end(); ++i) 1168 { 1077 { 1169 for(auto j=theImR.cbegin(); j!=theImR << 1078 >> 1079 for(std::vector<G4BCAction *>::iterator j = theImR.begin(); >> 1080 j!=theImR.end(); j++) 1170 { 1081 { 1171 // G4cout << "G4BinaryCascad 1082 // G4cout << "G4BinaryCascade::FindCollisions: at action " << *j << G4endl; 1172 const std::vector<G4CollisionInit 1083 const std::vector<G4CollisionInitialState *> & aCandList 1173 = (*j)->GetCollisions(*i, theTarg 1084 = (*j)->GetCollisions(*i, theTargetList, theCurrentTime); 1174 for(std::size_t count=0; count<aC << 1085 for(size_t count=0; count<aCandList.size(); count++) 1175 { 1086 { 1176 theCollisionMgr->AddCollision 1087 theCollisionMgr->AddCollision(aCandList[count]); 1177 //4cout << "================= 1088 //4cout << "====================== New Collision ================="<<G4endl; 1178 //theCollisionMgr->Print(); 1089 //theCollisionMgr->Print(); 1179 } 1090 } 1180 } 1091 } 1181 } 1092 } 1182 } 1093 } 1183 1094 1184 1095 1185 //------------------------------------------- 1096 //---------------------------------------------------------------------------- 1186 void G4BinaryCascade::FindDecayCollision(G4K 1097 void G4BinaryCascade::FindDecayCollision(G4KineticTrack * secondary) 1187 //------------------------------------------- 1098 //---------------------------------------------------------------------------- 1188 { 1099 { 1189 const auto& aCandList = theDecay->GetColl << 1100 const std::vector<G4CollisionInitialState *> & aCandList 1190 for(std::size_t count=0; count<aCandList. << 1101 = theDecay->GetCollisions(secondary, theTargetList, theCurrentTime); >> 1102 for(size_t count=0; count<aCandList.size(); count++) 1191 { 1103 { 1192 theCollisionMgr->AddCollision(aCandLi 1104 theCollisionMgr->AddCollision(aCandList[count]); 1193 } 1105 } 1194 } 1106 } 1195 1107 1196 //------------------------------------------- 1108 //---------------------------------------------------------------------------- 1197 void G4BinaryCascade::FindLateParticleCollis 1109 void G4BinaryCascade::FindLateParticleCollision(G4KineticTrack * secondary) 1198 //------------------------------------------- 1110 //---------------------------------------------------------------------------- 1199 { 1111 { 1200 1112 1201 G4double tin=0., tout=0.; 1113 G4double tin=0., tout=0.; 1202 if (((G4RKPropagation*)thePropagator)->Ge 1114 if (((G4RKPropagation*)thePropagator)->GetSphereIntersectionTimes(secondary,tin,tout)) 1203 { 1115 { 1204 if ( tin > 0 ) 1116 if ( tin > 0 ) 1205 { 1117 { 1206 secondary->SetState(G4KineticTrac 1118 secondary->SetState(G4KineticTrack::outside); 1207 } else if ( tout > 0 ) 1119 } else if ( tout > 0 ) 1208 { 1120 { 1209 secondary->SetState(G4KineticTrac 1121 secondary->SetState(G4KineticTrack::inside); 1210 } else { 1122 } else { 1211 //G4cout << "G4BC set miss , tin, 1123 //G4cout << "G4BC set miss , tin, tout " << tin << " , " << tout <<G4endl; 1212 secondary->SetState(G4KineticTrac 1124 secondary->SetState(G4KineticTrack::miss_nucleus); 1213 } 1125 } 1214 } else { 1126 } else { 1215 secondary->SetState(G4KineticTrack::m 1127 secondary->SetState(G4KineticTrack::miss_nucleus); 1216 //G4cout << "G4BC set miss ,no inters 1128 //G4cout << "G4BC set miss ,no intersect tin, tout " << tin << " , " << tout <<G4endl; 1217 } 1129 } 1218 1130 1219 1131 1220 #ifdef debug_BIC_FindCollision 1132 #ifdef debug_BIC_FindCollision 1221 G4cout << "FindLateP Particle, 4-mom, tim 1133 G4cout << "FindLateP Particle, 4-mom, times newState " 1222 << secondary->GetDefinition()->Ge 1134 << secondary->GetDefinition()->GetParticleName() << " " 1223 << secondary->Get4Momentum() 1135 << secondary->Get4Momentum() 1224 << " times " << tin << " " << to 1136 << " times " << tin << " " << tout << " " 1225 << secondary->GetState() << G4end 1137 << secondary->GetState() << G4endl; 1226 #endif 1138 #endif 1227 1139 1228 const auto& aCandList = theLateParticle-> << 1140 const std::vector<G4CollisionInitialState *> & aCandList 1229 for(std::size_t count=0; count<aCandList. << 1141 = theLateParticle->GetCollisions(secondary, theTargetList, theCurrentTime); >> 1142 for(size_t count=0; count<aCandList.size(); count++) 1230 { 1143 { 1231 #ifdef debug_BIC_FindCollision 1144 #ifdef debug_BIC_FindCollision 1232 G4cout << " Adding a late Col : " << 1145 G4cout << " Adding a late Col : " << aCandList[count] << G4endl; 1233 #endif 1146 #endif 1234 theCollisionMgr->AddCollision(aCandLi 1147 theCollisionMgr->AddCollision(aCandList[count]); 1235 } 1148 } 1236 } 1149 } 1237 1150 1238 1151 1239 //------------------------------------------- 1152 //---------------------------------------------------------------------------- 1240 G4bool G4BinaryCascade::ApplyCollision(G4Coll 1153 G4bool G4BinaryCascade::ApplyCollision(G4CollisionInitialState * collision) 1241 //------------------------------------------- 1154 //---------------------------------------------------------------------------- 1242 { 1155 { 1243 G4KineticTrack * primary = collision->Get 1156 G4KineticTrack * primary = collision->GetPrimary(); 1244 1157 1245 #ifdef debug_BIC_ApplyCollision 1158 #ifdef debug_BIC_ApplyCollision 1246 G4cerr << "G4BinaryCascade::ApplyCollisio 1159 G4cerr << "G4BinaryCascade::ApplyCollision start"<<G4endl; 1247 theCollisionMgr->Print(); 1160 theCollisionMgr->Print(); 1248 G4cout << "ApplyCollisions : projte 4mom 1161 G4cout << "ApplyCollisions : projte 4mom " << primary->GetTrackingMomentum()<< G4endl; 1249 #endif 1162 #endif 1250 1163 1251 G4KineticTrackVector target_collection=co 1164 G4KineticTrackVector target_collection=collision->GetTargetCollection(); 1252 G4bool haveTarget=target_collection.size( 1165 G4bool haveTarget=target_collection.size()>0; 1253 if( haveTarget && (primary->GetState() != 1166 if( haveTarget && (primary->GetState() != G4KineticTrack::inside) ) 1254 { 1167 { 1255 #ifdef debug_G4BinaryCascade 1168 #ifdef debug_G4BinaryCascade 1256 G4cout << "G4BinaryCasacde::ApplyColl 1169 G4cout << "G4BinaryCasacde::ApplyCollision(): StateError " << primary << G4endl; 1257 PrintKTVector(primary,std::string("pr 1170 PrintKTVector(primary,std::string("primay- ...")); 1258 PrintKTVector(&target_collection,std: 1171 PrintKTVector(&target_collection,std::string("... targets")); 1259 collision->Print(); 1172 collision->Print(); 1260 G4cout << G4endl; 1173 G4cout << G4endl; 1261 theCollisionMgr->Print(); 1174 theCollisionMgr->Print(); 1262 //*GF* throw G4HadronicException( 1175 //*GF* throw G4HadronicException(__FILE__, __LINE__, "G4BinaryCasacde::ApplyCollision()"); 1263 #endif 1176 #endif 1264 return false; 1177 return false; 1265 // } else { 1178 // } else { 1266 // G4cout << "G4BinaryCasacde::ApplyCol 1179 // G4cout << "G4BinaryCasacde::ApplyCollision(): decay " << G4endl; 1267 // PrintKTVector(primary,std::string("p 1180 // PrintKTVector(primary,std::string("primay- ...")); 1268 // G4double tin=0., tout=0.; 1181 // G4double tin=0., tout=0.; 1269 // if (((G4RKPropagation*)thePropagator 1182 // if (((G4RKPropagation*)thePropagator)->GetSphereIntersectionTimes(primary,tin,tout)) 1270 // { 1183 // { 1271 // G4cout << "tin tout: " << tin << 1184 // G4cout << "tin tout: " << tin << " " << tout << G4endl; 1272 // } 1185 // } 1273 1186 1274 } 1187 } 1275 1188 1276 G4LorentzVector mom4Primary=primary->Get4 1189 G4LorentzVector mom4Primary=primary->Get4Momentum(); 1277 1190 1278 G4int initialBaryon(0); 1191 G4int initialBaryon(0); 1279 G4int initialCharge(0); 1192 G4int initialCharge(0); 1280 if (primary->GetState() == G4KineticTrack 1193 if (primary->GetState() == G4KineticTrack::inside) 1281 { 1194 { 1282 initialBaryon = primary->GetDefinitio 1195 initialBaryon = primary->GetDefinition()->GetBaryonNumber(); 1283 initialCharge = G4lrint(primary->GetD 1196 initialCharge = G4lrint(primary->GetDefinition()->GetPDGCharge()/eplus); 1284 } 1197 } 1285 1198 1286 // for primary resonances, subtract neutr 1199 // for primary resonances, subtract neutron ( = proton) field ( ie. add std::abs(field)) 1287 G4double initial_Efermi=CorrectShortlived << 1200 G4double initial_Efermi=CorrectShortlivedPrimaryForFermi(primary,target_collection); 1288 //*************************************** 1201 //**************************************** 1289 1202 1290 1203 1291 G4KineticTrackVector * products = collisi 1204 G4KineticTrackVector * products = collision->GetFinalState(); 1292 1205 1293 #ifdef debug_BIC_ApplyCollision 1206 #ifdef debug_BIC_ApplyCollision 1294 DebugApplyCollisionFail(collision, produc 1207 DebugApplyCollisionFail(collision, products); 1295 #endif 1208 #endif 1296 1209 1297 // reset primary to initial state, in cas 1210 // reset primary to initial state, in case there is a veto... 1298 primary->Set4Momentum(mom4Primary); 1211 primary->Set4Momentum(mom4Primary); 1299 1212 1300 G4bool lateParticleCollision= (!haveTarge 1213 G4bool lateParticleCollision= (!haveTarget) && products && products->size() == 1; 1301 G4bool decayCollision= (!haveTarget) && p 1214 G4bool decayCollision= (!haveTarget) && products && products->size() > 1; 1302 G4bool Success(true); 1215 G4bool Success(true); 1303 1216 1304 1217 1305 #ifdef debug_G4BinaryCascade 1218 #ifdef debug_G4BinaryCascade 1306 G4int lateBaryon(0), lateCharge(0); 1219 G4int lateBaryon(0), lateCharge(0); 1307 #endif 1220 #endif 1308 1221 1309 if ( lateParticleCollision ) 1222 if ( lateParticleCollision ) 1310 { // for late particles, reset charges 1223 { // for late particles, reset charges 1311 //G4cout << "lateP, initial B C state 1224 //G4cout << "lateP, initial B C state " << initialBaryon << " " 1312 // << initialCharge<< " " << p 1225 // << initialCharge<< " " << primary->GetState() << " "<< primary->GetDefinition()->GetParticleName()<< G4endl; 1313 #ifdef debug_G4BinaryCascade 1226 #ifdef debug_G4BinaryCascade 1314 lateBaryon = initialBaryon; 1227 lateBaryon = initialBaryon; 1315 lateCharge = initialCharge; 1228 lateCharge = initialCharge; 1316 #endif 1229 #endif 1317 initialBaryon=initialCharge=0; 1230 initialBaryon=initialCharge=0; 1318 lateA -= primary->GetDefinition()->Ge 1231 lateA -= primary->GetDefinition()->GetBaryonNumber(); 1319 lateZ -= G4lrint(primary->GetDefiniti 1232 lateZ -= G4lrint(primary->GetDefinition()->GetPDGCharge()/eplus); 1320 } 1233 } 1321 1234 1322 initialBaryon += collision->GetTargetBary 1235 initialBaryon += collision->GetTargetBaryonNumber(); 1323 initialCharge += G4lrint(collision->GetTa 1236 initialCharge += G4lrint(collision->GetTargetCharge()); 1324 if (!lateParticleCollision) 1237 if (!lateParticleCollision) 1325 { 1238 { 1326 if( !products || products->size()==0 | 1239 if( !products || products->size()==0 || !CheckPauliPrinciple(products) ) 1327 { 1240 { 1328 #ifdef debug_BIC_ApplyCollision 1241 #ifdef debug_BIC_ApplyCollision 1329 if (products) G4cout << " ======Fai 1242 if (products) G4cout << " ======Failed Pauli =====" << G4endl; 1330 G4cerr << "G4BinaryCascade::ApplyCo 1243 G4cerr << "G4BinaryCascade::ApplyCollision blocked"<<G4endl; 1331 #endif 1244 #endif 1332 Success=false; 1245 Success=false; 1333 } 1246 } 1334 1247 1335 1248 1336 1249 1337 if (Success && primary->GetState() == 1250 if (Success && primary->GetState() == G4KineticTrack::inside ) { // if the primary was outside, nothing to correct 1338 if (! CorrectShortlivedFinalsForFer 1251 if (! CorrectShortlivedFinalsForFermi(products, initial_Efermi)){ 1339 Success=false; 1252 Success=false; 1340 } 1253 } 1341 } 1254 } 1342 } 1255 } 1343 1256 1344 #ifdef debug_BIC_ApplyCollision 1257 #ifdef debug_BIC_ApplyCollision 1345 DebugApplyCollision(collision, products); 1258 DebugApplyCollision(collision, products); 1346 #endif 1259 #endif 1347 1260 1348 if ( ! Success ){ 1261 if ( ! Success ){ 1349 if (products) ClearAndDestroy(product 1262 if (products) ClearAndDestroy(products); 1350 if ( decayCollision ) FindDecayCollis 1263 if ( decayCollision ) FindDecayCollision(primary); // for decay, sample new decay 1351 delete products; 1264 delete products; 1352 products=nullptr; << 1265 products=0; 1353 return false; 1266 return false; 1354 } 1267 } 1355 1268 1356 G4int finalBaryon(0); 1269 G4int finalBaryon(0); 1357 G4int finalCharge(0); 1270 G4int finalCharge(0); 1358 G4KineticTrackVector toFinalState; 1271 G4KineticTrackVector toFinalState; 1359 for(auto i=products->cbegin(); i!=product << 1272 for(std::vector<G4KineticTrack *>::iterator i =products->begin(); i != products->end(); i++) 1360 { 1273 { 1361 if ( ! lateParticleCollision ) 1274 if ( ! lateParticleCollision ) 1362 { 1275 { 1363 (*i)->SetState(primary->GetState( 1276 (*i)->SetState(primary->GetState()); // decay may be anywhere! 1364 if ( (*i)->GetState() == G4Kineti 1277 if ( (*i)->GetState() == G4KineticTrack::inside ){ 1365 finalBaryon+=(*i)->GetDefinit 1278 finalBaryon+=(*i)->GetDefinition()->GetBaryonNumber(); 1366 finalCharge+=G4lrint((*i)->Ge 1279 finalCharge+=G4lrint((*i)->GetDefinition()->GetPDGCharge()/eplus); 1367 } else { 1280 } else { 1368 G4double tin=0., tout=0.; 1281 G4double tin=0., tout=0.; 1369 if (((G4RKPropagation*)theProp 1282 if (((G4RKPropagation*)thePropagator)->GetSphereIntersectionTimes((*i),tin,tout) && 1370 tin < 0 && tout > 0 ) 1283 tin < 0 && tout > 0 ) 1371 { 1284 { 1372 PrintKTVector((*i),"particl 1285 PrintKTVector((*i),"particle inside marked not-inside"); 1373 G4cout << "tin tout: " << 1286 G4cout << "tin tout: " << tin << " " << tout << G4endl; 1374 } 1287 } 1375 } 1288 } 1376 } else { 1289 } else { 1377 G4double tin=0., tout=0.; 1290 G4double tin=0., tout=0.; 1378 if (((G4RKPropagation*)thePropaga 1291 if (((G4RKPropagation*)thePropagator)->GetSphereIntersectionTimes((*i),tin,tout)) 1379 { 1292 { 1380 //G4cout << "tin tout: " << t 1293 //G4cout << "tin tout: " << tin << " " << tout << G4endl; 1381 if ( tin > 0 ) 1294 if ( tin > 0 ) 1382 { 1295 { 1383 (*i)->SetState(G4KineticT 1296 (*i)->SetState(G4KineticTrack::outside); 1384 } 1297 } 1385 else if ( tout > 0 ) 1298 else if ( tout > 0 ) 1386 { 1299 { 1387 (*i)->SetState(G4KineticT 1300 (*i)->SetState(G4KineticTrack::inside); 1388 finalBaryon+=(*i)->GetDef 1301 finalBaryon+=(*i)->GetDefinition()->GetBaryonNumber(); 1389 finalCharge+=G4lrint((*i) 1302 finalCharge+=G4lrint((*i)->GetDefinition()->GetPDGCharge()/eplus); 1390 } 1303 } 1391 else 1304 else 1392 { 1305 { 1393 (*i)->SetState(G4KineticT 1306 (*i)->SetState(G4KineticTrack::gone_out); 1394 toFinalState.push_back((* 1307 toFinalState.push_back((*i)); 1395 } 1308 } 1396 } else 1309 } else 1397 { 1310 { 1398 (*i)->SetState(G4KineticTrack 1311 (*i)->SetState(G4KineticTrack::miss_nucleus); 1399 //G4cout << " G4BC - miss -la 1312 //G4cout << " G4BC - miss -late Part- no intersection found " << G4endl; 1400 toFinalState.push_back((*i)); 1313 toFinalState.push_back((*i)); 1401 } 1314 } 1402 1315 1403 } 1316 } 1404 } 1317 } 1405 if(!toFinalState.empty()) 1318 if(!toFinalState.empty()) 1406 { 1319 { 1407 theFinalState.insert(theFinalState.ce << 1320 theFinalState.insert(theFinalState.end(), 1408 toFinalState.cbegin(),toFinal << 1321 toFinalState.begin(),toFinalState.end()); 1409 std::vector<G4KineticTrack *>::iterat << 1322 std::vector<G4KineticTrack *>::iterator iter1, iter2; 1410 for(auto iter1=toFinalState.cbegin(); << 1323 for(iter1 = toFinalState.begin(); iter1 != toFinalState.end(); 1411 { << 1324 ++iter1) 1412 iter2 = std::find(products->begin << 1325 { 1413 if ( iter2 != products->cend() ) << 1326 iter2 = std::find(products->begin(), products->end(), >> 1327 *iter1); >> 1328 if ( iter2 != products->end() ) products->erase(iter2); 1414 } 1329 } 1415 theCollisionMgr->RemoveTracksCollisio 1330 theCollisionMgr->RemoveTracksCollisions(&toFinalState); 1416 } 1331 } 1417 1332 1418 //G4cout << " currentA, Z be4: " << curre 1333 //G4cout << " currentA, Z be4: " << currentA << " " << currentZ << G4endl; 1419 currentA += finalBaryon-initialBaryon; 1334 currentA += finalBaryon-initialBaryon; 1420 currentZ += finalCharge-initialCharge; 1335 currentZ += finalCharge-initialCharge; 1421 //G4cout << " ApplyCollision currentA, Z 1336 //G4cout << " ApplyCollision currentA, Z aft: " << currentA << " " << currentZ << G4endl; 1422 1337 1423 G4KineticTrackVector oldSecondaries; 1338 G4KineticTrackVector oldSecondaries; 1424 oldSecondaries.push_back(primary); 1339 oldSecondaries.push_back(primary); 1425 primary->Hit(); 1340 primary->Hit(); 1426 1341 1427 #ifdef debug_G4BinaryCascade 1342 #ifdef debug_G4BinaryCascade 1428 if ( (finalBaryon-initialBaryon-lateBaryo 1343 if ( (finalBaryon-initialBaryon-lateBaryon) != 0 || (finalCharge-initialCharge-lateCharge) != 0 ) 1429 { 1344 { 1430 G4cout << "G4BinaryCascade: Error in 1345 G4cout << "G4BinaryCascade: Error in Balancing: " << G4endl; 1431 G4cout << "initial/final baryon numbe 1346 G4cout << "initial/final baryon number, initial/final Charge " 1432 << initialBaryon <<" "<< fina 1347 << initialBaryon <<" "<< finalBaryon <<" " 1433 << initialCharge <<" "<< fina 1348 << initialCharge <<" "<< finalCharge <<" " 1434 << " in Collision type: "<< t 1349 << " in Collision type: "<< typeid(*collision->GetGenerator()).name() 1435 << ", with number of products 1350 << ", with number of products: "<< products->size() <<G4endl; 1436 G4cout << G4endl<<"Initial condition 1351 G4cout << G4endl<<"Initial condition are these:"<<G4endl; 1437 G4cout << "proj: "<<collision->GetPri 1352 G4cout << "proj: "<<collision->GetPrimary()->GetDefinition()->GetParticleName()<<G4endl; 1438 for(std::size_t it=0; it<collision->G << 1353 for(size_t it=0; it<collision->GetTargetCollection().size(); it++) 1439 { 1354 { 1440 G4cout << "targ: " 1355 G4cout << "targ: " 1441 <<collision->GetTargetCol 1356 <<collision->GetTargetCollection()[it]->GetDefinition()->GetParticleName()<<G4endl; 1442 } 1357 } 1443 PrintKTVector(&collision->GetTargetCo 1358 PrintKTVector(&collision->GetTargetCollection(),std::string(" Target particles")); 1444 G4cout << G4endl<<G4endl; 1359 G4cout << G4endl<<G4endl; 1445 } 1360 } 1446 #endif 1361 #endif 1447 1362 1448 G4KineticTrackVector oldTarget = collisio 1363 G4KineticTrackVector oldTarget = collision->GetTargetCollection(); 1449 for(std::size_t ii=0; ii< oldTarget.size( << 1364 for(size_t ii=0; ii< oldTarget.size(); ii++) 1450 { 1365 { 1451 oldTarget[ii]->Hit(); 1366 oldTarget[ii]->Hit(); 1452 } 1367 } 1453 1368 1454 UpdateTracksAndCollisions(&oldSecondaries 1369 UpdateTracksAndCollisions(&oldSecondaries, &oldTarget, products); 1455 std::for_each(oldSecondaries.cbegin(), ol << 1370 std::for_each(oldSecondaries.begin(), oldSecondaries.end(), Delete<G4KineticTrack>()); 1456 std::for_each(oldTarget.cbegin(), oldTarg << 1371 std::for_each(oldTarget.begin(), oldTarget.end(), Delete<G4KineticTrack>()); 1457 1372 1458 delete products; 1373 delete products; 1459 return true; 1374 return true; 1460 } 1375 } 1461 1376 1462 1377 1463 //------------------------------------------- 1378 //---------------------------------------------------------------------------- 1464 G4bool G4BinaryCascade::Absorb() 1379 G4bool G4BinaryCascade::Absorb() 1465 //------------------------------------------- 1380 //---------------------------------------------------------------------------- 1466 { 1381 { 1467 // Do it in two step: cannot change theSe 1382 // Do it in two step: cannot change theSecondaryList inside the first loop. 1468 G4Absorber absorber(theCutOnPAbsorb); 1383 G4Absorber absorber(theCutOnPAbsorb); 1469 1384 1470 // Build the vector of G4KineticTracks th 1385 // Build the vector of G4KineticTracks that must be absorbed 1471 G4KineticTrackVector absorbList; 1386 G4KineticTrackVector absorbList; 1472 std::vector<G4KineticTrack *>::const_iter << 1387 std::vector<G4KineticTrack *>::iterator iter; 1473 // PrintKTVector(&theSecondaryList, " te 1388 // PrintKTVector(&theSecondaryList, " testing for Absorb" ); 1474 for(iter = theSecondaryList.cbegin(); << 1389 for(iter = theSecondaryList.begin(); 1475 iter != theSecondaryList.cend(); << 1390 iter != theSecondaryList.end(); ++iter) 1476 { 1391 { 1477 G4KineticTrack * kt = *iter; 1392 G4KineticTrack * kt = *iter; 1478 if(kt->GetState() == G4KineticTrack:: 1393 if(kt->GetState() == G4KineticTrack::inside)// absorption happens only inside the nucleus 1479 { 1394 { 1480 if(absorber.WillBeAbsorbed(*kt)) 1395 if(absorber.WillBeAbsorbed(*kt)) 1481 { 1396 { 1482 absorbList.push_back(kt); 1397 absorbList.push_back(kt); 1483 } 1398 } 1484 } 1399 } 1485 } 1400 } 1486 1401 1487 if(absorbList.empty()) 1402 if(absorbList.empty()) 1488 return false; 1403 return false; 1489 1404 1490 G4KineticTrackVector toDelete; 1405 G4KineticTrackVector toDelete; 1491 for(iter = absorbList.cbegin(); iter != a << 1406 for(iter = absorbList.begin(); iter != absorbList.end(); ++iter) 1492 { 1407 { 1493 G4KineticTrack * kt = *iter; 1408 G4KineticTrack * kt = *iter; 1494 if(!absorber.FindAbsorbers(*kt, theTa 1409 if(!absorber.FindAbsorbers(*kt, theTargetList)) 1495 throw G4HadronicException(__FILE_ 1410 throw G4HadronicException(__FILE__, __LINE__, "G4BinaryCascade::Absorb(): Cannot absorb a particle."); 1496 1411 1497 if(!absorber.FindProducts(*kt)) 1412 if(!absorber.FindProducts(*kt)) 1498 throw G4HadronicException(__FILE_ 1413 throw G4HadronicException(__FILE__, __LINE__, "G4BinaryCascade::Absorb(): Cannot absorb a particle."); 1499 1414 1500 G4KineticTrackVector * products = abs 1415 G4KineticTrackVector * products = absorber.GetProducts(); 1501 G4int maxLoopCount = 1000; << 1416 // ------ debug 1502 while(!CheckPauliPrinciple(products) << 1417 G4int count = 0; >> 1418 // ------ end debug >> 1419 while(!CheckPauliPrinciple(products)) 1503 { 1420 { >> 1421 // ------ debug >> 1422 count++; >> 1423 // ------ end debug 1504 ClearAndDestroy(products); 1424 ClearAndDestroy(products); 1505 if(!absorber.FindProducts(*kt)) 1425 if(!absorber.FindProducts(*kt)) 1506 throw G4HadronicException(__F 1426 throw G4HadronicException(__FILE__, __LINE__, 1507 "G4BinaryCascade::Abs 1427 "G4BinaryCascade::Absorb(): Cannot absorb a particle."); 1508 } 1428 } 1509 if ( --maxLoopCount < 0 ) throw G4Hadro << 1510 // ------ debug 1429 // ------ debug 1511 // G4cerr << "Absorb CheckPauliPri << 1430 // G4cerr << "Absorb CheckPauliPrinciple count= " << count << G4endl; 1512 // ------ end debug 1431 // ------ end debug 1513 G4KineticTrackVector toRemove; // bu 1432 G4KineticTrackVector toRemove; // build a vector for UpdateTrack... 1514 toRemove.push_back(kt); 1433 toRemove.push_back(kt); 1515 toDelete.push_back(kt); // delete th 1434 toDelete.push_back(kt); // delete the track later 1516 G4KineticTrackVector * absorbers = ab 1435 G4KineticTrackVector * absorbers = absorber.GetAbsorbers(); 1517 UpdateTracksAndCollisions(&toRemove, 1436 UpdateTracksAndCollisions(&toRemove, absorbers, products); 1518 ClearAndDestroy(absorbers); 1437 ClearAndDestroy(absorbers); 1519 } 1438 } 1520 ClearAndDestroy(&toDelete); 1439 ClearAndDestroy(&toDelete); 1521 return true; 1440 return true; 1522 } 1441 } 1523 1442 1524 1443 1525 1444 1526 // Capture all p and n with Energy < theCutOn 1445 // Capture all p and n with Energy < theCutOnP 1527 //------------------------------------------- 1446 //---------------------------------------------------------------------------- 1528 G4bool G4BinaryCascade::Capture(G4bool verbos 1447 G4bool G4BinaryCascade::Capture(G4bool verbose) 1529 //------------------------------------------- 1448 //---------------------------------------------------------------------------- 1530 { 1449 { 1531 G4KineticTrackVector captured; 1450 G4KineticTrackVector captured; 1532 G4bool capture = false; 1451 G4bool capture = false; 1533 std::vector<G4KineticTrack *>::const_iter << 1452 std::vector<G4KineticTrack *>::iterator i; 1534 1453 1535 G4RKPropagation * RKprop=(G4RKPropagation 1454 G4RKPropagation * RKprop=(G4RKPropagation *)thePropagator; 1536 1455 1537 G4double capturedEnergy = 0; 1456 G4double capturedEnergy = 0; 1538 G4int particlesAboveCut=0; 1457 G4int particlesAboveCut=0; 1539 G4int particlesBelowCut=0; 1458 G4int particlesBelowCut=0; 1540 if ( verbose ) G4cout << " Capture: secon 1459 if ( verbose ) G4cout << " Capture: secondaries " << theSecondaryList.size() << G4endl; 1541 for(i = theSecondaryList.cbegin(); i != t << 1460 for(i = theSecondaryList.begin(); i != theSecondaryList.end(); ++i) 1542 { 1461 { 1543 G4KineticTrack * kt = *i; 1462 G4KineticTrack * kt = *i; 1544 if (verbose) G4cout << "Capture posit 1463 if (verbose) G4cout << "Capture position, radius, state " <<kt->GetPosition().mag()<<" "<<theOuterRadius<<" "<<kt->GetState()<<G4endl; 1545 if(kt->GetState() == G4KineticTrack:: 1464 if(kt->GetState() == G4KineticTrack::inside) // capture happens only inside the nucleus 1546 { 1465 { 1547 if((kt->GetDefinition() == G4Prot 1466 if((kt->GetDefinition() == G4Proton::Proton()) || 1548 (kt->GetDefinition() == G 1467 (kt->GetDefinition() == G4Neutron::Neutron())) 1549 { 1468 { 1550 //GF cut on kinetic energy 1469 //GF cut on kinetic energy if(kt->Get4Momentum().vect().mag() >= theCutOnP) 1551 G4double field=RKprop->GetFie 1470 G4double field=RKprop->GetField(kt->GetDefinition()->GetPDGEncoding(),kt->GetPosition()) 1552 -RKprop->GetBarr 1471 -RKprop->GetBarrier(kt->GetDefinition()->GetPDGEncoding()); 1553 G4double energy= kt->Get4Mome 1472 G4double energy= kt->Get4Momentum().e() - kt->GetActualMass() + field; 1554 if (verbose ) G4cout << "Capt 1473 if (verbose ) G4cout << "Capture: .e(), mass, field, energy" << kt->Get4Momentum().e() <<" "<<kt->GetActualMass()<<" "<<field<<" "<<energy<< G4endl; 1555 // if( energy < theCutOnP ) 1474 // if( energy < theCutOnP ) 1556 // { 1475 // { 1557 capturedEnergy+=energy; 1476 capturedEnergy+=energy; 1558 ++particlesBelowCut; 1477 ++particlesBelowCut; 1559 // } else 1478 // } else 1560 // { 1479 // { 1561 // ++particlesAboveCut; 1480 // ++particlesAboveCut; 1562 // } 1481 // } 1563 } 1482 } 1564 } 1483 } 1565 } 1484 } 1566 if (verbose) G4cout << "Capture particles 1485 if (verbose) G4cout << "Capture particlesAboveCut,particlesBelowCut, capturedEnergy,capturedEnergy/particlesBelowCut <? 0.2*theCutOnP " 1567 << particlesAboveCut << " " << pa 1486 << particlesAboveCut << " " << particlesBelowCut << " " << capturedEnergy 1568 << " " << G4endl; 1487 << " " << G4endl; 1569 // << " " << (particlesBelowCut>0) ? (capt 1488 // << " " << (particlesBelowCut>0) ? (capturedEnergy/particlesBelowCut) : (capturedEnergy) << " " << 0.2*theCutOnP << G4endl; 1570 // if(particlesAboveCut==0 && particlesB 1489 // if(particlesAboveCut==0 && particlesBelowCut>0 && capturedEnergy/particlesBelowCut<0.2*theCutOnP) 1571 if(particlesBelowCut>0 && capturedEnergy/ 1490 if(particlesBelowCut>0 && capturedEnergy/particlesBelowCut<0.2*theCutOnP) 1572 { 1491 { 1573 capture=true; 1492 capture=true; 1574 for(i = theSecondaryList.cbegin(); i << 1493 for(i = theSecondaryList.begin(); i != theSecondaryList.end(); ++i) 1575 { 1494 { 1576 G4KineticTrack * kt = *i; 1495 G4KineticTrack * kt = *i; 1577 if(kt->GetState() == G4KineticTra 1496 if(kt->GetState() == G4KineticTrack::inside) // capture happens only inside the nucleus 1578 { 1497 { 1579 if((kt->GetDefinition() == G4 1498 if((kt->GetDefinition() == G4Proton::Proton()) || 1580 (kt->GetDefinition() 1499 (kt->GetDefinition() == G4Neutron::Neutron())) 1581 { 1500 { 1582 captured.push_back(kt); 1501 captured.push_back(kt); 1583 kt->Hit(); // 1502 kt->Hit(); // 1584 theCapturedList.push_back 1503 theCapturedList.push_back(kt); 1585 } 1504 } 1586 } 1505 } 1587 } 1506 } 1588 UpdateTracksAndCollisions(&captured, << 1507 UpdateTracksAndCollisions(&captured, NULL, NULL); 1589 } 1508 } 1590 1509 1591 return capture; 1510 return capture; 1592 } 1511 } 1593 1512 1594 1513 1595 //------------------------------------------- 1514 //---------------------------------------------------------------------------- 1596 G4bool G4BinaryCascade::CheckPauliPrinciple(G 1515 G4bool G4BinaryCascade::CheckPauliPrinciple(G4KineticTrackVector * products) 1597 //------------------------------------------- 1516 //---------------------------------------------------------------------------- 1598 { 1517 { 1599 G4int A = the3DNucleus->GetMassNumber(); 1518 G4int A = the3DNucleus->GetMassNumber(); 1600 G4int Z = the3DNucleus->GetCharge(); 1519 G4int Z = the3DNucleus->GetCharge(); 1601 1520 1602 G4FermiMomentum fermiMom; 1521 G4FermiMomentum fermiMom; 1603 fermiMom.Init(A, Z); 1522 fermiMom.Init(A, Z); 1604 1523 1605 const G4VNuclearDensity * density=the3DNu 1524 const G4VNuclearDensity * density=the3DNucleus->GetNuclearDensity(); 1606 1525 1607 G4KineticTrackVector::const_iterator i; << 1526 G4KineticTrackVector::iterator i; 1608 const G4ParticleDefinition * definition; << 1527 G4ParticleDefinition * definition; 1609 1528 1610 // ------ debug 1529 // ------ debug 1611 G4bool myflag = true; 1530 G4bool myflag = true; 1612 // ------ end debug 1531 // ------ end debug 1613 // G4ThreeVector xpos(0); 1532 // G4ThreeVector xpos(0); 1614 for(i = products->cbegin(); i != products << 1533 for(i = products->begin(); i != products->end(); ++i) 1615 { 1534 { 1616 definition = (*i)->GetDefinition(); 1535 definition = (*i)->GetDefinition(); 1617 if((definition == G4Proton::Proton()) 1536 if((definition == G4Proton::Proton()) || 1618 (definition == G4Neutron::Neu 1537 (definition == G4Neutron::Neutron())) 1619 { 1538 { 1620 G4ThreeVector pos = (*i)->GetPosi 1539 G4ThreeVector pos = (*i)->GetPosition(); 1621 G4double d = density->GetDensity( 1540 G4double d = density->GetDensity(pos); 1622 // energy correspondiing to fermi 1541 // energy correspondiing to fermi momentum 1623 G4double eFermi = std::sqrt( sqr( 1542 G4double eFermi = std::sqrt( sqr(fermiMom.GetFermiMomentum(d)) + (*i)->Get4Momentum().mag2() ); 1624 if( definition == G4Proton::Proto 1543 if( definition == G4Proton::Proton() ) 1625 { 1544 { 1626 eFermi -= the3DNucleus->Coulo 1545 eFermi -= the3DNucleus->CoulombBarrier(); 1627 } 1546 } 1628 G4LorentzVector mom = (*i)->Get4M 1547 G4LorentzVector mom = (*i)->Get4Momentum(); 1629 // ------ debug 1548 // ------ debug 1630 /* 1549 /* 1631 * G4cout << "p:[" << (1/M 1550 * G4cout << "p:[" << (1/MeV)*mom.x() << " " << (1/MeV)*mom.y() << " " 1632 * << (1/MeV)*mom.z() 1551 * << (1/MeV)*mom.z() << "] |p3|:" 1633 * << (1/MeV)*mom.vect 1552 * << (1/MeV)*mom.vect().mag() << " E:" << (1/MeV)*mom.t() << " ] m: " 1634 * << (1/MeV)*mom.mag( 1553 * << (1/MeV)*mom.mag() << " pos[" 1635 * << (1/fermi)*pos.x( 1554 * << (1/fermi)*pos.x() << " "<< (1/fermi)*pos.y() << " " 1636 * << (1/fermi)*pos.z( 1555 * << (1/fermi)*pos.z() << "] |Dpos|: " 1637 * << (1/fermi)*(pos-x 1556 * << (1/fermi)*(pos-xpos).mag() << " Pfermi:" 1638 * << (1/MeV)*p << G4e 1557 * << (1/MeV)*p << G4endl; 1639 * xpos=pos; 1558 * xpos=pos; 1640 */ 1559 */ 1641 // ------ end debug 1560 // ------ end debug 1642 if(mom.e() < eFermi ) 1561 if(mom.e() < eFermi ) 1643 { 1562 { 1644 // ------ debug 1563 // ------ debug 1645 myflag = false; 1564 myflag = false; 1646 // ------ end debug 1565 // ------ end debug 1647 // return false; 1566 // return false; 1648 } 1567 } 1649 } 1568 } 1650 } 1569 } 1651 #ifdef debug_BIC_CheckPauli 1570 #ifdef debug_BIC_CheckPauli 1652 if ( myflag ) 1571 if ( myflag ) 1653 { 1572 { 1654 for(i = products->cbegin(); i != prod << 1573 for(i = products->begin(); i != products->end(); ++i) 1655 { 1574 { 1656 definition = (*i)->GetDefinition( 1575 definition = (*i)->GetDefinition(); 1657 if((definition == G4Proton::Proto 1576 if((definition == G4Proton::Proton()) || 1658 (definition == G4Neutron: 1577 (definition == G4Neutron::Neutron())) 1659 { 1578 { 1660 G4ThreeVector pos = (*i)->Get 1579 G4ThreeVector pos = (*i)->GetPosition(); 1661 G4double d = density->GetDens 1580 G4double d = density->GetDensity(pos); 1662 G4double pFermi = fermiMom.Ge 1581 G4double pFermi = fermiMom.GetFermiMomentum(d); 1663 G4LorentzVector mom = (*i)->G 1582 G4LorentzVector mom = (*i)->Get4Momentum(); 1664 G4double field =((G4RKPropaga 1583 G4double field =((G4RKPropagation*)thePropagator)->GetField(definition->GetPDGEncoding(),pos); 1665 if ( mom.e()-mom.mag()+field 1584 if ( mom.e()-mom.mag()+field > 160*MeV ) 1666 { 1585 { 1667 G4cout << "momentum probl 1586 G4cout << "momentum problem pFermi=" << pFermi 1668 << " mom, mom.m " 1587 << " mom, mom.m " << mom << " " << mom.mag() 1669 << " field " << f 1588 << " field " << field << G4endl; 1670 } 1589 } 1671 } 1590 } 1672 } 1591 } 1673 } 1592 } 1674 #endif 1593 #endif 1675 1594 1676 return myflag; 1595 return myflag; 1677 } 1596 } 1678 1597 1679 //------------------------------------------- 1598 //---------------------------------------------------------------------------- 1680 void G4BinaryCascade::StepParticlesOut() 1599 void G4BinaryCascade::StepParticlesOut() 1681 //------------------------------------------- 1600 //---------------------------------------------------------------------------- 1682 { 1601 { 1683 G4int counter=0; 1602 G4int counter=0; 1684 G4int countreset=0; 1603 G4int countreset=0; 1685 //G4cout << " nucl. Radius " << radius << 1604 //G4cout << " nucl. Radius " << radius << G4endl; 1686 // G4cerr <<"pre-while- theSecondaryList 1605 // G4cerr <<"pre-while- theSecondaryList "<<G4endl; 1687 while( theSecondaryList.size() > 0 ) << 1606 while( theSecondaryList.size() > 0 ) 1688 << 1689 { 1607 { >> 1608 G4int nsec=0; 1690 G4double minTimeStep = 1.e-12*ns; / 1609 G4double minTimeStep = 1.e-12*ns; // about 30*fermi/(0.1*c_light);1.e-12*ns 1691 / << 1610 // i.e. a big step 1692 for(auto i = theSecondaryList.cbegin( << 1611 std::vector<G4KineticTrack *>::iterator i; >> 1612 for(i = theSecondaryList.begin(); i != theSecondaryList.end(); ++i) 1693 { 1613 { 1694 G4KineticTrack * kt = *i; 1614 G4KineticTrack * kt = *i; 1695 if( kt->GetState() == G4KineticTr 1615 if( kt->GetState() == G4KineticTrack::inside ) 1696 { 1616 { >> 1617 nsec++; 1697 G4double tStep(0), tdummy(0); 1618 G4double tStep(0), tdummy(0); 1698 G4bool intersect = 1619 G4bool intersect = 1699 ((G4RKPropagation*)th 1620 ((G4RKPropagation*)thePropagator)->GetSphereIntersectionTimes(kt,tdummy,tStep); 1700 #ifdef debug_BIC_StepParticlesOut 1621 #ifdef debug_BIC_StepParticlesOut 1701 G4cout << " minTimeStep, tSte 1622 G4cout << " minTimeStep, tStep Particle " <<minTimeStep << " " <<tStep 1702 << " " <<kt->GetDefin 1623 << " " <<kt->GetDefinition()->GetParticleName() 1703 << " 4mom " << kt->Ge 1624 << " 4mom " << kt->GetTrackingMomentum()<<G4endl; 1704 if ( ! intersect ); 1625 if ( ! intersect ); 1705 { 1626 { 1706 PrintKTVector(&theSeconda 1627 PrintKTVector(&theSecondaryList, std::string(" state ERROR.....")); 1707 throw G4HadronicException 1628 throw G4HadronicException(__FILE__, __LINE__, "G4BinaryCascade::StepParticlesOut() particle not in nucleus"); 1708 } 1629 } 1709 #endif 1630 #endif 1710 if(intersect && tStep<minTime 1631 if(intersect && tStep<minTimeStep && tStep> 0 ) 1711 { 1632 { 1712 minTimeStep = tStep; 1633 minTimeStep = tStep; 1713 } 1634 } 1714 } else if ( kt->GetState() != G4K 1635 } else if ( kt->GetState() != G4KineticTrack::outside ){ 1715 PrintKTVector(&theSecondaryLi 1636 PrintKTVector(&theSecondaryList, std::string(" state ERROR.....")); 1716 throw G4HadronicException(__F 1637 throw G4HadronicException(__FILE__, __LINE__, "G4BinaryCascade::StepParticlesOut() particle not in nucleus"); 1717 } 1638 } 1718 } 1639 } 1719 minTimeStep *= 1.2; 1640 minTimeStep *= 1.2; >> 1641 // G4cerr << "CaptureCount = "<<counter<<" "<<nsec<<" "<<minTimeStep<<" "<<1*ns<<G4endl; 1720 G4double timeToCollision=DBL_MAX; 1642 G4double timeToCollision=DBL_MAX; 1721 G4CollisionInitialState * nextCollisi << 1643 G4CollisionInitialState * nextCollision=0; 1722 if(theCollisionMgr->Entries() > 0) 1644 if(theCollisionMgr->Entries() > 0) 1723 { 1645 { 1724 nextCollision = theCollisionMgr-> 1646 nextCollision = theCollisionMgr->GetNextCollision(); 1725 timeToCollision = nextCollision-> 1647 timeToCollision = nextCollision->GetCollisionTime()-theCurrentTime; 1726 // G4cout << " NextCollision * << 1648 G4cout << " NextCollision * , Time= " << nextCollision << " " >> 1649 <<timeToCollision<< G4endl; 1727 } 1650 } 1728 if ( timeToCollision > minTimeStep ) 1651 if ( timeToCollision > minTimeStep ) 1729 { 1652 { 1730 DoTimeStep(minTimeStep); 1653 DoTimeStep(minTimeStep); 1731 ++counter; 1654 ++counter; 1732 } else 1655 } else 1733 { 1656 { 1734 if (!DoTimeStep(timeToCollision) 1657 if (!DoTimeStep(timeToCollision) ) 1735 { 1658 { 1736 // Check if nextCollision is 1659 // Check if nextCollision is still valid, ie. partcile did not leave nucleus 1737 if (theCollisionMgr->GetNextC 1660 if (theCollisionMgr->GetNextCollision() != nextCollision ) 1738 { 1661 { 1739 nextCollision = nullptr; << 1662 nextCollision = 0; 1740 } 1663 } 1741 } 1664 } 1742 // G4cerr <<"post- DoTimeStep 3"< 1665 // G4cerr <<"post- DoTimeStep 3"<<G4endl; 1743 1666 1744 if(nextCollision) 1667 if(nextCollision) 1745 { 1668 { 1746 if ( ApplyCollision(nextColl 1669 if ( ApplyCollision(nextCollision)) 1747 { 1670 { 1748 // G4cout << "ApplyCollis 1671 // G4cout << "ApplyCollision sucess " << G4endl; 1749 } else 1672 } else 1750 { 1673 { 1751 theCollisionMgr->RemoveCo 1674 theCollisionMgr->RemoveCollision(nextCollision); 1752 } 1675 } 1753 } 1676 } 1754 } 1677 } 1755 1678 1756 if(countreset>100) 1679 if(countreset>100) 1757 { 1680 { 1758 #ifdef debug_G4BinaryCascade 1681 #ifdef debug_G4BinaryCascade 1759 G4cerr << "G4BinaryCascade.cc: Wa 1682 G4cerr << "G4BinaryCascade.cc: Warning - aborting looping particle(s)" << G4endl; 1760 PrintKTVector(&theSecondaryList," 1683 PrintKTVector(&theSecondaryList," looping particles added to theFinalState"); 1761 #endif 1684 #endif 1762 1685 1763 // add left secondaries to Final 1686 // add left secondaries to FinalSate 1764 for (auto iter=theSecondaryList.c << 1687 std::vector<G4KineticTrack *>::iterator iter; >> 1688 for ( iter =theSecondaryList.begin(); iter != theSecondaryList.end(); ++iter) 1765 { 1689 { 1766 theFinalState.push_back(*iter 1690 theFinalState.push_back(*iter); 1767 } 1691 } 1768 theSecondaryList.clear(); 1692 theSecondaryList.clear(); 1769 1693 1770 break; 1694 break; 1771 } 1695 } 1772 1696 1773 if(Absorb()) 1697 if(Absorb()) 1774 { 1698 { 1775 // haveProducts = true; 1699 // haveProducts = true; 1776 // G4cout << "Absorb sucess " << 1700 // G4cout << "Absorb sucess " << G4endl; 1777 } 1701 } 1778 1702 1779 if(Capture(false)) 1703 if(Capture(false)) 1780 { 1704 { 1781 // haveProducts = true; 1705 // haveProducts = true; 1782 #ifdef debug_BIC_StepParticlesOut 1706 #ifdef debug_BIC_StepParticlesOut 1783 G4cout << "Capture sucess " << G4 1707 G4cout << "Capture sucess " << G4endl; 1784 #endif 1708 #endif 1785 } 1709 } 1786 if ( counter > 100 && theCollisionMgr << 1710 if ( counter > 100 && theCollisionMgr->Entries() == 0) // no collision, and stepping a while.... 1787 { 1711 { 1788 #ifdef debug_BIC_StepParticlesOut 1712 #ifdef debug_BIC_StepParticlesOut 1789 PrintKTVector(&theSecondaryList,s 1713 PrintKTVector(&theSecondaryList,std::string("stepping 100 steps")); 1790 #endif 1714 #endif 1791 FindCollisions(&theSecondaryList) 1715 FindCollisions(&theSecondaryList); 1792 counter=0; 1716 counter=0; 1793 ++countreset; 1717 ++countreset; 1794 } 1718 } 1795 //G4cout << "currentZ @ end loop " << << 1796 if ( ! currentZ ){ << 1797 // nucleus completely destroyed, << 1798 // products = FillVoidNucleusProdu << 1799 #ifdef debug_BIC_return << 1800 G4cout << "return @ Z=0 after col << 1801 PrintKTVector(&theSecondaryList,s << 1802 G4cout << "theTargetList size: " << 1803 PrintKTVector(&theTargetList,std: << 1804 PrintKTVector(&theCapturedList,st << 1805 << 1806 G4cout << " ExcitE be4 Correct : << 1807 G4cout << " Mom Transfered to nuc << 1808 PrintKTVector(&theFinalState,std: << 1809 // G4cout << "returned products: " << 1810 #endif << 1811 } << 1812 << 1813 } 1719 } 1814 // G4cerr <<"Finished capture loop "<<G4 1720 // G4cerr <<"Finished capture loop "<<G4endl; 1815 1721 1816 //G4cerr <<"pre- DoTimeStep 4"<<G4endl; 1722 //G4cerr <<"pre- DoTimeStep 4"<<G4endl; 1817 DoTimeStep(DBL_MAX); 1723 DoTimeStep(DBL_MAX); 1818 //G4cerr <<"post- DoTimeStep 4"<<G4endl; 1724 //G4cerr <<"post- DoTimeStep 4"<<G4endl; >> 1725 >> 1726 1819 } 1727 } 1820 1728 1821 //------------------------------------------- 1729 //---------------------------------------------------------------------------- 1822 G4double G4BinaryCascade::CorrectShortlivedPr 1730 G4double G4BinaryCascade::CorrectShortlivedPrimaryForFermi( 1823 G4KineticTrack* primary,G4KineticTrac 1731 G4KineticTrack* primary,G4KineticTrackVector target_collection) 1824 //------------------------------------------- 1732 //---------------------------------------------------------------------------- 1825 { 1733 { 1826 G4double Efermi(0); 1734 G4double Efermi(0); 1827 if (primary->GetState() == G4KineticTrack 1735 if (primary->GetState() == G4KineticTrack::inside ) { 1828 G4int PDGcode=primary->GetDefinition( 1736 G4int PDGcode=primary->GetDefinition()->GetPDGEncoding(); 1829 Efermi=((G4RKPropagation *)thePropaga 1737 Efermi=((G4RKPropagation *)thePropagator)->GetField(PDGcode,primary->GetPosition()); 1830 1738 1831 if ( std::abs(PDGcode) > 1000 && PDGc << 1739 if ( std::abs(PDGcode > 1000) && PDGcode != 2112 && PDGcode != 2212 ) 1832 { 1740 { 1833 Efermi = ((G4RKPropagation *)theP 1741 Efermi = ((G4RKPropagation *)thePropagator)->GetField(G4Neutron::Neutron()->GetPDGEncoding(),primary->GetPosition()); 1834 G4LorentzVector mom4Primary=prima 1742 G4LorentzVector mom4Primary=primary->Get4Momentum(); 1835 primary->Update4Momentum(mom4Prim 1743 primary->Update4Momentum(mom4Primary.e() - Efermi); 1836 } 1744 } 1837 1745 1838 for (auto titer=target_collection.cbe << 1746 std::vector<G4KineticTrack *>::iterator titer; >> 1747 for ( titer=target_collection.begin() ; titer!=target_collection.end(); ++titer) 1839 { 1748 { 1840 const G4ParticleDefinition * aDef << 1749 G4ParticleDefinition * aDef=(*titer)->GetDefinition(); 1841 G4int aCode=aDef->GetPDGEncoding( 1750 G4int aCode=aDef->GetPDGEncoding(); 1842 G4ThreeVector aPos=(*titer)->GetP 1751 G4ThreeVector aPos=(*titer)->GetPosition(); 1843 Efermi+= ((G4RKPropagation *)theP 1752 Efermi+= ((G4RKPropagation *)thePropagator)->GetField(aCode, aPos); 1844 } 1753 } 1845 } 1754 } 1846 return Efermi; 1755 return Efermi; 1847 } 1756 } 1848 1757 1849 //------------------------------------------- 1758 //---------------------------------------------------------------------------- 1850 G4bool G4BinaryCascade::CorrectShortlivedFina 1759 G4bool G4BinaryCascade::CorrectShortlivedFinalsForFermi(G4KineticTrackVector * products, 1851 G4double initial_Efermi) 1760 G4double initial_Efermi) 1852 //------------------------------------------- 1761 //---------------------------------------------------------------------------- 1853 { 1762 { 1854 G4double final_Efermi(0); 1763 G4double final_Efermi(0); 1855 G4KineticTrackVector resonances; 1764 G4KineticTrackVector resonances; 1856 for (auto i =products->cbegin(); i != pro << 1765 for ( std::vector<G4KineticTrack *>::iterator i =products->begin(); i != products->end(); i++) 1857 { 1766 { 1858 G4int PDGcode=(*i)->GetDefinition()-> 1767 G4int PDGcode=(*i)->GetDefinition()->GetPDGEncoding(); 1859 // G4cout << " PDGcode, state " 1768 // G4cout << " PDGcode, state " << PDGcode << " " << (*i)->GetState()<<G4endl; 1860 final_Efermi+=((G4RKPropagation *)the 1769 final_Efermi+=((G4RKPropagation *)thePropagator)->GetField(PDGcode,(*i)->GetPosition()); 1861 if ( std::abs(PDGcode) > 1000 && PDGc 1770 if ( std::abs(PDGcode) > 1000 && PDGcode != 2112 && PDGcode != 2212 ) 1862 { 1771 { 1863 resonances.push_back(*i); 1772 resonances.push_back(*i); 1864 } 1773 } 1865 } 1774 } 1866 if ( resonances.size() > 0 ) 1775 if ( resonances.size() > 0 ) 1867 { 1776 { 1868 G4double delta_Fermi= (initial_Efermi 1777 G4double delta_Fermi= (initial_Efermi-final_Efermi)/resonances.size(); 1869 for (auto res=resonances.cbegin(); re << 1778 for (std::vector<G4KineticTrack *>::iterator res=resonances.begin(); res != resonances.end(); res++) 1870 { 1779 { 1871 G4LorentzVector mom=(*res)->Get4M 1780 G4LorentzVector mom=(*res)->Get4Momentum(); 1872 G4double mass2=mom.mag2(); 1781 G4double mass2=mom.mag2(); 1873 G4double newEnergy=mom.e() + delt 1782 G4double newEnergy=mom.e() + delta_Fermi; 1874 G4double newEnergy2= newEnergy*ne 1783 G4double newEnergy2= newEnergy*newEnergy; 1875 //G4cout << "mom = " << mom <<" n 1784 //G4cout << "mom = " << mom <<" newE " << newEnergy<< G4endl; 1876 if ( newEnergy2 < mass2 ) 1785 if ( newEnergy2 < mass2 ) 1877 { 1786 { 1878 return false; 1787 return false; 1879 } 1788 } >> 1789 // G4cout << " correct resonance from /to " << mom.e() << " / " << newEnergy<< G4endl; 1880 G4ThreeVector mom3=std::sqrt(newE 1790 G4ThreeVector mom3=std::sqrt(newEnergy2 - mass2) * mom.vect().unit(); 1881 (*res)->Set4Momentum(G4LorentzVec 1791 (*res)->Set4Momentum(G4LorentzVector(mom3,newEnergy)); 1882 //G4cout << " correct resonan << 1883 // " 3mom from/to " << << 1884 } 1792 } 1885 } 1793 } 1886 return true; 1794 return true; 1887 } 1795 } 1888 1796 1889 //------------------------------------------- 1797 //---------------------------------------------------------------------------- 1890 void G4BinaryCascade::CorrectFinalPandE() 1798 void G4BinaryCascade::CorrectFinalPandE() 1891 //------------------------------------------- 1799 //---------------------------------------------------------------------------- 1892 // 1800 // 1893 // Modify momenta of outgoing particles. 1801 // Modify momenta of outgoing particles. 1894 // Assume two body decay, nucleus(@nominal 1802 // Assume two body decay, nucleus(@nominal mass) + sum of final state particles(SFSP). 1895 // momentum of SFSP shall be less than mome 1803 // momentum of SFSP shall be less than momentum for two body decay. 1896 // 1804 // 1897 { 1805 { 1898 #ifdef debug_BIC_CorrectFinalPandE 1806 #ifdef debug_BIC_CorrectFinalPandE 1899 G4cerr << "BIC: -CorrectFinalPandE called 1807 G4cerr << "BIC: -CorrectFinalPandE called" << G4endl; 1900 #endif 1808 #endif 1901 1809 1902 if ( theFinalState.size() == 0 ) return; 1810 if ( theFinalState.size() == 0 ) return; 1903 1811 1904 G4KineticTrackVector::const_iterator i; << 1812 G4KineticTrackVector::iterator i; 1905 G4LorentzVector pNucleus=GetFinal4Momentu 1813 G4LorentzVector pNucleus=GetFinal4Momentum(); 1906 if ( pNucleus.e() == 0 ) return; // ch 1814 if ( pNucleus.e() == 0 ) return; // check against explicit 0 from GetNucleus4Momentum() 1907 #ifdef debug_BIC_CorrectFinalPandE 1815 #ifdef debug_BIC_CorrectFinalPandE 1908 G4cerr << " -CorrectFinalPandE 3" << G4en 1816 G4cerr << " -CorrectFinalPandE 3" << G4endl; 1909 #endif 1817 #endif 1910 G4LorentzVector pFinals(0); 1818 G4LorentzVector pFinals(0); 1911 for(i = theFinalState.cbegin(); i != theF << 1819 G4int nFinals(0); >> 1820 for(i = theFinalState.begin(); i != theFinalState.end(); ++i) 1912 { 1821 { 1913 pFinals += (*i)->Get4Momentum(); 1822 pFinals += (*i)->Get4Momentum(); >> 1823 ++nFinals; 1914 #ifdef debug_BIC_CorrectFinalPandE 1824 #ifdef debug_BIC_CorrectFinalPandE 1915 G4cout <<"CorrectFinalPandE a final " 1825 G4cout <<"CorrectFinalPandE a final " << (*i)->GetDefinition()->GetParticleName() 1916 << " 4mom " << (*i)->G 1826 << " 4mom " << (*i)->Get4Momentum()<< G4endl; 1917 #endif 1827 #endif 1918 } 1828 } 1919 #ifdef debug_BIC_CorrectFinalPandE 1829 #ifdef debug_BIC_CorrectFinalPandE 1920 G4cout << "CorrectFinalPandE pN pF: " <<p 1830 G4cout << "CorrectFinalPandE pN pF: " <<pNucleus << " " <<pFinals << G4endl; 1921 #endif 1831 #endif 1922 G4LorentzVector pCM=pNucleus + pFinals; 1832 G4LorentzVector pCM=pNucleus + pFinals; 1923 1833 1924 G4LorentzRotation toCMS(-pCM.boostVector( 1834 G4LorentzRotation toCMS(-pCM.boostVector()); 1925 pFinals *=toCMS; 1835 pFinals *=toCMS; 1926 #ifdef debug_BIC_CorrectFinalPandE 1836 #ifdef debug_BIC_CorrectFinalPandE 1927 G4cout << "CorrectFinalPandE pCM, CMS pCM 1837 G4cout << "CorrectFinalPandE pCM, CMS pCM " << pCM << " " <<toCMS*pCM<< G4endl; 1928 G4cout << "CorrectFinal CMS pN pF " <<toC 1838 G4cout << "CorrectFinal CMS pN pF " <<toCMS*pNucleus << " " 1929 <<pFinals << G4endl 1839 <<pFinals << G4endl 1930 << " nucleus initial mass : " <<G 1840 << " nucleus initial mass : " <<GetFinal4Momentum().mag() 1931 <<" massInNucleus m(nucleus) m(fi 1841 <<" massInNucleus m(nucleus) m(finals) std::sqrt(s): " << massInNucleus << " " <<pNucleus.mag()<< " " 1932 << pFinals.mag() << " " << pCM.ma 1842 << pFinals.mag() << " " << pCM.mag() << G4endl; 1933 #endif 1843 #endif 1934 1844 1935 G4LorentzRotation toLab = toCMS.inverse() 1845 G4LorentzRotation toLab = toCMS.inverse(); 1936 1846 1937 G4double s0 = pCM.mag2(); 1847 G4double s0 = pCM.mag2(); 1938 G4double m10 = GetIonMass(currentZ,curren 1848 G4double m10 = GetIonMass(currentZ,currentA); 1939 G4double m20 = pFinals.mag(); 1849 G4double m20 = pFinals.mag(); 1940 if( s0-(m10+m20)*(m10+m20) < 0 ) 1850 if( s0-(m10+m20)*(m10+m20) < 0 ) 1941 { 1851 { 1942 #ifdef debug_BIC_CorrectFinalPandE 1852 #ifdef debug_BIC_CorrectFinalPandE 1943 G4cout << "G4BinaryCascade::CorrectFi 1853 G4cout << "G4BinaryCascade::CorrectFinalPandE() : error! " << G4endl; 1944 1854 1945 G4cout << "not enough mass to correct << 1855 G4cout << "not enough mass to correct: mass, A,Z, mass(nucl), mass(finals) " 1946 << (s0-(m10+m20)*(m10+m20)) < << 1856 << std::sqrt(s0-(m10+m20)*(m10+m20)) << " " 1947 << currentA << " " << current 1857 << currentA << " " << currentZ << " " 1948 << m10 << " " << m20 1858 << m10 << " " << m20 1949 << G4endl; 1859 << G4endl; 1950 G4cerr << " -CorrectFinalPandE 4" << 1860 G4cerr << " -CorrectFinalPandE 4" << G4endl; 1951 1861 1952 PrintKTVector(&theFinalState," mass p 1862 PrintKTVector(&theFinalState," mass problem"); 1953 #endif 1863 #endif 1954 return; 1864 return; 1955 } 1865 } 1956 1866 1957 // Three momentum in cm system 1867 // Three momentum in cm system 1958 G4double pInCM = std::sqrt((s0-(m10+m20)* 1868 G4double pInCM = std::sqrt((s0-(m10+m20)*(m10+m20))*(s0-(m10-m20)*(m10-m20))/(4.*s0)); 1959 #ifdef debug_BIC_CorrectFinalPandE 1869 #ifdef debug_BIC_CorrectFinalPandE 1960 G4cout <<" CorrectFinalPandE pInCM new, 1870 G4cout <<" CorrectFinalPandE pInCM new, CURRENT, ratio : " << pInCM 1961 << " " << (pFinals).vect().mag()< 1871 << " " << (pFinals).vect().mag()<< " " << pInCM/(pFinals).vect().mag() << G4endl; 1962 #endif 1872 #endif 1963 if ( pFinals.vect().mag() > pInCM ) 1873 if ( pFinals.vect().mag() > pInCM ) 1964 { 1874 { 1965 G4ThreeVector p3finals=pInCM*pFinals. 1875 G4ThreeVector p3finals=pInCM*pFinals.vect().unit(); 1966 1876 >> 1877 // G4ThreeVector deltap=(p3finals - pFinals.vect() ) / nFinals; 1967 G4double factor=std::max(0.98,pInCM/p 1878 G4double factor=std::max(0.98,pInCM/pFinals.vect().mag()); // small correction 1968 G4LorentzVector qFinals(0); 1879 G4LorentzVector qFinals(0); 1969 for(i = theFinalState.cbegin(); i != << 1880 for(i = theFinalState.begin(); i != theFinalState.end(); ++i) 1970 { 1881 { 1971 // G4ThreeVector p3((toCMS*( 1882 // G4ThreeVector p3((toCMS*(*i)->Get4Momentum()).vect() + deltap); 1972 G4ThreeVector p3(factor*(toCMS*(* 1883 G4ThreeVector p3(factor*(toCMS*(*i)->Get4Momentum()).vect()); 1973 G4LorentzVector p(p3,std::sqrt((* 1884 G4LorentzVector p(p3,std::sqrt((*i)->Get4Momentum().mag2() + p3.mag2())); 1974 qFinals += p; 1885 qFinals += p; 1975 p *= toLab; 1886 p *= toLab; 1976 #ifdef debug_BIC_CorrectFinalPandE 1887 #ifdef debug_BIC_CorrectFinalPandE 1977 G4cout << " final p corrected: " 1888 G4cout << " final p corrected: " << p << G4endl; 1978 #endif 1889 #endif 1979 (*i)->Set4Momentum(p); 1890 (*i)->Set4Momentum(p); 1980 } 1891 } 1981 #ifdef debug_BIC_CorrectFinalPandE 1892 #ifdef debug_BIC_CorrectFinalPandE 1982 G4cout << "CorrectFinalPandE nucleus 1893 G4cout << "CorrectFinalPandE nucleus corrected mass : " << GetFinal4Momentum() << " " 1983 <<GetFinal4Momentum().mag() < 1894 <<GetFinal4Momentum().mag() << G4endl 1984 << " CMS pFinals , mag, 3.mag 1895 << " CMS pFinals , mag, 3.mag : " << qFinals << " " << qFinals.mag() << " " << qFinals.vect().mag()<< G4endl; 1985 G4cerr << " -CorrectFinalPandE 5 " << 1896 G4cerr << " -CorrectFinalPandE 5 " << factor << G4endl; 1986 #endif 1897 #endif 1987 } 1898 } 1988 #ifdef debug_BIC_CorrectFinalPandE 1899 #ifdef debug_BIC_CorrectFinalPandE 1989 else { G4cerr << " -CorrectFinalPandE 6 - 1900 else { G4cerr << " -CorrectFinalPandE 6 - no correction done" << G4endl; } 1990 #endif 1901 #endif 1991 1902 1992 } 1903 } 1993 1904 1994 //------------------------------------------- 1905 //---------------------------------------------------------------------------- 1995 void G4BinaryCascade::UpdateTracksAndCollisio 1906 void G4BinaryCascade::UpdateTracksAndCollisions( 1996 //----------------------------------- 1907 //---------------------------------------------------------------------------- 1997 G4KineticTrackVector * oldSecondaries 1908 G4KineticTrackVector * oldSecondaries, 1998 G4KineticTrackVector * oldTarget, 1909 G4KineticTrackVector * oldTarget, 1999 G4KineticTrackVector * newSecondaries 1910 G4KineticTrackVector * newSecondaries) 2000 { 1911 { 2001 std::vector<G4KineticTrack *>::const_iter << 1912 std::vector<G4KineticTrack *>::iterator iter1, iter2; 2002 1913 2003 // remove old secondaries from the second 1914 // remove old secondaries from the secondary list 2004 if(oldSecondaries) 1915 if(oldSecondaries) 2005 { 1916 { 2006 if(!oldSecondaries->empty()) 1917 if(!oldSecondaries->empty()) 2007 { 1918 { 2008 for(auto iter1=oldSecondaries->cb << 1919 for(iter1 = oldSecondaries->begin(); iter1 != oldSecondaries->end(); >> 1920 ++iter1) 2009 { 1921 { 2010 iter2 = std::find(theSecondar << 1922 iter2 = std::find(theSecondaryList.begin(), theSecondaryList.end(), 2011 if ( iter2 != theSecondaryLis << 1923 *iter1); >> 1924 if ( iter2 != theSecondaryList.end() ) theSecondaryList.erase(iter2); 2012 } 1925 } 2013 theCollisionMgr->RemoveTracksColl 1926 theCollisionMgr->RemoveTracksCollisions(oldSecondaries); 2014 } 1927 } 2015 } 1928 } 2016 1929 2017 // remove old target from the target list 1930 // remove old target from the target list 2018 if(oldTarget) 1931 if(oldTarget) 2019 { 1932 { 2020 // G4cout << "################## Debu 1933 // G4cout << "################## Debugging 0 "<<G4endl; 2021 if(oldTarget->size()!=0) 1934 if(oldTarget->size()!=0) 2022 { 1935 { 2023 1936 2024 // G4cout << "################## 1937 // G4cout << "################## Debugging 1 "<<oldTarget->size()<<G4endl; 2025 for(auto iter1 = oldTarget->cbegi << 1938 for(iter1 = oldTarget->begin(); iter1 != oldTarget->end(); ++iter1) 2026 { 1939 { 2027 iter2 = std::find(theTargetLi << 1940 iter2 = std::find(theTargetList.begin(), theTargetList.end(), >> 1941 *iter1); 2028 theTargetList.erase(iter2); 1942 theTargetList.erase(iter2); 2029 } 1943 } 2030 theCollisionMgr->RemoveTracksColl 1944 theCollisionMgr->RemoveTracksCollisions(oldTarget); 2031 } 1945 } 2032 } 1946 } 2033 1947 2034 if(newSecondaries) 1948 if(newSecondaries) 2035 { 1949 { 2036 if(!newSecondaries->empty()) 1950 if(!newSecondaries->empty()) 2037 { 1951 { 2038 // insert new secondaries in the 1952 // insert new secondaries in the secondary list 2039 for(auto iter1 = newSecondaries-> << 1953 for(iter1 = newSecondaries->begin(); iter1 != newSecondaries->end(); >> 1954 ++iter1) 2040 { 1955 { 2041 theSecondaryList.push_back(*i 1956 theSecondaryList.push_back(*iter1); 2042 if ((*iter1)->GetState() == G 1957 if ((*iter1)->GetState() == G4KineticTrack::undefined) 2043 { 1958 { 2044 PrintKTVector(*iter1, "und 1959 PrintKTVector(*iter1, "undefined in FindCollisions"); 2045 } 1960 } 2046 1961 2047 1962 2048 } 1963 } 2049 // look for collisions of new sec 1964 // look for collisions of new secondaries 2050 FindCollisions(newSecondaries); 1965 FindCollisions(newSecondaries); 2051 } 1966 } 2052 } 1967 } 2053 // G4cout << "Exiting ... "<<oldTarget<<G 1968 // G4cout << "Exiting ... "<<oldTarget<<G4endl; 2054 } 1969 } 2055 1970 2056 1971 2057 class SelectFromKTV 1972 class SelectFromKTV 2058 { 1973 { 2059 private: 1974 private: 2060 G4KineticTrackVector * ktv; 1975 G4KineticTrackVector * ktv; 2061 G4KineticTrack::CascadeState wanted_state 1976 G4KineticTrack::CascadeState wanted_state; 2062 public: 1977 public: 2063 SelectFromKTV(G4KineticTrackVector * out, 1978 SelectFromKTV(G4KineticTrackVector * out, G4KineticTrack::CascadeState astate) 2064 : 1979 : 2065 ktv(out), wanted_state(astate) 1980 ktv(out), wanted_state(astate) 2066 {}; 1981 {}; 2067 void operator() (G4KineticTrack *& kt) co 1982 void operator() (G4KineticTrack *& kt) const 2068 { 1983 { 2069 if ( (kt)->GetState() == wanted_state 1984 if ( (kt)->GetState() == wanted_state ) ktv->push_back(kt); 2070 }; 1985 }; 2071 }; 1986 }; 2072 1987 2073 1988 2074 1989 2075 //------------------------------------------- 1990 //---------------------------------------------------------------------------- 2076 G4bool G4BinaryCascade::DoTimeStep(G4double t 1991 G4bool G4BinaryCascade::DoTimeStep(G4double theTimeStep) 2077 //------------------------------------------- 1992 //---------------------------------------------------------------------------- 2078 { 1993 { 2079 1994 2080 #ifdef debug_BIC_DoTimeStep 1995 #ifdef debug_BIC_DoTimeStep 2081 G4ping debug("debug_G4BinaryCascade"); 1996 G4ping debug("debug_G4BinaryCascade"); 2082 debug.push_back("======> DoTimeStep 1"); 1997 debug.push_back("======> DoTimeStep 1"); debug.dump(); 2083 G4cerr <<"G4BinaryCascade::DoTimeStep: en 1998 G4cerr <<"G4BinaryCascade::DoTimeStep: enter step="<< theTimeStep 2084 << " , time="<<theCurrentTime << 1999 << " , time="<<theCurrentTime << G4endl; 2085 PrintKTVector(&theSecondaryList, std::str 2000 PrintKTVector(&theSecondaryList, std::string("DoTimeStep - theSecondaryList")); 2086 //PrintKTVector(&theTargetList, std::stri 2001 //PrintKTVector(&theTargetList, std::string("DoTimeStep - theTargetList")); 2087 #endif 2002 #endif 2088 2003 2089 G4bool success=true; 2004 G4bool success=true; 2090 std::vector<G4KineticTrack *>::const_iter << 2005 std::vector<G4KineticTrack *>::iterator iter; 2091 2006 2092 G4KineticTrackVector * kt_outside = new G 2007 G4KineticTrackVector * kt_outside = new G4KineticTrackVector; 2093 std::for_each( theSecondaryList.begin(),t 2008 std::for_each( theSecondaryList.begin(),theSecondaryList.end(), 2094 SelectFromKTV(kt_outside,G4Kineti 2009 SelectFromKTV(kt_outside,G4KineticTrack::outside)); 2095 //PrintKTVector(kt_outside, std::string(" 2010 //PrintKTVector(kt_outside, std::string("DoTimeStep - found outside")); 2096 2011 2097 G4KineticTrackVector * kt_inside = new G4 2012 G4KineticTrackVector * kt_inside = new G4KineticTrackVector; 2098 std::for_each( theSecondaryList.begin(),t 2013 std::for_each( theSecondaryList.begin(),theSecondaryList.end(), 2099 SelectFromKTV(kt_inside, G4Kineti 2014 SelectFromKTV(kt_inside, G4KineticTrack::inside)); 2100 // PrintKTVector(kt_inside, std::string( 2015 // PrintKTVector(kt_inside, std::string("DoTimeStep - found inside")); 2101 //----- 2016 //----- 2102 G4KineticTrackVector dummy; // needed f 2017 G4KineticTrackVector dummy; // needed for re-usability 2103 #ifdef debug_BIC_DoTimeStep 2018 #ifdef debug_BIC_DoTimeStep 2104 G4cout << "NOW WE ARE ENTERING THE TRANSP 2019 G4cout << "NOW WE ARE ENTERING THE TRANSPORT"<<G4endl; 2105 #endif 2020 #endif 2106 2021 2107 // =================== Here we move the p 2022 // =================== Here we move the particles =================== 2108 2023 2109 thePropagator->Transport(theSecondaryList 2024 thePropagator->Transport(theSecondaryList, dummy, theTimeStep); 2110 2025 2111 // =================== Here we move the p 2026 // =================== Here we move the particles =================== 2112 2027 2113 //------ 2028 //------ 2114 2029 2115 theMomentumTransfer += thePropagator->Get 2030 theMomentumTransfer += thePropagator->GetMomentumTransfer(); 2116 #ifdef debug_BIC_DoTimeStep 2031 #ifdef debug_BIC_DoTimeStep 2117 G4cout << "DoTimeStep : theMomentumTransf 2032 G4cout << "DoTimeStep : theMomentumTransfer = " << theMomentumTransfer << G4endl; 2118 PrintKTVector(&theSecondaryList, std::str 2033 PrintKTVector(&theSecondaryList, std::string("DoTimeStep - secondaries aft trsprt")); 2119 #endif 2034 #endif 2120 2035 2121 //_DebugEpConservation(" after stepping") << 2122 << 2123 // Partclies which went INTO nucleus 2036 // Partclies which went INTO nucleus 2124 2037 2125 G4KineticTrackVector * kt_gone_in = new G 2038 G4KineticTrackVector * kt_gone_in = new G4KineticTrackVector; 2126 std::for_each( kt_outside->begin(),kt_out 2039 std::for_each( kt_outside->begin(),kt_outside->end(), 2127 SelectFromKTV(kt_gone_in,G4Kineti 2040 SelectFromKTV(kt_gone_in,G4KineticTrack::inside)); 2128 // PrintKTVector(kt_gone_in, std::string 2041 // PrintKTVector(kt_gone_in, std::string("DoTimeStep - gone in")); 2129 2042 2130 2043 2131 // Partclies which went OUT OF nucleus 2044 // Partclies which went OUT OF nucleus 2132 G4KineticTrackVector * kt_gone_out = new 2045 G4KineticTrackVector * kt_gone_out = new G4KineticTrackVector; 2133 std::for_each( kt_inside->begin(),kt_insi 2046 std::for_each( kt_inside->begin(),kt_inside->end(), 2134 SelectFromKTV(kt_gone_out, G4Kine 2047 SelectFromKTV(kt_gone_out, G4KineticTrack::gone_out)); 2135 2048 2136 // PrintKTVector(kt_gone_out, std::strin 2049 // PrintKTVector(kt_gone_out, std::string("DoTimeStep - gone out")); 2137 2050 2138 G4KineticTrackVector *fail=CorrectBarions 2051 G4KineticTrackVector *fail=CorrectBarionsOnBoundary(kt_gone_in,kt_gone_out); 2139 2052 2140 if ( fail ) 2053 if ( fail ) 2141 { 2054 { 2142 // some particle(s) supposed to enter 2055 // some particle(s) supposed to enter/leave were miss_nucleus/captured by the correction 2143 kt_gone_in->clear(); 2056 kt_gone_in->clear(); 2144 std::for_each( kt_outside->begin(),kt 2057 std::for_each( kt_outside->begin(),kt_outside->end(), 2145 SelectFromKTV(kt_gone_in,G4Ki 2058 SelectFromKTV(kt_gone_in,G4KineticTrack::inside)); 2146 2059 2147 kt_gone_out->clear(); 2060 kt_gone_out->clear(); 2148 std::for_each( kt_inside->begin(),kt_ 2061 std::for_each( kt_inside->begin(),kt_inside->end(), 2149 SelectFromKTV(kt_gone_out, G4 2062 SelectFromKTV(kt_gone_out, G4KineticTrack::gone_out)); 2150 2063 2151 #ifdef debug_BIC_DoTimeStep 2064 #ifdef debug_BIC_DoTimeStep 2152 PrintKTVector(fail,std::string(" Fail 2065 PrintKTVector(fail,std::string(" Failed to go in/out -> miss_nucleus/captured")); 2153 PrintKTVector(kt_gone_in, std::string 2066 PrintKTVector(kt_gone_in, std::string("recreated kt_gone_in")); 2154 PrintKTVector(kt_gone_out, std::strin 2067 PrintKTVector(kt_gone_out, std::string("recreated kt_gone_out")); 2155 #endif 2068 #endif 2156 delete fail; 2069 delete fail; 2157 } 2070 } 2158 2071 2159 // Add tracks missing nucleus and tracks 2072 // Add tracks missing nucleus and tracks going straight though to addFinals 2160 std::for_each( kt_outside->begin(),kt_out 2073 std::for_each( kt_outside->begin(),kt_outside->end(), 2161 SelectFromKTV(kt_gone_out,G4Kinet 2074 SelectFromKTV(kt_gone_out,G4KineticTrack::miss_nucleus)); 2162 //PrintKTVector(kt_gone_out, std::string( 2075 //PrintKTVector(kt_gone_out, std::string("miss to append to final state..")); 2163 std::for_each( kt_outside->begin(),kt_out 2076 std::for_each( kt_outside->begin(),kt_outside->end(), 2164 SelectFromKTV(kt_gone_out,G4Kinet 2077 SelectFromKTV(kt_gone_out,G4KineticTrack::gone_out)); 2165 2078 2166 #ifdef debug_BIC_DoTimeStep 2079 #ifdef debug_BIC_DoTimeStep 2167 PrintKTVector(kt_gone_out, std::string("a 2080 PrintKTVector(kt_gone_out, std::string("append gone_outs to final state.. theFinalState")); 2168 #endif 2081 #endif 2169 2082 2170 theFinalState.insert(theFinalState.end(), 2083 theFinalState.insert(theFinalState.end(), 2171 kt_gone_out->begin(),kt_gone_out- 2084 kt_gone_out->begin(),kt_gone_out->end()); 2172 2085 2173 // Partclies which could not leave nucleu 2086 // Partclies which could not leave nucleus, captured... 2174 G4KineticTrackVector * kt_captured = new 2087 G4KineticTrackVector * kt_captured = new G4KineticTrackVector; 2175 std::for_each( theSecondaryList.begin(),t 2088 std::for_each( theSecondaryList.begin(),theSecondaryList.end(), 2176 SelectFromKTV(kt_captured, G4Kine 2089 SelectFromKTV(kt_captured, G4KineticTrack::captured)); 2177 2090 2178 // Check no track is part in next collisi 2091 // Check no track is part in next collision, ie. 2179 // this step was to far, and collisions 2092 // this step was to far, and collisions should not occur any more 2180 2093 2181 if ( theCollisionMgr->Entries()> 0 ) 2094 if ( theCollisionMgr->Entries()> 0 ) 2182 { 2095 { 2183 if (kt_gone_out->size() ) 2096 if (kt_gone_out->size() ) 2184 { 2097 { 2185 G4KineticTrack * nextPrimary = th 2098 G4KineticTrack * nextPrimary = theCollisionMgr->GetNextCollision()->GetPrimary(); 2186 iter = std::find(kt_gone_out->beg 2099 iter = std::find(kt_gone_out->begin(),kt_gone_out->end(),nextPrimary); 2187 if ( iter != kt_gone_out->cend() << 2100 if ( iter != kt_gone_out->end() ) 2188 { 2101 { 2189 success=false; 2102 success=false; 2190 #ifdef debug_BIC_DoTimeStep 2103 #ifdef debug_BIC_DoTimeStep 2191 G4cout << " DoTimeStep - WARN 2104 G4cout << " DoTimeStep - WARNING: deleting current collision!" << G4endl; 2192 #endif 2105 #endif 2193 } 2106 } 2194 } 2107 } 2195 if ( kt_captured->size() ) 2108 if ( kt_captured->size() ) 2196 { 2109 { 2197 G4KineticTrack * nextPrimary = th 2110 G4KineticTrack * nextPrimary = theCollisionMgr->GetNextCollision()->GetPrimary(); 2198 iter = std::find(kt_captured->beg 2111 iter = std::find(kt_captured->begin(),kt_captured->end(),nextPrimary); 2199 if ( iter != kt_captured->cend() << 2112 if ( iter != kt_captured->end() ) 2200 { 2113 { 2201 success=false; 2114 success=false; 2202 #ifdef debug_BIC_DoTimeStep 2115 #ifdef debug_BIC_DoTimeStep 2203 G4cout << " DoTimeStep - WARN 2116 G4cout << " DoTimeStep - WARNING: deleting current collision!" << G4endl; 2204 #endif 2117 #endif 2205 } 2118 } 2206 } 2119 } 2207 2120 2208 } 2121 } 2209 // PrintKTVector(kt_gone_out," kt_gone_ou 2122 // PrintKTVector(kt_gone_out," kt_gone_out be4 updatetrack..."); 2210 UpdateTracksAndCollisions(kt_gone_out,0 , 2123 UpdateTracksAndCollisions(kt_gone_out,0 ,0); 2211 2124 2212 2125 2213 if ( kt_captured->size() ) 2126 if ( kt_captured->size() ) 2214 { 2127 { 2215 theCapturedList.insert(theCapturedLis 2128 theCapturedList.insert(theCapturedList.end(), 2216 kt_captured->begin(),kt_captu 2129 kt_captured->begin(),kt_captured->end()); 2217 //should be std::for_each(kt_cap 2130 //should be std::for_each(kt_captured->begin(),kt_captured->end(), 2218 // std::mem_fun(&G4Kinet 2131 // std::mem_fun(&G4KineticTrack::Hit)); 2219 // but VC 6 requires: 2132 // but VC 6 requires: 2220 for(auto i_captured=kt_captured->cbeg << 2133 std::vector<G4KineticTrack *>::iterator i_captured; >> 2134 for(i_captured=kt_captured->begin();i_captured!=kt_captured->end();i_captured++) 2221 { 2135 { 2222 (*i_captured)->Hit(); 2136 (*i_captured)->Hit(); 2223 } 2137 } 2224 // PrintKTVector(kt_captured," kt 2138 // PrintKTVector(kt_captured," kt_captured be4 updatetrack..."); 2225 UpdateTracksAndCollisions(kt_captured << 2139 UpdateTracksAndCollisions(kt_captured, NULL, NULL); 2226 } 2140 } 2227 2141 2228 #ifdef debug_G4BinaryCascade 2142 #ifdef debug_G4BinaryCascade 2229 delete kt_inside; 2143 delete kt_inside; 2230 kt_inside = new G4KineticTrackVector; 2144 kt_inside = new G4KineticTrackVector; 2231 std::for_each( theSecondaryList.begin(),t 2145 std::for_each( theSecondaryList.begin(),theSecondaryList.end(), 2232 SelectFromKTV(kt_inside, G4Kineti 2146 SelectFromKTV(kt_inside, G4KineticTrack::inside)); 2233 if ( currentZ != (GetTotalCharge(theTarge 2147 if ( currentZ != (GetTotalCharge(theTargetList) 2234 + GetTotalCharge(theCapturedList) 2148 + GetTotalCharge(theCapturedList) 2235 + GetTotalCharge(*kt_inside)) ) 2149 + GetTotalCharge(*kt_inside)) ) 2236 { 2150 { 2237 G4cout << " error-DoTimeStep aft, A, 2151 G4cout << " error-DoTimeStep aft, A, Z: " << currentA << " " << currentZ 2238 << " sum(tgt,capt,active) " 2152 << " sum(tgt,capt,active) " 2239 << GetTotalCharge(theTargetLi 2153 << GetTotalCharge(theTargetList) + GetTotalCharge(theCapturedList) + GetTotalCharge(*kt_inside) 2240 << " targets: " << GetTotalC 2154 << " targets: " << GetTotalCharge(theTargetList) 2241 << " captured: " << GetTotalC 2155 << " captured: " << GetTotalCharge(theCapturedList) 2242 << " active: " << GetTotalC 2156 << " active: " << GetTotalCharge(*kt_inside) 2243 << G4endl; 2157 << G4endl; 2244 } 2158 } 2245 #endif 2159 #endif 2246 2160 2247 delete kt_inside; 2161 delete kt_inside; 2248 delete kt_outside; 2162 delete kt_outside; 2249 delete kt_captured; 2163 delete kt_captured; 2250 delete kt_gone_in; 2164 delete kt_gone_in; 2251 delete kt_gone_out; 2165 delete kt_gone_out; 2252 2166 2253 // G4cerr <<"G4BinaryCascade::DoTimeStep 2167 // G4cerr <<"G4BinaryCascade::DoTimeStep: exit "<<G4endl; 2254 theCurrentTime += theTimeStep; 2168 theCurrentTime += theTimeStep; 2255 2169 2256 //debug.push_back("======> DoTimeStep 2") 2170 //debug.push_back("======> DoTimeStep 2"); debug.dump(); 2257 return success; 2171 return success; 2258 2172 2259 } 2173 } 2260 2174 2261 //------------------------------------------- 2175 //---------------------------------------------------------------------------- 2262 G4KineticTrackVector* G4BinaryCascade::Correc 2176 G4KineticTrackVector* G4BinaryCascade::CorrectBarionsOnBoundary( 2263 G4KineticTrackVector *in, 2177 G4KineticTrackVector *in, 2264 G4KineticTrackVector *out) 2178 G4KineticTrackVector *out) 2265 //------------------------------------------- 2179 //---------------------------------------------------------------------------- 2266 { 2180 { 2267 G4KineticTrackVector * kt_fail(nullptr); << 2181 G4KineticTrackVector * kt_fail(0); 2268 std::vector<G4KineticTrack *>::const_iter << 2182 std::vector<G4KineticTrack *>::iterator iter; 2269 // G4cout << "CorrectBarionsOnBoundary,c 2183 // G4cout << "CorrectBarionsOnBoundary,currentZ,currentA," 2270 // << currentZ << " "<< currentA 2184 // << currentZ << " "<< currentA << G4endl; 2271 if (in->size()) 2185 if (in->size()) 2272 { 2186 { 2273 G4int secondaries_in(0); 2187 G4int secondaries_in(0); 2274 G4int secondaryBarions_in(0); 2188 G4int secondaryBarions_in(0); 2275 G4int secondaryCharge_in(0); 2189 G4int secondaryCharge_in(0); 2276 G4double secondaryMass_in(0); 2190 G4double secondaryMass_in(0); 2277 2191 2278 for ( iter =in->cbegin(); iter != in- << 2192 for ( iter =in->begin(); iter != in->end(); ++iter) 2279 { 2193 { 2280 ++secondaries_in; 2194 ++secondaries_in; 2281 secondaryCharge_in += G4lrint((*i 2195 secondaryCharge_in += G4lrint((*iter)->GetDefinition()->GetPDGCharge()/eplus); 2282 if ((*iter)->GetDefinition()->Get 2196 if ((*iter)->GetDefinition()->GetBaryonNumber()!=0 ) 2283 { 2197 { 2284 secondaryBarions_in += (*iter 2198 secondaryBarions_in += (*iter)->GetDefinition()->GetBaryonNumber(); 2285 if((*iter)->GetDefinition() = 2199 if((*iter)->GetDefinition() == G4Neutron::Neutron() || 2286 (*iter)->GetDefinitio 2200 (*iter)->GetDefinition() == G4Proton::Proton() ) 2287 { 2201 { 2288 secondaryMass_in += (*ite 2202 secondaryMass_in += (*iter)->GetDefinition()->GetPDGMass(); 2289 } else { 2203 } else { 2290 secondaryMass_in += G4Pro 2204 secondaryMass_in += G4Proton::Proton()->GetPDGMass(); 2291 } 2205 } 2292 } 2206 } 2293 } 2207 } 2294 G4double mass_initial= GetIonMass(cur 2208 G4double mass_initial= GetIonMass(currentZ,currentA); 2295 2209 2296 currentZ += secondaryCharge_in; 2210 currentZ += secondaryCharge_in; 2297 currentA += secondaryBarions_in; 2211 currentA += secondaryBarions_in; 2298 2212 2299 // G4cout << "CorrectBarionsOnBounda 2213 // G4cout << "CorrectBarionsOnBoundary,secondaryCharge_in, secondaryBarions_in " 2300 // << secondaryCharge_in < 2214 // << secondaryCharge_in << " "<< secondaryBarions_in << G4endl; 2301 2215 2302 G4double mass_final= GetIonMass(curre 2216 G4double mass_final= GetIonMass(currentZ,currentA); 2303 2217 2304 G4double correction= secondaryMass_in 2218 G4double correction= secondaryMass_in + mass_initial - mass_final; 2305 if (secondaries_in>1) 2219 if (secondaries_in>1) 2306 {correction /= secondaries_in;} 2220 {correction /= secondaries_in;} 2307 2221 2308 #ifdef debug_BIC_CorrectBarionsOnBoundary 2222 #ifdef debug_BIC_CorrectBarionsOnBoundary 2309 G4cout << "CorrectBarionsOnBoundary,c 2223 G4cout << "CorrectBarionsOnBoundary,currentZ,currentA," 2310 << "secondaryCharge_in,second 2224 << "secondaryCharge_in,secondaryBarions_in," 2311 << "energy correction,m_secon 2225 << "energy correction,m_secondry,m_nucl_init,m_nucl_final " 2312 << currentZ << " "<< currentA 2226 << currentZ << " "<< currentA <<" " 2313 << secondaryCharge_in<<" "<<s 2227 << secondaryCharge_in<<" "<<secondaryBarions_in<<" " 2314 << correction << " " 2228 << correction << " " 2315 << secondaryMass_in << " " 2229 << secondaryMass_in << " " 2316 << mass_initial << " " 2230 << mass_initial << " " 2317 << mass_final << " " 2231 << mass_final << " " 2318 << G4endl; 2232 << G4endl; 2319 PrintKTVector(in,std::string("in be4 2233 PrintKTVector(in,std::string("in be4 correction")); 2320 #endif 2234 #endif 2321 for ( iter = in->cbegin(); iter != in << 2235 >> 2236 for ( iter = in->begin(); iter != in->end(); ++iter) 2322 { 2237 { 2323 if ((*iter)->GetTrackingMomentum( 2238 if ((*iter)->GetTrackingMomentum().e()+correction > (*iter)->GetActualMass()) 2324 { 2239 { 2325 (*iter)->UpdateTrackingMoment 2240 (*iter)->UpdateTrackingMomentum((*iter)->GetTrackingMomentum().e() + correction); 2326 } else { 2241 } else { 2327 //particle cannot go in, put 2242 //particle cannot go in, put to miss_nucleus 2328 G4RKPropagation * RKprop=(G4R 2243 G4RKPropagation * RKprop=(G4RKPropagation *)thePropagator; 2329 (*iter)->SetState(G4KineticTr 2244 (*iter)->SetState(G4KineticTrack::miss_nucleus); 2330 // Undo correction for Colomb 2245 // Undo correction for Colomb Barrier 2331 G4double barrier=RKprop->GetB 2246 G4double barrier=RKprop->GetBarrier((*iter)->GetDefinition()->GetPDGEncoding()); 2332 (*iter)->UpdateTrackingMoment 2247 (*iter)->UpdateTrackingMomentum((*iter)->GetTrackingMomentum().e() + barrier); 2333 if ( ! kt_fail ) kt_fail=new 2248 if ( ! kt_fail ) kt_fail=new G4KineticTrackVector; 2334 kt_fail->push_back(*iter); 2249 kt_fail->push_back(*iter); 2335 currentZ -= G4lrint((*iter)-> 2250 currentZ -= G4lrint((*iter)->GetDefinition()->GetPDGCharge()/eplus); 2336 currentA -= (*iter)->GetDefin 2251 currentA -= (*iter)->GetDefinition()->GetBaryonNumber(); 2337 2252 2338 } 2253 } 2339 2254 2340 } 2255 } 2341 #ifdef debug_BIC_CorrectBarionsOnBoundary 2256 #ifdef debug_BIC_CorrectBarionsOnBoundary 2342 G4cout << " CorrectBarionsOnBoundary, << 2257 G4cout << " CorrectBarionsOnBoundary, aft, A, Z, sec-Z,A,m,m_in_nucleus " 2343 << currentZ << " " << current << 2258 << currentA << " " << currentZ << " " 2344 << secondaryCharge_in << " " 2259 << secondaryCharge_in << " " << secondaryBarions_in << " " 2345 << secondaryMass_in << " " 2260 << secondaryMass_in << " " 2346 << G4endl; 2261 << G4endl; 2347 PrintKTVector(in,std::string("in AFT 2262 PrintKTVector(in,std::string("in AFT correction")); 2348 #endif 2263 #endif 2349 2264 2350 } 2265 } 2351 //--------------------------------------- 2266 //---------------------------------------------- 2352 if (out->size()) 2267 if (out->size()) 2353 { 2268 { 2354 G4int secondaries_out(0); 2269 G4int secondaries_out(0); 2355 G4int secondaryBarions_out(0); 2270 G4int secondaryBarions_out(0); 2356 G4int secondaryCharge_out(0); 2271 G4int secondaryCharge_out(0); 2357 G4double secondaryMass_out(0); 2272 G4double secondaryMass_out(0); 2358 2273 2359 for ( iter = out->cbegin(); iter != o << 2274 for ( iter =out->begin(); iter != out->end(); ++iter) 2360 { 2275 { 2361 ++secondaries_out; 2276 ++secondaries_out; 2362 secondaryCharge_out += G4lrint((* 2277 secondaryCharge_out += G4lrint((*iter)->GetDefinition()->GetPDGCharge()/eplus); 2363 if ((*iter)->GetDefinition()->Get 2278 if ((*iter)->GetDefinition()->GetBaryonNumber() !=0 ) 2364 { 2279 { 2365 secondaryBarions_out += (*ite 2280 secondaryBarions_out += (*iter)->GetDefinition()->GetBaryonNumber(); 2366 if((*iter)->GetDefinition() = 2281 if((*iter)->GetDefinition() == G4Neutron::Neutron() || 2367 (*iter)->GetDefinitio 2282 (*iter)->GetDefinition() == G4Proton::Proton() ) 2368 { 2283 { 2369 secondaryMass_out += (*it 2284 secondaryMass_out += (*iter)->GetDefinition()->GetPDGMass(); 2370 } else { 2285 } else { 2371 secondaryMass_out += G4Ne 2286 secondaryMass_out += G4Neutron::Neutron()->GetPDGMass(); 2372 } 2287 } 2373 } 2288 } 2374 } 2289 } 2375 2290 2376 G4double mass_initial= GetIonMass(cu 2291 G4double mass_initial= GetIonMass(currentZ,currentA); 2377 currentA -=secondaryBarions_out; 2292 currentA -=secondaryBarions_out; 2378 currentZ -=secondaryCharge_out; 2293 currentZ -=secondaryCharge_out; 2379 2294 2380 // G4cout << "CorrectBarionsOnBounda << 2295 // G4cout << "CorrectBarionsOnBoundary,secondaryCharge_out, secondaryBarions_out" 2381 // << secondaryCharge_out 2296 // << secondaryCharge_out << " "<< secondaryBarions_out << G4endl; 2382 2297 2383 // a delta min 2298 // a delta minus will do currentZ < 0 in light nuclei 2384 // if (currentA < 0 || currentZ < 2299 // if (currentA < 0 || currentZ < 0 ) 2385 if (currentA < 0 ) 2300 if (currentA < 0 ) 2386 { 2301 { 2387 G4cerr << "G4BinaryCascade - seco 2302 G4cerr << "G4BinaryCascade - secondaryBarions_out,secondaryCharge_out " << 2388 secondaryBarions_out << " 2303 secondaryBarions_out << " " << secondaryCharge_out << G4endl; 2389 PrintKTVector(&theTargetList,"Cor 2304 PrintKTVector(&theTargetList,"CorrectBarionsOnBoundary Target"); 2390 PrintKTVector(&theCapturedList,"C 2305 PrintKTVector(&theCapturedList,"CorrectBarionsOnBoundary Captured"); 2391 PrintKTVector(&theSecondaryList," 2306 PrintKTVector(&theSecondaryList,"CorrectBarionsOnBoundary Secondaries"); 2392 G4cerr << "G4BinaryCascade - curr 2307 G4cerr << "G4BinaryCascade - currentA, currentZ " << currentA << " " << currentZ << G4endl; 2393 throw G4HadronicException(__FILE_ 2308 throw G4HadronicException(__FILE__, __LINE__, "G4BinaryCascade::CorrectBarionsOnBoundary() - fatal error"); 2394 } 2309 } 2395 G4double mass_final=GetIonMass(curren 2310 G4double mass_final=GetIonMass(currentZ,currentA); 2396 G4double correction= mass_initial - m 2311 G4double correction= mass_initial - mass_final - secondaryMass_out; 2397 // G4cout << "G4BinaryCascade::Correc << 2398 2312 2399 if (secondaries_out>1) correction /= 2313 if (secondaries_out>1) correction /= secondaries_out; 2400 #ifdef debug_BIC_CorrectBarionsOnBoundary 2314 #ifdef debug_BIC_CorrectBarionsOnBoundary 2401 G4cout << "DoTimeStep,(current Z,A)," << 2315 G4cout << "DoTimeStep,currentZ,currentA," 2402 << "(secondaries out,Charge,B << 2316 << "secondaries_out," 2403 <<"* energy correction,(m_sec << 2317 <<"secondaryCharge_out,secondaryBarions_out," 2404 << "("<< currentZ << ","<< cu << 2318 <<"energy correction,m_secondry,m_nucl_init,m_nucl_final " 2405 << secondaries_out << "," << 2319 << " "<< currentZ << " "<< currentA <<" " 2406 << secondaryCharge_out<<","<< << 2320 << secondaries_out << " " 2407 << correction << " (" << 2321 << secondaryCharge_out<<" "<<secondaryBarions_out<<" " 2408 << secondaryMass_out << ", " << 2322 << correction << " " 2409 << mass_initial << ", " << 2323 << secondaryMass_out << " " 2410 << mass_final << ")" << 2324 << mass_initial << " " >> 2325 << mass_final << " " 2411 << G4endl; 2326 << G4endl; 2412 PrintKTVector(out,std::string("out be 2327 PrintKTVector(out,std::string("out be4 correction")); 2413 #endif 2328 #endif 2414 2329 2415 for ( iter = out->cbegin(); iter != o << 2330 for ( iter = out->begin(); iter != out->end(); ++iter) 2416 { 2331 { 2417 if ((*iter)->GetTrackingMomentum( 2332 if ((*iter)->GetTrackingMomentum().e()+correction > (*iter)->GetActualMass()) 2418 { 2333 { 2419 (*iter)->UpdateTrackingMoment 2334 (*iter)->UpdateTrackingMomentum((*iter)->GetTrackingMomentum().e() + correction); 2420 } else 2335 } else 2421 { 2336 { 2422 // particle cannot go out due 2337 // particle cannot go out due to change of nuclear potential! 2423 // capture protons and neutr 2338 // capture protons and neutrons; 2424 if(((*iter)->GetDefinition() 2339 if(((*iter)->GetDefinition() == G4Proton::Proton()) || 2425 ((*iter)->GetDefiniti 2340 ((*iter)->GetDefinition() == G4Neutron::Neutron())) 2426 { 2341 { 2427 (*iter)->SetState(G4Kinet 2342 (*iter)->SetState(G4KineticTrack::captured); 2428 // Undo correction for Co 2343 // Undo correction for Colomb Barrier 2429 G4double barrier=((G4RKPr 2344 G4double barrier=((G4RKPropagation *)thePropagator)->GetBarrier((*iter)->GetDefinition()->GetPDGEncoding()); 2430 (*iter)->UpdateTrackingMo 2345 (*iter)->UpdateTrackingMomentum((*iter)->GetTrackingMomentum().e() - barrier); 2431 if ( kt_fail == 0 ) kt_fa 2346 if ( kt_fail == 0 ) kt_fail=new G4KineticTrackVector; 2432 kt_fail->push_back(*iter) 2347 kt_fail->push_back(*iter); 2433 currentZ += G4lrint((*ite 2348 currentZ += G4lrint((*iter)->GetDefinition()->GetPDGCharge()/eplus); 2434 currentA += (*iter)->GetD 2349 currentA += (*iter)->GetDefinition()->GetBaryonNumber(); 2435 } 2350 } 2436 #ifdef debug_BIC_CorrectBarionsOnBoundary 2351 #ifdef debug_BIC_CorrectBarionsOnBoundary 2437 else 2352 else 2438 { 2353 { 2439 G4cout << "Not correcting 2354 G4cout << "Not correcting outgoing " << *iter << " " 2440 << (*iter)->GetDe 2355 << (*iter)->GetDefinition()->GetPDGEncoding() << " " 2441 << (*iter)->GetDe 2356 << (*iter)->GetDefinition()->GetParticleName() << G4endl; 2442 PrintKTVector(out,std::st 2357 PrintKTVector(out,std::string("outgoing, one not corrected")); 2443 } 2358 } 2444 #endif 2359 #endif 2445 } 2360 } 2446 } 2361 } 2447 2362 2448 #ifdef debug_BIC_CorrectBarionsOnBoundary 2363 #ifdef debug_BIC_CorrectBarionsOnBoundary 2449 PrintKTVector(out,std::string("out AF 2364 PrintKTVector(out,std::string("out AFTER correction")); 2450 G4cout << " DoTimeStep, nucl-update, 2365 G4cout << " DoTimeStep, nucl-update, A, Z, sec-Z,A,m,m_in_nucleus, table-mass, delta " 2451 << currentA << " "<< currentZ 2366 << currentA << " "<< currentZ << " " 2452 << secondaryCharge_out << " " 2367 << secondaryCharge_out << " "<< secondaryBarions_out << " "<< 2453 secondaryMass_out << " " 2368 secondaryMass_out << " " 2454 << massInNucleus << " " 2369 << massInNucleus << " " 2455 << GetIonMass(currentZ,curren << 2370 << G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(currentZ,currentA) 2456 << " " << massInNucleus - Get << 2371 << " " << massInNucleus -G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(currentZ,currentA) 2457 << G4endl; 2372 << G4endl; 2458 #endif 2373 #endif 2459 } 2374 } 2460 2375 2461 return kt_fail; 2376 return kt_fail; 2462 } 2377 } 2463 2378 2464 2379 2465 //------------------------------------------- 2380 //---------------------------------------------------------------------------- 2466 2381 2467 G4Fragment * G4BinaryCascade::FindFragments() 2382 G4Fragment * G4BinaryCascade::FindFragments() 2468 //------------------------------------------- 2383 //---------------------------------------------------------------------------- 2469 { 2384 { 2470 2385 2471 #ifdef debug_BIC_FindFragments 2386 #ifdef debug_BIC_FindFragments 2472 G4cout << "target, captured, secondary: " 2387 G4cout << "target, captured, secondary: " 2473 << theTargetList.size() << " " 2388 << theTargetList.size() << " " 2474 << theCapturedList.size()<< " " 2389 << theCapturedList.size()<< " " 2475 << theSecondaryList.size() 2390 << theSecondaryList.size() 2476 << G4endl; 2391 << G4endl; 2477 #endif 2392 #endif 2478 2393 2479 G4int a = G4int(theTargetList.size()+theC << 2394 G4int a = theTargetList.size()+theCapturedList.size(); 2480 G4int zTarget = 0; 2395 G4int zTarget = 0; 2481 for(auto i = theTargetList.cbegin(); i != << 2396 G4KineticTrackVector::iterator i; >> 2397 for(i = theTargetList.begin(); i != theTargetList.end(); ++i) 2482 { 2398 { 2483 if(G4lrint((*i)->GetDefinition()->Get 2399 if(G4lrint((*i)->GetDefinition()->GetPDGCharge()/eplus) == 1 ) 2484 { 2400 { 2485 zTarget++; 2401 zTarget++; 2486 } 2402 } 2487 } 2403 } 2488 2404 2489 G4int zCaptured = 0; 2405 G4int zCaptured = 0; 2490 G4LorentzVector CapturedMomentum(0.,0.,0. 2406 G4LorentzVector CapturedMomentum(0.,0.,0.,0.); 2491 for(auto i = theCapturedList.cbegin(); i << 2407 for(i = theCapturedList.begin(); i != theCapturedList.end(); ++i) 2492 { 2408 { 2493 CapturedMomentum += (*i)->Get4Momentu 2409 CapturedMomentum += (*i)->Get4Momentum(); 2494 if(G4lrint((*i)->GetDefinition()->Get 2410 if(G4lrint((*i)->GetDefinition()->GetPDGCharge()/eplus) == 1 ) 2495 { 2411 { 2496 zCaptured++; 2412 zCaptured++; 2497 } 2413 } 2498 } 2414 } 2499 2415 2500 G4int z = zTarget+zCaptured; 2416 G4int z = zTarget+zCaptured; 2501 2417 2502 #ifdef debug_G4BinaryCascade 2418 #ifdef debug_G4BinaryCascade 2503 if ( z != (GetTotalCharge(theTargetList) 2419 if ( z != (GetTotalCharge(theTargetList) + GetTotalCharge(theCapturedList)) ) 2504 { 2420 { 2505 G4cout << " FindFragment Counting err 2421 G4cout << " FindFragment Counting error z a " << z << " " <<a << " " 2506 << GetTotalCharge(theTargetLi 2422 << GetTotalCharge(theTargetList) << " " << GetTotalCharge(theCapturedList)<< 2507 G4endl; 2423 G4endl; 2508 PrintKTVector(&theTargetList, std::st 2424 PrintKTVector(&theTargetList, std::string("theTargetList")); 2509 PrintKTVector(&theCapturedList, std:: 2425 PrintKTVector(&theCapturedList, std::string("theCapturedList")); 2510 } 2426 } 2511 #endif 2427 #endif 2512 //debug 2428 //debug 2513 /* 2429 /* 2514 * G4cout << " Fragment mass table / re 2430 * G4cout << " Fragment mass table / real " 2515 * << GetIonMass(z, a) << 2431 * << G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(z, a) 2516 * << " / " << GetFinal4Momentum().mag( 2432 * << " / " << GetFinal4Momentum().mag() 2517 * << " difference " 2433 * << " difference " 2518 * << GetFinal4Momentum().mag() - GetI << 2434 * << GetFinal4Momentum().mag() - G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(z, a) 2519 * << G4endl; 2435 * << G4endl; 2520 */ 2436 */ 2521 // 2437 // 2522 // if(fBCDEBUG) G4cerr << "Fragment A, Z << 2438 // if(getenv("BCDEBUG") ) G4cerr << "Fragment A, Z "<< a <<" "<< z<<G4endl; 2523 if ( z < 1 ) return 0; 2439 if ( z < 1 ) return 0; 2524 2440 2525 G4int holes = G4int(the3DNucleus->GetMass << 2441 G4int holes = the3DNucleus->GetMassNumber() - theTargetList.size(); 2526 G4int excitons = (G4int)theCapturedList.s << 2442 G4int excitons = theCapturedList.size(); 2527 #ifdef debug_BIC_FindFragments 2443 #ifdef debug_BIC_FindFragments 2528 G4cout << "Fragment: a= " << a << " z= " 2444 G4cout << "Fragment: a= " << a << " z= " << z << " particles= " << excitons 2529 << " Charged= " << zCaptured << " 2445 << " Charged= " << zCaptured << " holes= " << holes 2530 << " excitE= " <<GetExcitationEne 2446 << " excitE= " <<GetExcitationEnergy() 2531 << " Final4Momentum= " << GetFina 2447 << " Final4Momentum= " << GetFinalNucleusMomentum() << " capturMomentum= " << CapturedMomentum 2532 << G4endl; 2448 << G4endl; 2533 #endif 2449 #endif 2534 2450 2535 G4Fragment * fragment = new G4Fragment(a, 2451 G4Fragment * fragment = new G4Fragment(a,z,GetFinalNucleusMomentum()); 2536 fragment->SetNumberOfHoles(holes); 2452 fragment->SetNumberOfHoles(holes); 2537 2453 2538 //GF fragment->SetNumberOfParticles(exci 2454 //GF fragment->SetNumberOfParticles(excitons-holes); 2539 fragment->SetNumberOfParticles(excitons); 2455 fragment->SetNumberOfParticles(excitons); 2540 fragment->SetNumberOfCharged(zCaptured); 2456 fragment->SetNumberOfCharged(zCaptured); 2541 fragment->SetCreatorModelID(theBIC_ID); << 2457 G4ParticleDefinition * aIonDefinition = >> 2458 G4ParticleTable::GetParticleTable()->FindIon(a,z,0,z); >> 2459 fragment->SetParticleDefinition(aIonDefinition); 2542 2460 2543 return fragment; 2461 return fragment; 2544 } 2462 } 2545 2463 2546 //------------------------------------------- 2464 //---------------------------------------------------------------------------- 2547 2465 2548 G4LorentzVector G4BinaryCascade::GetFinal4Mom 2466 G4LorentzVector G4BinaryCascade::GetFinal4Momentum() 2549 //------------------------------------------- 2467 //---------------------------------------------------------------------------- 2550 // Return momentum of reminder nulceus; 2468 // Return momentum of reminder nulceus; 2551 // ie. difference of (initial state(primarie 2469 // ie. difference of (initial state(primaries+nucleus) - final state) particles, ignoring remnant nucleus 2552 { 2470 { 2553 G4LorentzVector final4Momentum = theIniti 2471 G4LorentzVector final4Momentum = theInitial4Mom + theProjectile4Momentum; 2554 G4LorentzVector finals(0,0,0,0); 2472 G4LorentzVector finals(0,0,0,0); 2555 for(auto i = theFinalState.cbegin(); i != << 2473 for(G4KineticTrackVector::iterator i = theFinalState.begin(); i != theFinalState.end(); ++i) 2556 { 2474 { 2557 final4Momentum -= (*i)->Get4Momentum( 2475 final4Momentum -= (*i)->Get4Momentum(); 2558 finals += (*i)->Get4Momentum(); 2476 finals += (*i)->Get4Momentum(); 2559 } 2477 } 2560 2478 2561 if(final4Momentum.e()> 0 && (final4Moment << 2479 if((final4Momentum.vect()/final4Momentum.e()).mag()>1.0 && currentA > 0) 2562 { 2480 { 2563 #ifdef debug_BIC_Final4Momentum 2481 #ifdef debug_BIC_Final4Momentum 2564 G4cerr << G4endl; 2482 G4cerr << G4endl; 2565 G4cerr << "G4BinaryCascade::GetFinal4 2483 G4cerr << "G4BinaryCascade::GetFinal4Momentum - Fatal"<<G4endl; 2566 G4KineticTrackVector::iterator i; 2484 G4KineticTrackVector::iterator i; 2567 G4cerr <<"Total initial 4-momentum " 2485 G4cerr <<"Total initial 4-momentum " << theProjectile4Momentum << G4endl; 2568 G4cerr <<" GetFinal4Momentum: Initial 2486 G4cerr <<" GetFinal4Momentum: Initial nucleus "<<theInitial4Mom<<G4endl; 2569 for(i = theFinalState.begin(); i != t 2487 for(i = theFinalState.begin(); i != theFinalState.end(); ++i) 2570 { 2488 { 2571 G4cerr <<" Final state: "<<(*i)-> 2489 G4cerr <<" Final state: "<<(*i)->Get4Momentum()<<(*i)->GetDefinition()->GetParticleName()<<G4endl; 2572 } 2490 } 2573 G4cerr << "Sum( 4-mom ) finals " << f << 2491 G4cerr << "Sum( 4-mon ) finals " << finals << G4endl; 2574 G4cerr<< " Final4Momentum = "<<final4 2492 G4cerr<< " Final4Momentum = "<<final4Momentum <<" "<<final4Momentum.m()<<G4endl; 2575 G4cerr <<" current A, Z = "<< current 2493 G4cerr <<" current A, Z = "<< currentA<<", "<<currentZ<<G4endl; 2576 G4cerr << G4endl; 2494 G4cerr << G4endl; 2577 #endif 2495 #endif 2578 2496 2579 final4Momentum=G4LorentzVector(0,0,0, 2497 final4Momentum=G4LorentzVector(0,0,0,0); 2580 } 2498 } 2581 return final4Momentum; 2499 return final4Momentum; 2582 } 2500 } 2583 2501 2584 //------------------------------------------- 2502 //---------------------------------------------------------------------------- 2585 G4LorentzVector G4BinaryCascade::GetFinalNucl 2503 G4LorentzVector G4BinaryCascade::GetFinalNucleusMomentum() 2586 //------------------------------------------- 2504 //---------------------------------------------------------------------------- 2587 { 2505 { 2588 // return momentum of nucleus for use wit 2506 // return momentum of nucleus for use with precompound model; also keep transformation to 2589 // apply to precompoud products. 2507 // apply to precompoud products. 2590 2508 2591 G4LorentzVector CapturedMomentum(0,0,0,0) 2509 G4LorentzVector CapturedMomentum(0,0,0,0); >> 2510 G4KineticTrackVector::iterator i; 2592 // G4cout << "GetFinalNucleusMomentum Ca 2511 // G4cout << "GetFinalNucleusMomentum Captured size: " <<theCapturedList.size() << G4endl; 2593 for(auto i = theCapturedList.cbegin(); i << 2512 for(i = theCapturedList.begin(); i != theCapturedList.end(); ++i) 2594 { 2513 { 2595 CapturedMomentum += (*i)->Get4Momentu 2514 CapturedMomentum += (*i)->Get4Momentum(); 2596 } 2515 } 2597 //G4cout << "GetFinalNucleusMomentum Capt 2516 //G4cout << "GetFinalNucleusMomentum CapturedMomentum= " <<CapturedMomentum << G4endl; 2598 // G4cerr << "it 9"<<G4endl; 2517 // G4cerr << "it 9"<<G4endl; 2599 2518 2600 G4LorentzVector NucleusMomentum = GetFina 2519 G4LorentzVector NucleusMomentum = GetFinal4Momentum(); 2601 if ( NucleusMomentum.e() > 0 ) 2520 if ( NucleusMomentum.e() > 0 ) 2602 { 2521 { 2603 // G4cout << "GetFinalNucleusMomentum 2522 // G4cout << "GetFinalNucleusMomentum GetFinal4Momentum= " <<NucleusMomentum <<" "<<NucleusMomentum.mag()<<G4endl; 2604 // boost nucleus to a frame such that 2523 // boost nucleus to a frame such that the momentum of nucleus == momentum of Captured 2605 G4ThreeVector boost= (NucleusMomentum 2524 G4ThreeVector boost= (NucleusMomentum.vect() -CapturedMomentum.vect())/NucleusMomentum.e(); 2606 if(boost.mag2()>1.0) 2525 if(boost.mag2()>1.0) 2607 { 2526 { 2608 # ifdef debug_BIC_FinalNucleusMomentum 2527 # ifdef debug_BIC_FinalNucleusMomentum 2609 G4cerr << "G4BinaryCascade::GetFi 2528 G4cerr << "G4BinaryCascade::GetFinalNucleusMomentum - Fatal"<<G4endl; 2610 G4cerr << "it 0"<<boost <<G4endl; 2529 G4cerr << "it 0"<<boost <<G4endl; 2611 G4cerr << "it 01"<<NucleusMomentu 2530 G4cerr << "it 01"<<NucleusMomentum<<" "<<CapturedMomentum<<" "<<G4endl; 2612 G4cout <<" testing boost "<<boost 2531 G4cout <<" testing boost "<<boost<<" "<<boost.mag()<<G4endl; 2613 # endif 2532 # endif 2614 boost=G4ThreeVector(0,0,0); 2533 boost=G4ThreeVector(0,0,0); 2615 NucleusMomentum=G4LorentzVector(0 2534 NucleusMomentum=G4LorentzVector(0,0,0,0); 2616 } 2535 } 2617 G4LorentzRotation nucleusBoost( -boo 2536 G4LorentzRotation nucleusBoost( -boost ); 2618 precompoundLorentzboost.set( boost ); 2537 precompoundLorentzboost.set( boost ); 2619 #ifdef debug_debug_BIC_FinalNucleusMomentum 2538 #ifdef debug_debug_BIC_FinalNucleusMomentum 2620 G4cout << "GetFinalNucleusMomentum be 2539 G4cout << "GetFinalNucleusMomentum be4 boostNucleusMomentum, CapturedMomentum"<<NucleusMomentum<<" "<<CapturedMomentum<<" "<<G4endl; 2621 #endif 2540 #endif 2622 NucleusMomentum *= nucleusBoost; 2541 NucleusMomentum *= nucleusBoost; 2623 #ifdef debug_BIC_FinalNucleusMomentum 2542 #ifdef debug_BIC_FinalNucleusMomentum 2624 G4cout << "GetFinalNucleusMomentum af 2543 G4cout << "GetFinalNucleusMomentum aft boost GetFinal4Momentum= " <<NucleusMomentum <<G4endl; 2625 #endif 2544 #endif 2626 } 2545 } 2627 return NucleusMomentum; 2546 return NucleusMomentum; 2628 } 2547 } 2629 2548 2630 //------------------------------------------- 2549 //---------------------------------------------------------------------------- 2631 G4ReactionProductVector * G4BinaryCascade::Pr 2550 G4ReactionProductVector * G4BinaryCascade::Propagate1H1( 2632 //----------------------------------- 2551 //---------------------------------------------------------------------------- 2633 G4KineticTrackVector * secondaries, G 2552 G4KineticTrackVector * secondaries, G4V3DNucleus * nucleus) 2634 { 2553 { 2635 G4ReactionProductVector * products = new 2554 G4ReactionProductVector * products = new G4ReactionProductVector; 2636 const G4ParticleDefinition * aHTarg = G4P << 2555 G4ParticleDefinition * aHTarg = G4Proton::ProtonDefinition(); 2637 if (nucleus->GetCharge() == 0) aHTarg = G << 2638 G4double mass = aHTarg->GetPDGMass(); 2556 G4double mass = aHTarg->GetPDGMass(); 2639 G4KineticTrackVector * secs = nullptr; << 2557 if (nucleus->GetCharge() == 0) aHTarg = G4Neutron::NeutronDefinition(); >> 2558 mass = aHTarg->GetPDGMass(); >> 2559 G4KineticTrackVector * secs = 0; 2640 G4ThreeVector pos(0,0,0); 2560 G4ThreeVector pos(0,0,0); 2641 G4LorentzVector mom(mass); 2561 G4LorentzVector mom(mass); 2642 G4KineticTrack aTarget(aHTarg, 0., pos, m 2562 G4KineticTrack aTarget(aHTarg, 0., pos, mom); 2643 G4bool done(false); 2563 G4bool done(false); >> 2564 std::vector<G4KineticTrack *>::iterator iter, jter; 2644 // data member static G4Scatterer theH 2565 // data member static G4Scatterer theH1Scatterer; 2645 //G4cout << " start 1H1 for " << (*second 2566 //G4cout << " start 1H1 for " << (*secondaries).front()->GetDefinition()->GetParticleName() 2646 // << " on " << aHTarg->GetParticle 2567 // << " on " << aHTarg->GetParticleName() << G4endl; 2647 G4int tryCount(0); 2568 G4int tryCount(0); 2648 while(!done && tryCount++ <200) << 2569 while(!done && tryCount++ <200) 2649 { 2570 { 2650 if(secs) 2571 if(secs) 2651 { 2572 { 2652 std::for_each(secs->begin(), secs 2573 std::for_each(secs->begin(), secs->end(), DeleteKineticTrack()); 2653 delete secs; 2574 delete secs; 2654 } 2575 } 2655 secs = theH1Scatterer->Scatter(*(*sec 2576 secs = theH1Scatterer->Scatter(*(*secondaries).front(), aTarget); 2656 #ifdef debug_H1_BinaryCascade 2577 #ifdef debug_H1_BinaryCascade 2657 PrintKTVector(secs," From Scatter"); 2578 PrintKTVector(secs," From Scatter"); 2658 #endif 2579 #endif 2659 for(std::size_t ss=0; secs && ss<secs << 2580 for(size_t ss=0; secs && ss<secs->size(); ss++) 2660 { 2581 { 2661 // must have one resonance in fin 2582 // must have one resonance in final state, or it was elastic, not allowed here. 2662 if((*secs)[ss]->GetDefinition()-> 2583 if((*secs)[ss]->GetDefinition()->IsShortLived()) done = true; 2663 } 2584 } 2664 } 2585 } 2665 2586 >> 2587 size_t current(0); 2666 ClearAndDestroy(&theFinalState); 2588 ClearAndDestroy(&theFinalState); 2667 ClearAndDestroy(secondaries); << 2589 for(current=0; secs && current<secs->size(); current++) 2668 delete secondaries; << 2669 << 2670 for(std::size_t current=0; secs && curren << 2671 { 2590 { 2672 if((*secs)[current]->GetDefinition()- 2591 if((*secs)[current]->GetDefinition()->IsShortLived()) 2673 { 2592 { 2674 done = true; // must have one re 2593 done = true; // must have one resonance in final state, elastic not allowed here! 2675 G4KineticTrackVector * dec = (*se 2594 G4KineticTrackVector * dec = (*secs)[current]->Decay(); 2676 for(auto jter=dec->cbegin(); jter << 2595 for(jter=dec->begin(); jter != dec->end(); jter++) 2677 { 2596 { 2678 //G4cout << "Decay"<<G4endl; 2597 //G4cout << "Decay"<<G4endl; 2679 secs->push_back(*jter); 2598 secs->push_back(*jter); 2680 //G4cout << "decay "<<(*jter) 2599 //G4cout << "decay "<<(*jter)->GetDefinition()->GetParticleName()<<G4endl; 2681 } 2600 } 2682 delete (*secs)[current]; 2601 delete (*secs)[current]; 2683 delete dec; 2602 delete dec; 2684 } 2603 } 2685 else 2604 else 2686 { 2605 { 2687 //G4cout << "FS "<<G4endl; 2606 //G4cout << "FS "<<G4endl; 2688 //G4cout << "FS "<<(*secs)[curren 2607 //G4cout << "FS "<<(*secs)[current]->GetDefinition()->GetParticleName()<<G4endl; 2689 theFinalState.push_back((*secs)[c 2608 theFinalState.push_back((*secs)[current]); 2690 } 2609 } 2691 } 2610 } 2692 2611 2693 delete secs; 2612 delete secs; 2694 #ifdef debug_H1_BinaryCascade 2613 #ifdef debug_H1_BinaryCascade 2695 PrintKTVector(&theFinalState," FinalState 2614 PrintKTVector(&theFinalState," FinalState"); 2696 #endif 2615 #endif 2697 for(auto iter = theFinalState.cbegin(); i << 2616 for(iter = theFinalState.begin(); iter != theFinalState.end(); ++iter) 2698 { 2617 { 2699 G4KineticTrack * kt = *iter; 2618 G4KineticTrack * kt = *iter; 2700 G4ReactionProduct * aNew = new G4Reac 2619 G4ReactionProduct * aNew = new G4ReactionProduct(kt->GetDefinition()); 2701 aNew->SetMomentum(kt->Get4Momentum(). 2620 aNew->SetMomentum(kt->Get4Momentum().vect()); 2702 aNew->SetTotalEnergy(kt->Get4Momentum 2621 aNew->SetTotalEnergy(kt->Get4Momentum().e()); 2703 aNew->SetCreatorModelID(theBIC_ID); << 2704 aNew->SetParentResonanceDef(kt->GetPa << 2705 aNew->SetParentResonanceID(kt->GetPar << 2706 products->push_back(aNew); 2622 products->push_back(aNew); 2707 #ifdef debug_H1_BinaryCascade 2623 #ifdef debug_H1_BinaryCascade 2708 if (! kt->GetDefinition()->GetPDGStab 2624 if (! kt->GetDefinition()->GetPDGStable() ) 2709 { 2625 { 2710 if (kt->GetDefinition()->IsShortL 2626 if (kt->GetDefinition()->IsShortLived()) 2711 { 2627 { 2712 G4cout << "final shortlived : 2628 G4cout << "final shortlived : "; 2713 } else 2629 } else 2714 { 2630 { 2715 G4cout << "final un stable : 2631 G4cout << "final un stable : "; 2716 } 2632 } 2717 G4cout <<kt->GetDefinition()->Get 2633 G4cout <<kt->GetDefinition()->GetParticleName()<< G4endl; 2718 } 2634 } 2719 #endif 2635 #endif 2720 delete kt; 2636 delete kt; 2721 } 2637 } 2722 theFinalState.clear(); 2638 theFinalState.clear(); 2723 return products; 2639 return products; 2724 2640 2725 } 2641 } 2726 2642 2727 //------------------------------------------- 2643 //---------------------------------------------------------------------------- 2728 G4ThreeVector G4BinaryCascade::GetSpherePoint 2644 G4ThreeVector G4BinaryCascade::GetSpherePoint( 2729 G4double r, const G4LorentzVector & m 2645 G4double r, const G4LorentzVector & mom4) 2730 //------------------------------------------- 2646 //---------------------------------------------------------------------------- 2731 { 2647 { 2732 // Get a point outside radius. 2648 // Get a point outside radius. 2733 // point is random in plane (circle o 2649 // point is random in plane (circle of radius r) orthogonal to mom, 2734 // plus -1*r*mom->vect()->unit(); 2650 // plus -1*r*mom->vect()->unit(); 2735 G4ThreeVector o1, o2; 2651 G4ThreeVector o1, o2; 2736 G4ThreeVector mom = mom4.vect(); 2652 G4ThreeVector mom = mom4.vect(); 2737 2653 2738 o1= mom.orthogonal(); // we simply need 2654 o1= mom.orthogonal(); // we simply need any vector non parallel 2739 o2= mom.cross(o1); // o2 is now orth 2655 o2= mom.cross(o1); // o2 is now orthogonal to mom and o1, ie. o1 and o2 define plane. 2740 2656 2741 G4double x2, x1; 2657 G4double x2, x1; 2742 2658 2743 do 2659 do 2744 { 2660 { 2745 x1=(G4UniformRand()-.5)*2; 2661 x1=(G4UniformRand()-.5)*2; 2746 x2=(G4UniformRand()-.5)*2; 2662 x2=(G4UniformRand()-.5)*2; 2747 } while (sqr(x1) +sqr(x2) > 1.); << 2663 } while (sqr(x1) +sqr(x2) > 1.); 2748 2664 2749 return G4ThreeVector(r*(x1*o1.unit() + x2 2665 return G4ThreeVector(r*(x1*o1.unit() + x2*o2.unit() - 1.5* mom.unit())); 2750 2666 2751 2667 2752 2668 2753 /* 2669 /* 2754 * // Get a point uniformly distributed o 2670 * // Get a point uniformly distributed on the surface of a sphere, 2755 * // with z < 0. 2671 * // with z < 0. 2756 * G4double b = r*G4UniformRand(); // 2672 * G4double b = r*G4UniformRand(); // impact parameter 2757 * G4double phi = G4UniformRand()*2*pi; 2673 * G4double phi = G4UniformRand()*2*pi; 2758 * G4double x = b*std::cos(phi); 2674 * G4double x = b*std::cos(phi); 2759 * G4double y = b*std::sin(phi); 2675 * G4double y = b*std::sin(phi); 2760 * G4double z = -std::sqrt(r*r-b*b); 2676 * G4double z = -std::sqrt(r*r-b*b); 2761 * z *= 1.001; // Get position a little 2677 * z *= 1.001; // Get position a little bit out of the sphere... 2762 * point.setX(x); 2678 * point.setX(x); 2763 * point.setY(y); 2679 * point.setY(y); 2764 * point.setZ(z); 2680 * point.setZ(z); 2765 */ 2681 */ 2766 } 2682 } 2767 2683 2768 //------------------------------------------- 2684 //---------------------------------------------------------------------------- 2769 2685 2770 void G4BinaryCascade::ClearAndDestroy(G4Kinet 2686 void G4BinaryCascade::ClearAndDestroy(G4KineticTrackVector * ktv) 2771 //------------------------------------------- 2687 //---------------------------------------------------------------------------- 2772 { 2688 { 2773 for(auto i = ktv->cbegin(); i != ktv->cen << 2689 std::vector<G4KineticTrack *>::iterator i; >> 2690 for(i = ktv->begin(); i != ktv->end(); ++i) 2774 delete (*i); 2691 delete (*i); 2775 ktv->clear(); 2692 ktv->clear(); 2776 } 2693 } 2777 2694 2778 //------------------------------------------- 2695 //---------------------------------------------------------------------------- 2779 void G4BinaryCascade::ClearAndDestroy(G4React 2696 void G4BinaryCascade::ClearAndDestroy(G4ReactionProductVector * rpv) 2780 //------------------------------------------- 2697 //---------------------------------------------------------------------------- 2781 { 2698 { 2782 for(auto i = rpv->cbegin(); i != rpv->cen << 2699 std::vector<G4ReactionProduct *>::iterator i; >> 2700 for(i = rpv->begin(); i != rpv->end(); ++i) 2783 delete (*i); 2701 delete (*i); 2784 rpv->clear(); 2702 rpv->clear(); 2785 } 2703 } 2786 2704 2787 //------------------------------------------- 2705 //---------------------------------------------------------------------------- 2788 void G4BinaryCascade::PrintKTVector(G4Kinetic 2706 void G4BinaryCascade::PrintKTVector(G4KineticTrackVector * ktv, std::string comment) 2789 //------------------------------------------- 2707 //---------------------------------------------------------------------------- 2790 { 2708 { 2791 if (comment.size() > 0 ) G4cout << "G4Bin 2709 if (comment.size() > 0 ) G4cout << "G4BinaryCascade::PrintKTVector() " << comment << G4endl; 2792 if (ktv) { 2710 if (ktv) { 2793 G4cout << " vector: " << ktv << ", n 2711 G4cout << " vector: " << ktv << ", number of tracks: " << ktv->size() 2794 << G4endl; 2712 << G4endl; 2795 std::vector<G4KineticTrack *>::const_ << 2713 std::vector<G4KineticTrack *>::iterator i; 2796 G4int count; 2714 G4int count; 2797 2715 2798 for(count = 0, i = ktv->cbegin(); i ! << 2716 for(count = 0, i = ktv->begin(); i != ktv->end(); ++i, ++count) 2799 { 2717 { 2800 G4KineticTrack * kt = *i; 2718 G4KineticTrack * kt = *i; 2801 G4cout << " track n. " << count; 2719 G4cout << " track n. " << count; 2802 PrintKTVector(kt); 2720 PrintKTVector(kt); 2803 } 2721 } 2804 } else { 2722 } else { 2805 G4cout << "G4BinaryCascade::PrintKTVe 2723 G4cout << "G4BinaryCascade::PrintKTVector():No KineticTrackVector given " << G4endl; 2806 } 2724 } 2807 } 2725 } 2808 //------------------------------------------- 2726 //---------------------------------------------------------------------------- 2809 void G4BinaryCascade::PrintKTVector(G4Kinetic 2727 void G4BinaryCascade::PrintKTVector(G4KineticTrack * kt, std::string comment) 2810 //------------------------------------------- 2728 //---------------------------------------------------------------------------- 2811 { 2729 { 2812 if (comment.size() > 0 ) G4cout << "G4Bin 2730 if (comment.size() > 0 ) G4cout << "G4BinaryCascade::PrintKTVector() "<< comment << G4endl; 2813 if ( kt ){ 2731 if ( kt ){ 2814 G4cout << ", id: " << kt << G4endl; 2732 G4cout << ", id: " << kt << G4endl; 2815 G4ThreeVector pos = kt->GetPosition() 2733 G4ThreeVector pos = kt->GetPosition(); 2816 G4LorentzVector mom = kt->Get4Momentu 2734 G4LorentzVector mom = kt->Get4Momentum(); 2817 G4LorentzVector tmom = kt->GetTrackin 2735 G4LorentzVector tmom = kt->GetTrackingMomentum(); 2818 const G4ParticleDefinition * definiti << 2736 G4ParticleDefinition * definition = kt->GetDefinition(); 2819 G4cout << " definition: " << defin 2737 G4cout << " definition: " << definition->GetPDGEncoding() << " pos: " 2820 << 1/fermi*pos << " R: " << 1 2738 << 1/fermi*pos << " R: " << 1/fermi*pos.mag() << " 4mom: " 2821 << 1/MeV*mom <<"Tr_mom" << 1 2739 << 1/MeV*mom <<"Tr_mom" << 1/MeV*tmom << " P: " << 1/MeV*mom.vect().mag() 2822 << " M: " << 1/MeV*mom.mag() 2740 << " M: " << 1/MeV*mom.mag() << G4endl; 2823 G4cout <<" trackstatus: "<<kt->Get << 2741 G4cout <<" trackstatus: "<<kt->GetState()<<G4endl; 2824 } else { 2742 } else { 2825 G4cout << "G4BinaryCascade::PrintKTVe 2743 G4cout << "G4BinaryCascade::PrintKTVector(): No Kinetictrack given" << G4endl; 2826 } 2744 } 2827 } 2745 } 2828 2746 2829 2747 2830 //------------------------------------------- 2748 //---------------------------------------------------------------------------- 2831 G4double G4BinaryCascade::GetIonMass(G4int Z, 2749 G4double G4BinaryCascade::GetIonMass(G4int Z, G4int A) 2832 //------------------------------------------- 2750 //---------------------------------------------------------------------------- 2833 { 2751 { 2834 G4double mass(0); 2752 G4double mass(0); 2835 if ( Z > 0 && A >= Z ) 2753 if ( Z > 0 && A >= Z ) 2836 { 2754 { 2837 mass = G4ParticleTable::GetParticleTa 2755 mass = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(Z,A); 2838 2756 2839 } else if ( A > 0 && Z>0 ) 2757 } else if ( A > 0 && Z>0 ) 2840 { 2758 { 2841 // charge Z > A; will happen for ligh 2759 // charge Z > A; will happen for light nuclei with pions involved. 2842 mass = G4ParticleTable::GetParticleTa 2760 mass = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(A,A); 2843 2761 2844 } else if ( A >= 0 && Z<=0 ) 2762 } else if ( A >= 0 && Z<=0 ) 2845 { 2763 { 2846 // all neutral, or empty nucleus 2764 // all neutral, or empty nucleus 2847 mass = A * G4Neutron::Neutron()->GetP 2765 mass = A * G4Neutron::Neutron()->GetPDGMass(); 2848 2766 2849 } else if ( A == 0 ) << 2767 } else if ( A == 0 && std::abs(Z)<2 ) 2850 { 2768 { 2851 // empty nucleus, except maybe pions 2769 // empty nucleus, except maybe pions 2852 mass = 0; 2770 mass = 0; 2853 2771 2854 } else 2772 } else 2855 { 2773 { 2856 G4cerr << "G4BinaryCascade::GetIonMas 2774 G4cerr << "G4BinaryCascade::GetIonMass() - invalid (A,Z) = (" 2857 << A << "," << Z << ")" <<G4e 2775 << A << "," << Z << ")" <<G4endl; 2858 throw G4HadronicException(__FILE__, _ 2776 throw G4HadronicException(__FILE__, __LINE__, "G4BinaryCascade::GetIonMass() - giving up"); 2859 2777 2860 } 2778 } 2861 // G4cout << "G4BinaryCascade::GetIonMass << 2862 return mass; 2779 return mass; 2863 } 2780 } 2864 G4ReactionProductVector * G4BinaryCascade::Fi 2781 G4ReactionProductVector * G4BinaryCascade::FillVoidNucleusProducts(G4ReactionProductVector * products) 2865 { 2782 { 2866 // return product when nucleus is destroy << 2783 // return product when nucleus is destroyed, i.e. charge=0 2867 G4double Esecondaries(0.); 2784 G4double Esecondaries(0.); 2868 G4LorentzVector psecondaries; << 2785 std::vector<G4KineticTrack *>::iterator iter; 2869 std::vector<G4KineticTrack *>::const_iter << 2870 std::vector<G4ReactionProduct *>::const_i << 2871 decayKTV.Decay(&theFinalState); 2786 decayKTV.Decay(&theFinalState); 2872 << 2787 //PrintKTVector(&theFinalState, " theFinalState @ void Nucl"); 2873 for(iter = theFinalState.cbegin(); iter ! << 2788 for(iter = theFinalState.begin(); iter != theFinalState.end(); ++iter) 2874 { 2789 { 2875 G4ReactionProduct * aNew = new G4Reac 2790 G4ReactionProduct * aNew = new G4ReactionProduct((*iter)->GetDefinition()); 2876 aNew->SetMomentum((*iter)->Get4Moment 2791 aNew->SetMomentum((*iter)->Get4Momentum().vect()); 2877 aNew->SetTotalEnergy((*iter)->Get4Mom 2792 aNew->SetTotalEnergy((*iter)->Get4Momentum().e()); 2878 aNew->SetCreatorModelID(theBIC_ID); << 2879 aNew->SetParentResonanceDef((*iter)->GetPar << 2880 aNew->SetParentResonanceID((*iter)->GetPare << 2881 Esecondaries +=(*iter)->Get4Momentum( 2793 Esecondaries +=(*iter)->Get4Momentum().e(); 2882 psecondaries +=(*iter)->Get4Momentum( << 2883 aNew->SetNewlyAdded(true); 2794 aNew->SetNewlyAdded(true); 2884 //G4cout << " Particle Ekin " << aNew 2795 //G4cout << " Particle Ekin " << aNew->GetKineticEnergy() << G4endl; 2885 products->push_back(aNew); 2796 products->push_back(aNew); 2886 } 2797 } 2887 2798 >> 2799 //PrintKTVector(&theCapturedList, " theCapturedList @ void Nucl"); >> 2800 for(iter = theCapturedList.begin(); iter != theCapturedList.end(); ++iter) >> 2801 { >> 2802 G4ReactionProduct * aNew = new G4ReactionProduct((*iter)->GetDefinition()); >> 2803 aNew->SetMomentum((*iter)->Get4Momentum().vect()); >> 2804 aNew->SetTotalEnergy((*iter)->Get4Momentum().e()); >> 2805 Esecondaries +=(*iter)->Get4Momentum().e(); >> 2806 aNew->SetNewlyAdded(true); >> 2807 //G4cout << " Particle Ekin " << aNew->GetKineticEnergy() << G4endl; >> 2808 products->push_back(aNew); >> 2809 } 2888 // pull out late particles from collision 2810 // pull out late particles from collisions 2889 //theCollisionMgr->Print(); 2811 //theCollisionMgr->Print(); 2890 while(theCollisionMgr->Entries() > 0) << 2812 while(theCollisionMgr->Entries() > 0) 2891 { 2813 { 2892 G4CollisionInitialState * 2814 G4CollisionInitialState * 2893 collision = theCollisionMgr->GetNextC 2815 collision = theCollisionMgr->GetNextCollision(); 2894 2816 2895 if ( ! collision->GetTargetCollection 2817 if ( ! collision->GetTargetCollection().size() ){ 2896 G4KineticTrackVector * lates = co 2818 G4KineticTrackVector * lates = collision->GetFinalState(); 2897 if ( lates->size() == 1 ) { 2819 if ( lates->size() == 1 ) { 2898 G4KineticTrack * atrack=*(lat 2820 G4KineticTrack * atrack=*(lates->begin()); 2899 //PrintKTVector(atrack, " lat << 2900 << 2901 G4ReactionProduct * aNew = ne 2821 G4ReactionProduct * aNew = new G4ReactionProduct(atrack->GetDefinition()); 2902 aNew->SetMomentum(atrack->Get 2822 aNew->SetMomentum(atrack->Get4Momentum().vect()); 2903 aNew->SetTotalEnergy(atrack-> 2823 aNew->SetTotalEnergy(atrack->Get4Momentum().e()); 2904 aNew->SetCreatorModelID(atrac << 2905 aNew->SetParentResonanceDef(atrack->GetPa << 2906 aNew->SetParentResonanceID(atrack->GetPar << 2907 Esecondaries +=atrack->Get4Mo 2824 Esecondaries +=atrack->Get4Momentum().e(); 2908 psecondaries +=atrack->Get4Mo << 2909 aNew->SetNewlyAdded(true); 2825 aNew->SetNewlyAdded(true); >> 2826 //G4cout << " Particle Ekin " << aNew->GetKineticEnergy() << G4endl; 2910 products->push_back(aNew); 2827 products->push_back(aNew); 2911 2828 2912 } 2829 } 2913 } 2830 } 2914 theCollisionMgr->RemoveCollision(coll 2831 theCollisionMgr->RemoveCollision(collision); 2915 2832 2916 } 2833 } 2917 2834 2918 // decay must be after loop on Collisions 2835 // decay must be after loop on Collisions, and Decay() will delete entries in theSecondaryList, refered 2919 // to by Collisions. 2836 // to by Collisions. 2920 decayKTV.Decay(&theSecondaryList); 2837 decayKTV.Decay(&theSecondaryList); 2921 << 2838 //PrintKTVector(&theSecondaryList, " secondaires @ void Nucl"); 2922 // Correct for momentum transfered to Nuc << 2839 for(iter = theSecondaryList.begin(); iter != theSecondaryList.end(); ++iter) 2923 G4ThreeVector transferCorrection(0); << 2924 if ( (theSecondaryList.size() + theCaptur << 2925 { << 2926 transferCorrection= theMomentumTransfer << 2927 } << 2928 << 2929 for(iter = theSecondaryList.cbegin(); ite << 2930 { 2840 { 2931 G4ReactionProduct * aNew = new G4Reac 2841 G4ReactionProduct * aNew = new G4ReactionProduct((*iter)->GetDefinition()); 2932 (*iter)->Update4Momentum((*iter)->Get << 2933 aNew->SetMomentum((*iter)->Get4Moment 2842 aNew->SetMomentum((*iter)->Get4Momentum().vect()); 2934 aNew->SetTotalEnergy((*iter)->Get4Mom 2843 aNew->SetTotalEnergy((*iter)->Get4Momentum().e()); 2935 aNew->SetCreatorModelID(theBIC_ID); << 2936 aNew->SetParentResonanceDef((*iter)->GetPar << 2937 aNew->SetParentResonanceID((*iter)->GetPare << 2938 Esecondaries +=(*iter)->Get4Momentum( 2844 Esecondaries +=(*iter)->Get4Momentum().e(); 2939 psecondaries +=(*iter)->Get4Momentum( << 2940 if ( (*iter)->IsParticipant() ) aNew- << 2941 products->push_back(aNew); << 2942 } << 2943 << 2944 for(iter = theCapturedList.cbegin(); iter << 2945 { << 2946 G4ReactionProduct * aNew = new G4Reac << 2947 (*iter)->Update4Momentum((*iter)->Get << 2948 aNew->SetMomentum((*iter)->Get4Moment << 2949 aNew->SetTotalEnergy((*iter)->Get4Mom << 2950 aNew->SetCreatorModelID(theBIC_ID); << 2951 aNew->SetParentResonanceDef((*iter)->GetPar << 2952 aNew->SetParentResonanceID((*iter)->GetPare << 2953 Esecondaries +=(*iter)->Get4Momentum( << 2954 psecondaries +=(*iter)->Get4Momentum( << 2955 aNew->SetNewlyAdded(true); 2845 aNew->SetNewlyAdded(true); >> 2846 //G4cout << " Particle Ekin " << aNew->GetKineticEnergy() << G4endl; 2956 products->push_back(aNew); 2847 products->push_back(aNew); 2957 } 2848 } 2958 2849 2959 G4double SumMassNucleons(0.); 2850 G4double SumMassNucleons(0.); 2960 G4LorentzVector pNucleons(0.); << 2851 for(iter = theTargetList.begin(); iter != theTargetList.end(); ++iter) 2961 for(iter = theTargetList.cbegin(); iter ! << 2962 { 2852 { 2963 SumMassNucleons += (*iter)->GetDefini 2853 SumMassNucleons += (*iter)->GetDefinition()->GetPDGMass(); 2964 pNucleons += (*iter)->Get4Momentum(); << 2965 } 2854 } 2966 2855 2967 G4double Ekinetic=theProjectile4Momentum. 2856 G4double Ekinetic=theProjectile4Momentum.e() + initial_nuclear_mass - Esecondaries - SumMassNucleons; 2968 #ifdef debug_BIC_FillVoidnucleus << 2857 //G4cout << "Void nucleus nucleons : "<<theTargetList.size() << " , T: " << Ekinetic << G4endl; 2969 G4LorentzVector deltaP=theProjectile4 << 2858 if (Ekinetic > 0.){ 2970 psecondaries - pNucle << 2859 if (theTargetList.size() ) Ekinetic /= theTargetList.size(); 2971 //G4cout << "BIC::FillVoidNucleus() n << 2972 // ", deltaP " << deltaP << " de << 2973 #endif << 2974 if (Ekinetic > 0. && theTargetList.size() << 2975 Ekinetic /= theTargetList.size(); << 2976 } else { 2860 } else { 2977 G4double Ekineticrdm(0); << 2861 //G4cout << " zero or negative kinetic energy, use random: " << Ekinetic<< G4endl; 2978 if (theTargetList.size()) Ekineticrdm << 2862 Ekinetic= ( 0.1 + G4UniformRand()*5.) * MeV; // violate Energy conservation by small amount here 2979 G4double TotalEkin(Ekineticrdm); << 2863 } 2980 for (rpiter=products->cbegin(); rpite << 2864 //G4cout << " using T per nucleon: " << Ekinetic << G4endl; 2981 TotalEkin+=(*rpiter)->GetKineticEne << 2982 } << 2983 G4double correction(1.); << 2984 if ( std::abs(Ekinetic) < 20*perCent << 2985 correction=1. + (Ekinetic-Ekineticr << 2986 } << 2987 #ifdef debug_G4BinaryCascade << 2988 else { << 2989 G4cout << "BLIC::FillVoidNucleus() fa << 2990 } << 2991 #endif << 2992 << 2993 for (rpiter=products->cbegin(); rpite << 2994 (*rpiter)->SetKineticEnergy((*rpiter) << 2995 (*rpiter)->SetMomentum((*rpiter)->G << 2996 << 2997 } << 2998 << 2999 Ekinetic=Ekineticrdm*correction; << 3000 if (theTargetList.size())Ekinetic /= << 3001 << 3002 } << 3003 2865 3004 for(iter = theTargetList.cbegin(); iter ! << 2866 for(iter = theTargetList.begin(); iter != theTargetList.end(); ++iter) 3005 // set Nucleon it to be hit - as it is << 2867 { 3006 (*iter)->Hit(); << 3007 G4ReactionProduct * aNew = new G4Reac 2868 G4ReactionProduct * aNew = new G4ReactionProduct((*iter)->GetDefinition()); 3008 aNew->SetKineticEnergy(Ekinetic); << 2869 G4double mass=(*iter)->GetDefinition()->GetPDGMass(); 3009 aNew->SetMomentum(aNew->GetTotalMomen << 2870 G4double p=std::sqrt(sqr(Ekinetic) + 2.*Ekinetic*mass); >> 2871 aNew->SetMomentum(p*(*iter)->Get4Momentum().vect().unit()); >> 2872 aNew->SetTotalEnergy(mass+Ekinetic); 3010 aNew->SetNewlyAdded(true); 2873 aNew->SetNewlyAdded(true); 3011 aNew->SetCreatorModelID(theBIC_ID); << 2874 //G4cout << " Particle Ekin " << aNew->GetKineticEnergy() << G4endl; 3012 aNew->SetParentResonanceDef((*iter)->GetPar << 3013 aNew->SetParentResonanceID((*iter)->GetPare << 3014 products->push_back(aNew); 2875 products->push_back(aNew); 3015 Esecondaries += aNew->GetTotalEnergy( << 3016 psecondaries += G4LorentzVector(aNew- << 3017 } << 3018 psecondaries=G4LorentzVector(0); << 3019 for (rpiter=products->cbegin(); rpiter!=p << 3020 psecondaries += G4LorentzVector((*rpite << 3021 } 2876 } 3022 2877 3023 G4LorentzVector initial4Mom=theProjectile << 3024 << 3025 //G4cout << "::FillVoidNucleus()final e/p << 3026 // << " final " << psecondaries << " del << 3027 << 3028 G4ThreeVector SumMom=psecondaries.vect(); << 3029 << 3030 SumMom=initial4Mom.vect()-SumMom; << 3031 G4int loopcount(0); << 3032 << 3033 // reverse_iterator reverse - start to co << 3034 while ( SumMom.mag() > 0.1*MeV && loopcou << 3035 { << 3036 G4int index=(G4int)products->size(); << 3037 for (auto reverse=products->crbegin(); << 3038 SumMom=initial4Mom.vect(); << 3039 for (rpiter=products->cbegin(); rpite << 3040 SumMom-=(*rpiter)->GetMomentum(); << 3041 } << 3042 G4double p=((*reverse)->GetMomentum() << 3043 (*reverse)->SetMomentum( p*(((*rever << 3044 } << 3045 } << 3046 return products; 2878 return products; 3047 } 2879 } 3048 << 3049 G4ReactionProductVector * G4BinaryCascade::Hi 2880 G4ReactionProductVector * G4BinaryCascade::HighEnergyModelFSProducts(G4ReactionProductVector * products, 3050 G4KineticTrackVector * secondaries) 2881 G4KineticTrackVector * secondaries) 3051 { 2882 { 3052 for(auto iter = secondaries->cbegin(); it << 2883 std::vector<G4KineticTrack *>::iterator iter; >> 2884 for(iter = secondaries->begin(); iter != secondaries->end(); ++iter) 3053 { 2885 { 3054 G4ReactionProduct * aNew = new G4Reac 2886 G4ReactionProduct * aNew = new G4ReactionProduct((*iter)->GetDefinition()); 3055 aNew->SetMomentum((*iter)->Get4Moment 2887 aNew->SetMomentum((*iter)->Get4Momentum().vect()); 3056 aNew->SetTotalEnergy((*iter)->Get4Mom 2888 aNew->SetTotalEnergy((*iter)->Get4Momentum().e()); 3057 aNew->SetNewlyAdded(true); 2889 aNew->SetNewlyAdded(true); 3058 aNew->SetCreatorModelID((*iter)->GetC << 3059 aNew->SetParentResonanceDef((*iter)-> << 3060 aNew->SetParentResonanceID((*iter)->G << 3061 //G4cout << " Particle Ekin " << aNew 2890 //G4cout << " Particle Ekin " << aNew->GetKineticEnergy() << G4endl; 3062 products->push_back(aNew); 2891 products->push_back(aNew); 3063 } 2892 } 3064 const G4ParticleDefinition* fragment = 0; << 2893 G4ParticleDefinition* fragment = 0; 3065 if (currentA == 1 && currentZ == 0) { 2894 if (currentA == 1 && currentZ == 0) { 3066 fragment = G4Neutron::NeutronDefiniti 2895 fragment = G4Neutron::NeutronDefinition(); 3067 } else if (currentA == 1 && currentZ == 1 2896 } else if (currentA == 1 && currentZ == 1) { 3068 fragment = G4Proton::ProtonDefinition 2897 fragment = G4Proton::ProtonDefinition(); 3069 } else if (currentA == 2 && currentZ == 1 2898 } else if (currentA == 2 && currentZ == 1) { 3070 fragment = G4Deuteron::DeuteronDefini 2899 fragment = G4Deuteron::DeuteronDefinition(); 3071 } else if (currentA == 3 && currentZ == 1 2900 } else if (currentA == 3 && currentZ == 1) { 3072 fragment = G4Triton::TritonDefinition 2901 fragment = G4Triton::TritonDefinition(); 3073 } else if (currentA == 3 && currentZ == 2 2902 } else if (currentA == 3 && currentZ == 2) { 3074 fragment = G4He3::He3Definition(); 2903 fragment = G4He3::He3Definition(); 3075 } else if (currentA == 4 && currentZ == 2 2904 } else if (currentA == 4 && currentZ == 2) { 3076 fragment = G4Alpha::AlphaDefinition() 2905 fragment = G4Alpha::AlphaDefinition();; 3077 } else { 2906 } else { 3078 fragment = 2907 fragment = 3079 G4ParticleTable::GetParticleT 2908 G4ParticleTable::GetParticleTable()->GetIonTable()->GetIon(currentZ,currentA,0.0); 3080 } 2909 } 3081 if (fragment != 0) { 2910 if (fragment != 0) { 3082 G4ReactionProduct * theNew = new G4Re 2911 G4ReactionProduct * theNew = new G4ReactionProduct(fragment); 3083 theNew->SetMomentum(G4ThreeVector(0,0 2912 theNew->SetMomentum(G4ThreeVector(0,0,0)); 3084 theNew->SetTotalEnergy(massInNucleus) 2913 theNew->SetTotalEnergy(massInNucleus); 3085 theNew->SetCreatorModelID(theBIC_ID); << 3086 //theNew->SetFormationTime(??0.??); 2914 //theNew->SetFormationTime(??0.??); 3087 //G4cout << " Nucleus (" << currentZ 2915 //G4cout << " Nucleus (" << currentZ << ","<< currentA << "), mass "<< massInNucleus << G4endl; 3088 products->push_back(theNew); 2916 products->push_back(theNew); 3089 } 2917 } 3090 return products; 2918 return products; 3091 } 2919 } 3092 2920 3093 void G4BinaryCascade::PrintWelcomeMessage() 2921 void G4BinaryCascade::PrintWelcomeMessage() 3094 { 2922 { 3095 G4cout <<"Thank you for using G4BinaryCas 2923 G4cout <<"Thank you for using G4BinaryCascade. "<<G4endl; 3096 } 2924 } 3097 2925 3098 //------------------------------------------- 2926 //---------------------------------------------------------------------------- 3099 void G4BinaryCascade::DebugApplyCollisionFail 2927 void G4BinaryCascade::DebugApplyCollisionFail(G4CollisionInitialState * collision, 3100 G4KineticTrackVector * products) 2928 G4KineticTrackVector * products) 3101 { 2929 { 3102 G4bool havePion=false; 2930 G4bool havePion=false; 3103 if (products) 2931 if (products) 3104 { 2932 { 3105 for ( auto i =products->cbegin(); i != << 2933 for ( std::vector<G4KineticTrack *>::iterator i =products->begin(); i != products->end(); i++) 3106 { 2934 { 3107 G4int PDGcode=std::abs((*i)->GetDefi 2935 G4int PDGcode=std::abs((*i)->GetDefinition()->GetPDGEncoding()); 3108 if (std::abs(PDGcode)==211 || PDGcod 2936 if (std::abs(PDGcode)==211 || PDGcode==111 ) havePion=true; 3109 } 2937 } 3110 } 2938 } 3111 if ( !products || havePion) 2939 if ( !products || havePion) 3112 { 2940 { 3113 const G4BCAction &action= *collision->G << 2941 G4cout << " Collision " << collision << ", type: "<< typeid(*collision->GetGenerator()).name() 3114 G4cout << " Collision " << collision << << 3115 << ", with 2942 << ", with NO products! " <<G4endl; 3116 G4cout << G4endl<<"Initial condition ar 2943 G4cout << G4endl<<"Initial condition are these:"<<G4endl; 3117 G4cout << "proj: "<<collision->GetPrima 2944 G4cout << "proj: "<<collision->GetPrimary()->GetDefinition()->GetParticleName()<<G4endl; 3118 PrintKTVector(collision->GetPrimary()); 2945 PrintKTVector(collision->GetPrimary()); 3119 for(std::size_t it=0; it<collision->Get << 2946 for(size_t it=0; it<collision->GetTargetCollection().size(); it++) 3120 { 2947 { 3121 G4cout << "targ: " 2948 G4cout << "targ: " 3122 <<collision->GetTargetCollecti 2949 <<collision->GetTargetCollection()[it]->GetDefinition()->GetParticleName()<<G4endl; 3123 } 2950 } 3124 PrintKTVector(&collision->GetTargetColl 2951 PrintKTVector(&collision->GetTargetCollection(),std::string(" Target particles")); 3125 } 2952 } 3126 // if ( lateParticleCollision ) G4cout << 2953 // if ( lateParticleCollision ) G4cout << " Added late particle--------------------------"<<G4endl; 3127 // if ( lateParticleCollision && products 2954 // if ( lateParticleCollision && products ) PrintKTVector(products, " reaction products"); 3128 } 2955 } 3129 2956 3130 //------------------------------------------- 2957 //---------------------------------------------------------------------------- 3131 2958 3132 G4bool G4BinaryCascade::CheckChargeAndBaryonN << 2959 G4bool G4BinaryCascade::CheckChargeAndBaryonNumber(G4String where) 3133 { 2960 { 3134 static G4int lastdA(0), lastdZ(0); 2961 static G4int lastdA(0), lastdZ(0); 3135 G4int iStateA = the3DNucleus->GetMassNumbe 2962 G4int iStateA = the3DNucleus->GetMassNumber() + projectileA; 3136 G4int iStateZ = the3DNucleus->GetCharge() 2963 G4int iStateZ = the3DNucleus->GetCharge() + projectileZ; 3137 2964 3138 G4int fStateA(0); 2965 G4int fStateA(0); 3139 G4int fStateZ(0); 2966 G4int fStateZ(0); 3140 2967 >> 2968 std::vector<G4KineticTrack *>::iterator i; 3141 G4int CapturedA(0), CapturedZ(0); 2969 G4int CapturedA(0), CapturedZ(0); 3142 G4int secsA(0), secsZ(0); 2970 G4int secsA(0), secsZ(0); 3143 for (auto i=theCapturedList.cbegin(); i!=t << 2971 for ( i=theCapturedList.begin(); i!=theCapturedList.end(); ++i) { 3144 CapturedA += (*i)->GetDefinition()->Get 2972 CapturedA += (*i)->GetDefinition()->GetBaryonNumber(); 3145 CapturedZ += G4lrint((*i)->GetDefinitio 2973 CapturedZ += G4lrint((*i)->GetDefinition()->GetPDGCharge()/eplus); 3146 } 2974 } 3147 2975 3148 for (auto i=theSecondaryList.cbegin(); i!= << 2976 for ( i=theSecondaryList.begin(); i!=theSecondaryList.end(); ++i) { 3149 if ( (*i)->GetState() != G4KineticTrack 2977 if ( (*i)->GetState() != G4KineticTrack::inside ) { 3150 secsA += (*i)->GetDefinition()->GetB 2978 secsA += (*i)->GetDefinition()->GetBaryonNumber(); 3151 secsZ += G4lrint((*i)->GetDefinition 2979 secsZ += G4lrint((*i)->GetDefinition()->GetPDGCharge()/eplus); 3152 } 2980 } 3153 } 2981 } 3154 2982 3155 for (auto i=theFinalState.cbegin(); i!=the << 2983 for ( i=theFinalState.begin(); i!=theFinalState.end(); ++i) { 3156 fStateA += (*i)->GetDefinition()->GetBa 2984 fStateA += (*i)->GetDefinition()->GetBaryonNumber(); 3157 fStateZ += G4lrint((*i)->GetDefinition( 2985 fStateZ += G4lrint((*i)->GetDefinition()->GetPDGCharge()/eplus); 3158 } 2986 } 3159 2987 3160 G4int deltaA= iStateA - secsA - fStateA - 2988 G4int deltaA= iStateA - secsA - fStateA -currentA - lateA; 3161 G4int deltaZ= iStateZ - secsZ - fStateZ - 2989 G4int deltaZ= iStateZ - secsZ - fStateZ -currentZ - lateZ; 3162 2990 3163 #ifdef debugCheckChargeAndBaryonNumberverbose << 3164 G4cout << where <<" A: iState= "<< iStateA< << 3165 G4cout << where <<" Z: iState= "<< iStateZ< << 3166 #endif << 3167 << 3168 if (deltaA != 0 || deltaZ!=0 ) { 2991 if (deltaA != 0 || deltaZ!=0 ) { 3169 if (deltaA != lastdA || deltaZ != lastd 2992 if (deltaA != lastdA || deltaZ != lastdZ ) { 3170 G4cout << "baryon/charge imbalance - 2993 G4cout << "baryon/charge imbalance - " << where << G4endl 3171 << "deltaA " <<deltaA<<", iSta 2994 << "deltaA " <<deltaA<<", iStateA "<<iStateA<< ", CapturedA "<<CapturedA <<", secsA "<<secsA 3172 << ", fStateA "<<fStateA << ", 2995 << ", fStateA "<<fStateA << ", currentA "<<currentA << ", lateA "<<lateA << G4endl 3173 << "deltaZ "<<deltaZ<<", iStat 2996 << "deltaZ "<<deltaZ<<", iStateZ "<<iStateZ<< ", CapturedZ "<<CapturedZ <<", secsZ "<<secsZ 3174 << ", fStateZ "<<fStateZ << ", 2997 << ", fStateZ "<<fStateZ << ", currentZ "<<currentZ << ", lateZ "<<lateZ << G4endl<< G4endl; 3175 lastdA=deltaA; 2998 lastdA=deltaA; 3176 lastdZ=deltaZ; 2999 lastdZ=deltaZ; 3177 } 3000 } 3178 } else { lastdA=lastdZ=0;} 3001 } else { lastdA=lastdZ=0;} 3179 3002 3180 return true; 3003 return true; 3181 } 3004 } 3182 //------------------------------------------- 3005 //---------------------------------------------------------------------------- 3183 void G4BinaryCascade::DebugApplyCollision(G4C 3006 void G4BinaryCascade::DebugApplyCollision(G4CollisionInitialState * collision, 3184 G4KineticTrackVector * products) 3007 G4KineticTrackVector * products) 3185 { 3008 { 3186 3009 3187 PrintKTVector(collision->GetPrimary(),std 3010 PrintKTVector(collision->GetPrimary(),std::string(" Primary particle")); 3188 PrintKTVector(&collision->GetTargetCollec 3011 PrintKTVector(&collision->GetTargetCollection(),std::string(" Target particles")); 3189 PrintKTVector(products,std::string(" Scat 3012 PrintKTVector(products,std::string(" Scatterer products")); 3190 3013 3191 #ifdef dontUse 3014 #ifdef dontUse 3192 G4double thisExcitation(0); 3015 G4double thisExcitation(0); 3193 // excitation energy from this collision 3016 // excitation energy from this collision 3194 // initial state: 3017 // initial state: 3195 G4double initial(0); 3018 G4double initial(0); 3196 G4KineticTrack * kt=collision->GetPrimary 3019 G4KineticTrack * kt=collision->GetPrimary(); 3197 initial += kt->Get4Momentum().e(); 3020 initial += kt->Get4Momentum().e(); 3198 3021 3199 G4RKPropagation * RKprop=(G4RKPropagation 3022 G4RKPropagation * RKprop=(G4RKPropagation *)thePropagator; 3200 3023 3201 initial += RKprop->GetField(kt->GetDefin 3024 initial += RKprop->GetField(kt->GetDefinition()->GetPDGEncoding(),kt->GetPosition()); 3202 initial -= RKprop->GetBarrier(kt->GetDef 3025 initial -= RKprop->GetBarrier(kt->GetDefinition()->GetPDGEncoding()); 3203 G4cout << "prim. E/field/Barr/Sum " << kt 3026 G4cout << "prim. E/field/Barr/Sum " << kt->Get4Momentum().e() 3204 << " " << RKprop->GetFi 3027 << " " << RKprop->GetField(kt->GetDefinition()->GetPDGEncoding(),kt->GetPosition()) 3205 << " " << RKprop->GetBa 3028 << " " << RKprop->GetBarrier(kt->GetDefinition()->GetPDGEncoding()) 3206 << " " << initial << G4 3029 << " " << initial << G4endl;; 3207 3030 3208 G4KineticTrackVector ktv=collision->GetTa 3031 G4KineticTrackVector ktv=collision->GetTargetCollection(); 3209 for ( unsigned int it=0; it < ktv.size(); << 3032 for ( unsigned int it=0; it < ktv.size(); it++) 3210 { 3033 { 3211 kt=ktv[it]; 3034 kt=ktv[it]; 3212 initial += kt->Get4Momentum().e(); 3035 initial += kt->Get4Momentum().e(); 3213 thisExcitation += kt->GetDefinition() 3036 thisExcitation += kt->GetDefinition()->GetPDGMass() 3214 - kt->Get4Momentum() 3037 - kt->Get4Momentum().e() 3215 - RKprop->GetField(k 3038 - RKprop->GetField(kt->GetDefinition()->GetPDGEncoding(),kt->GetPosition()); 3216 // initial += RKprop->GetField(k 3039 // initial += RKprop->GetField(kt->GetDefinition()->GetPDGEncoding(),kt->GetPosition()); 3217 // initial -= RKprop->GetBarrier 3040 // initial -= RKprop->GetBarrier(kt->GetDefinition()->GetPDGEncoding()); 3218 G4cout << "Targ. def/E/field/Barr/Sum 3041 G4cout << "Targ. def/E/field/Barr/Sum " << kt->GetDefinition()->GetPDGEncoding() 3219 << " " << kt->Get4Momentum( 3042 << " " << kt->Get4Momentum().e() 3220 << " " << RKprop->GetField( 3043 << " " << RKprop->GetField(kt->GetDefinition()->GetPDGEncoding(),kt->GetPosition()) 3221 << " " << RKprop->GetBarrie 3044 << " " << RKprop->GetBarrier(kt->GetDefinition()->GetPDGEncoding()) 3222 << " " << initial <<" Excit 3045 << " " << initial <<" Excit " << thisExcitation << G4endl;; 3223 } 3046 } 3224 3047 3225 G4double fin(0); << 3048 G4double final(0); 3226 G4double mass_out(0); 3049 G4double mass_out(0); 3227 G4int product_barions(0); 3050 G4int product_barions(0); 3228 if ( products ) 3051 if ( products ) 3229 { 3052 { 3230 for ( unsigned int it=0; it < product << 3053 for ( unsigned int it=0; it < products->size(); it++) 3231 { 3054 { 3232 kt=(*products)[it]; 3055 kt=(*products)[it]; 3233 fin += kt->Get4Momentum().e(); << 3056 final += kt->Get4Momentum().e(); 3234 fin += RKprop->GetField(kt->GetD << 3057 final += RKprop->GetField(kt->GetDefinition()->GetPDGEncoding(),kt->GetPosition()); 3235 fin += RKprop->GetBarrier(kt->Ge << 3058 final += RKprop->GetBarrier(kt->GetDefinition()->GetPDGEncoding()); 3236 if ( kt->GetDefinition()->GetBary 3059 if ( kt->GetDefinition()->GetBaryonNumber()==1 ) product_barions++; 3237 mass_out += kt->GetDefinition()-> 3060 mass_out += kt->GetDefinition()->GetPDGMass(); 3238 G4cout << "sec. def/E/field/Barr/ 3061 G4cout << "sec. def/E/field/Barr/Sum " << kt->GetDefinition()->GetPDGEncoding() 3239 << " " << kt->Get4Moment 3062 << " " << kt->Get4Momentum().e() 3240 << " " << RKprop->GetFie 3063 << " " << RKprop->GetField(kt->GetDefinition()->GetPDGEncoding(),kt->GetPosition()) 3241 << " " << RKprop->GetBar 3064 << " " << RKprop->GetBarrier(kt->GetDefinition()->GetPDGEncoding()) 3242 << " " << fin << G4endl; << 3065 << " " << final << G4endl;; 3243 } 3066 } 3244 } 3067 } 3245 3068 3246 3069 3247 G4int finalA = currentA; 3070 G4int finalA = currentA; 3248 G4int finalZ = currentZ; 3071 G4int finalZ = currentZ; 3249 if ( products ) 3072 if ( products ) 3250 { 3073 { 3251 finalA -= product_barions; 3074 finalA -= product_barions; 3252 finalZ -= GetTotalCharge(*products); 3075 finalZ -= GetTotalCharge(*products); 3253 } 3076 } 3254 G4double delta = GetIonMass(currentZ,curr << 3077 G4double delta = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(currentZ,currentA) >> 3078 - (G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(finalZ,finalA) >> 3079 + mass_out); 3255 G4cout << " current/final a,z " << curren 3080 G4cout << " current/final a,z " << currentA << " " << currentZ << " "<< finalA<< " "<< finalZ 3256 << " delta-mass " << delta<<G4en 3081 << " delta-mass " << delta<<G4endl; 3257 fin+=delta; << 3082 final+=delta; 3258 mass_out = GetIonMass(finalZ,finalA); << 3083 mass_out = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(finalZ,finalA); 3259 G4cout << " initE/ E_out/ Mfinal/ Excit " 3084 G4cout << " initE/ E_out/ Mfinal/ Excit " << currentInitialEnergy 3260 << " " << fin << " " << 3085 << " " << final << " " 3261 << mass_out<<" " 3086 << mass_out<<" " 3262 << currentInitialEnergy - fin - << 3087 << currentInitialEnergy - final - mass_out 3263 << G4endl; 3088 << G4endl; 3264 currentInitialEnergy-=fin; << 3089 currentInitialEnergy-=final; 3265 #endif 3090 #endif 3266 } 3091 } 3267 3092 3268 //------------------------------------------- 3093 //---------------------------------------------------------------------------- 3269 G4bool G4BinaryCascade::DebugFinalEpConservat << 3094 G4bool G4BinaryCascade::DebugEpConservation(const G4HadProjectile & aTrack, 3270 G4ReactionProductVector* products) 3095 G4ReactionProductVector* products) 3271 //------------------------------------------- 3096 //---------------------------------------------------------------------------- 3272 { 3097 { >> 3098 G4ReactionProductVector::iterator iter; 3273 G4double Efinal(0); 3099 G4double Efinal(0); 3274 G4ThreeVector pFinal(0); 3100 G4ThreeVector pFinal(0); 3275 if (std::abs(theParticleChange.GetWeightC 3101 if (std::abs(theParticleChange.GetWeightChange() -1 ) > 1e-5 ) 3276 { 3102 { 3277 G4cout <<" BIC-weight change " << the 3103 G4cout <<" BIC-weight change " << theParticleChange.GetWeightChange()<< G4endl; 3278 } 3104 } 3279 3105 3280 for(auto iter = products->cbegin(); iter << 3106 for(iter = products->begin(); iter != products->end(); ++iter) 3281 { 3107 { 3282 3108 3283 G4cout << " Secondary E - Ekin / p " << 3109 // G4cout << " Secondary E - Ekin / p " << 3284 (*iter)->GetDefinition()->GetPa << 3110 // (*iter)->GetDefinition()->GetParticleName() << " " << 3285 (*iter)->GetTotalEnergy() << " << 3111 // (*iter)->GetTotalEnergy() << " - " << 3286 (*iter)->GetKineticEnergy()<< " << 3112 // (*iter)->GetKineticEnergy()<< " / " << 3287 (*iter)->GetMomentum().x() << " << 3113 // (*iter)->GetMomentum().x() << " " << 3288 (*iter)->GetMomentum().y() << " << 3114 // (*iter)->GetMomentum().y() << " " << 3289 (*iter)->GetMomentum().z() << G << 3115 // (*iter)->GetMomentum().z() << G4endl; 3290 Efinal += (*iter)->GetTotalEnergy(); 3116 Efinal += (*iter)->GetTotalEnergy(); 3291 pFinal += (*iter)->GetMomentum(); 3117 pFinal += (*iter)->GetMomentum(); 3292 } 3118 } 3293 3119 3294 G4cout << "e outgoing/ total : " << Efi << 3120 // G4cout << "e outgoing/ total : " << Efinal << " " << Efinal+GetFinal4Momentum().e()<< G4endl; 3295 G4cout << "BIC E/p delta " << << 3121 G4cout << "BIC E/p delta " << 3296 (aTrack.Get4Momentum().e()+theIni << 3122 (aTrack.Get4Momentum().e()+ the3DNucleus->GetMass() - Efinal)/MeV << 3297 " MeV / mom " << (aTrack.Get4Mome 3123 " MeV / mom " << (aTrack.Get4Momentum() - pFinal ) /MeV << G4endl; 3298 3124 3299 return (aTrack.Get4Momentum().e() + theIn << 3125 return (aTrack.Get4Momentum().e() + the3DNucleus->GetMass() - Efinal)/aTrack.Get4Momentum().e() < perCent; 3300 << 3301 } << 3302 //------------------------------------------- << 3303 G4bool G4BinaryCascade::DebugEpConservation(c << 3304 //------------------------------------------- << 3305 { << 3306 G4cout << where << G4endl; << 3307 G4LorentzVector psecs, ptgts, pcpts << 3308 if (std::abs(theParticleChange.GetWeightC << 3309 { << 3310 G4cout <<" BIC-weight change " << the << 3311 } << 3312 << 3313 std::vector<G4KineticTrack *>::const_iter << 3314 for(ktiter = theSecondaryList.cbegin(); k << 3315 { << 3316 << 3317 G4cout << " Secondary E - Ekin << 3318 (*ktiter)->GetDefinition()-> << 3319 (*ktiter)->Get4Momentum().e( << 3320 (*ktiter)->Get4Momentum().e( << 3321 (*ktiter)->Get4Momentum().ve << 3322 psecs += (*ktiter)->Get4Momentum() << 3323 } << 3324 << 3325 for(ktiter = theTargetList.cbegin(); ktit << 3326 { << 3327 << 3328 G4cout << " Target E - Ekin / p << 3329 (*ktiter)->GetDefinition()-> << 3330 (*ktiter)->Get4Momentum().e( << 3331 (*ktiter)->Get4Momentum().e( << 3332 (*ktiter)->Get4Momentum().ve << 3333 ptgts += (*ktiter)->Get4Momentum() << 3334 } << 3335 << 3336 for(ktiter = theCapturedList.cbegin(); kt << 3337 { << 3338 << 3339 G4cout << " Captured E - Ekin << 3340 (*ktiter)->GetDefinition()- << 3341 (*ktiter)->Get4Momentum().e << 3342 (*ktiter)->Get4Momentum().e << 3343 (*ktiter)->Get4Momentum().v << 3344 pcpts += (*ktiter)->Get4Momentum( << 3345 } << 3346 << 3347 for(ktiter = theFinalState.cbegin(); ktit << 3348 { << 3349 << 3350 G4cout << " Finals E - Ekin / << 3351 (*ktiter)->GetDefinition()- << 3352 (*ktiter)->Get4Momentum().e << 3353 (*ktiter)->Get4Momentum().e << 3354 (*ktiter)->Get4Momentum().v << 3355 pfins += (*ktiter)->Get4Momentum( << 3356 } << 3357 << 3358 G4cout << " Secondaries " << psecs << " << 3359 <<" Captured " << pcpts << ", Fi << 3360 <<" Sum " << psecs + ptgts + pcpts << 3361 <<" Sum+PTransfer " << psecs + ptgt << 3362 << G4endl<< G4endl; << 3363 << 3364 << 3365 return true; << 3366 3126 3367 } 3127 } 3368 3128 3369 //------------------------------------------- 3129 //---------------------------------------------------------------------------- 3370 G4ReactionProductVector * G4BinaryCascade::Pr 3130 G4ReactionProductVector * G4BinaryCascade::ProductsAddFakeGamma(G4ReactionProductVector *products ) 3371 //------------------------------------------- 3131 //---------------------------------------------------------------------------- 3372 { 3132 { 3373 // else 3133 // else 3374 // { 3134 // { 3375 // G4ReactionProduct * aNew=0; 3135 // G4ReactionProduct * aNew=0; 3376 // // return nucleus e and p 3136 // // return nucleus e and p 3377 // if (fragment != 0 ) { 3137 // if (fragment != 0 ) { 3378 // aNew = new G4ReactionProduct(G4 3138 // aNew = new G4ReactionProduct(G4Gamma::GammaDefinition()); // we only want to pass e/p 3379 // aNew->SetMomentum(fragment->Get 3139 // aNew->SetMomentum(fragment->GetMomentum().vect()); 3380 // aNew->SetTotalEnergy(fragment-> 3140 // aNew->SetTotalEnergy(fragment->GetMomentum().e()); 3381 // delete fragment; 3141 // delete fragment; 3382 // fragment=0; 3142 // fragment=0; 3383 // } else if (products->size() == 0) { 3143 // } else if (products->size() == 0) { 3384 // // FixMe GF: for testing withou 3144 // // FixMe GF: for testing without precompound, return 1 gamma of 0.01 MeV in +x 3385 //#include "G4Gamma.hh" 3145 //#include "G4Gamma.hh" 3386 // aNew = new G4ReactionProduct(G4 3146 // aNew = new G4ReactionProduct(G4Gamma::GammaDefinition()); 3387 // aNew->SetMomentum(G4ThreeVector 3147 // aNew->SetMomentum(G4ThreeVector(0.01*MeV,0,0)); 3388 // aNew->SetTotalEnergy(0.01*MeV); 3148 // aNew->SetTotalEnergy(0.01*MeV); 3389 // } 3149 // } 3390 // if ( aNew != 0 ) products->push_bac 3150 // if ( aNew != 0 ) products->push_back(aNew); 3391 // } 3151 // } 3392 return products; 3152 return products; 3393 } 3153 } 3394 3154 3395 //------------------------------------------- 3155 //---------------------------------------------------------------------------- 3396 3156