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 #include <algorithm> 26 #include <algorithm> 27 #include <vector> 27 #include <vector> 28 #include <cmath> 28 #include <cmath> 29 #include <numeric> 29 #include <numeric> 30 30 31 #include "G4BinaryLightIonReaction.hh" 31 #include "G4BinaryLightIonReaction.hh" 32 #include "G4PhysicalConstants.hh" 32 #include "G4PhysicalConstants.hh" 33 #include "G4SystemOfUnits.hh" 33 #include "G4SystemOfUnits.hh" 34 #include "G4LorentzVector.hh" 34 #include "G4LorentzVector.hh" 35 #include "G4LorentzRotation.hh" 35 #include "G4LorentzRotation.hh" 36 #include "G4ReactionProductVector.hh" 36 #include "G4ReactionProductVector.hh" 37 #include "G4ping.hh" 37 #include "G4ping.hh" 38 #include "G4Delete.hh" 38 #include "G4Delete.hh" 39 #include "G4Neutron.hh" 39 #include "G4Neutron.hh" 40 #include "G4VNuclearDensity.hh" 40 #include "G4VNuclearDensity.hh" 41 #include "G4FermiMomentum.hh" 41 #include "G4FermiMomentum.hh" 42 #include "G4HadTmpUtil.hh" 42 #include "G4HadTmpUtil.hh" 43 #include "G4PreCompoundModel.hh" 43 #include "G4PreCompoundModel.hh" 44 #include "G4HadronicInteractionRegistry.hh" 44 #include "G4HadronicInteractionRegistry.hh" 45 #include "G4Log.hh" 45 #include "G4Log.hh" 46 #include "G4PhysicsModelCatalog.hh" << 47 #include "G4HadronicParameters.hh" << 48 46 49 G4int G4BinaryLightIonReaction::theBLIR_ID = - << 47 #ifdef G4MULTITHREADED >> 48 G4Mutex G4BinaryLightIonReaction::BLIRMutex = G4MUTEX_INITIALIZER; >> 49 #endif >> 50 G4int G4BinaryLightIonReaction::theBLIR_ID = -1; >> 51 50 52 51 //#define debug_G4BinaryLightIonReaction 53 //#define debug_G4BinaryLightIonReaction 52 //#define debug_BLIR_finalstate 54 //#define debug_BLIR_finalstate 53 //#define debug_BLIR_result 55 //#define debug_BLIR_result 54 56 55 G4BinaryLightIonReaction::G4BinaryLightIonReac 57 G4BinaryLightIonReaction::G4BinaryLightIonReaction(G4VPreCompoundModel* ptr) 56 : G4HadronicInteraction("Binary Light Ion Casc 58 : G4HadronicInteraction("Binary Light Ion Cascade"), 57 theProjectileFragmentation(ptr), 59 theProjectileFragmentation(ptr), 58 pA(0),pZ(0), tA(0),tZ(0),spectatorA(0),spect 60 pA(0),pZ(0), tA(0),tZ(0),spectatorA(0),spectatorZ(0), 59 projectile3dNucleus(0),target3dNucleus(0) 61 projectile3dNucleus(0),target3dNucleus(0) 60 { 62 { 61 if(!ptr) { 63 if(!ptr) { 62 G4HadronicInteraction* p = 64 G4HadronicInteraction* p = 63 G4HadronicInteractionRegistry::Instance( 65 G4HadronicInteractionRegistry::Instance()->FindModel("PRECO"); 64 G4VPreCompoundModel* pre = static_cast<G4V 66 G4VPreCompoundModel* pre = static_cast<G4VPreCompoundModel*>(p); 65 if(!pre) { pre = new G4PreCompoundModel(); 67 if(!pre) { pre = new G4PreCompoundModel(); } 66 theProjectileFragmentation = pre; 68 theProjectileFragmentation = pre; 67 } 69 } 68 theModel = new G4BinaryCascade(theProjectile 70 theModel = new G4BinaryCascade(theProjectileFragmentation); 69 theHandler = theProjectileFragmentation->Get 71 theHandler = theProjectileFragmentation->GetExcitationHandler(); 70 theBLIR_ID = G4PhysicsModelCatalog::GetM << 72 if ( theBLIR_ID == -1 ) { 71 debug_G4BinaryLightIonReactionResults = G4Ha << 73 #ifdef G4MULTITHREADED >> 74 G4MUTEXLOCK(&G4BinaryLightIonReaction::BLIRMutex); >> 75 if ( theBLIR_ID == -1 ) { >> 76 #endif >> 77 theBLIR_ID = G4PhysicsModelCatalog::Register("Binary Light Ion Reaction"); >> 78 #ifdef G4MULTITHREADED >> 79 } >> 80 G4MUTEXUNLOCK(&G4BinaryLightIonReaction::BLIRMutex); >> 81 #endif >> 82 } >> 83 >> 84 >> 85 debug_G4BinaryLightIonReactionResults=getenv("debug_G4BinaryLightIonReactionResults")!=0; 72 } 86 } 73 87 74 G4BinaryLightIonReaction::~G4BinaryLightIonRea 88 G4BinaryLightIonReaction::~G4BinaryLightIonReaction() 75 {} 89 {} 76 90 77 void G4BinaryLightIonReaction::ModelDescriptio 91 void G4BinaryLightIonReaction::ModelDescription(std::ostream& outFile) const 78 { 92 { 79 outFile << "G4Binary Light Ion Cascade is an 93 outFile << "G4Binary Light Ion Cascade is an intra-nuclear cascade model\n" 80 << "using G4BinaryCasacde to model the i 94 << "using G4BinaryCasacde to model the interaction of a light\n" 81 << "nucleus with a nucleus.\n" 95 << "nucleus with a nucleus.\n" 82 << "The lighter of the two nuclei is tre 96 << "The lighter of the two nuclei is treated like a set of projectiles\n" 83 << "which are transported simultaneously << 97 << "which are transported simultanously through the heavier nucleus.\n"; 84 } 98 } 85 99 86 //-------------------------------------------- 100 //-------------------------------------------------------------------------------- 87 struct ReactionProduct4Mom 101 struct ReactionProduct4Mom 88 { 102 { 89 G4LorentzVector operator()(G4LorentzVector 103 G4LorentzVector operator()(G4LorentzVector a,G4ReactionProduct* b) {return a + G4LorentzVector(b->GetMomentum(), b->GetTotalEnergy() );} 90 }; 104 }; 91 105 92 G4HadFinalState *G4BinaryLightIonReaction:: 106 G4HadFinalState *G4BinaryLightIonReaction:: 93 ApplyYourself(const G4HadProjectile &aTrack, G 107 ApplyYourself(const G4HadProjectile &aTrack, G4Nucleus & targetNucleus ) 94 { 108 { 95 if(debug_G4BinaryLightIonReactionResults) G4 << 109 if(getenv("BLICDEBUG") ) G4cerr << " ######### Binary Light Ion Reaction starts ######### " << G4endl; 96 G4ping debug("debug_G4BinaryLightIonReaction 110 G4ping debug("debug_G4BinaryLightIonReaction"); 97 pA=aTrack.GetDefinition()->GetBaryonNumber() 111 pA=aTrack.GetDefinition()->GetBaryonNumber(); 98 pZ=G4lrint(aTrack.GetDefinition()->GetPDGCha 112 pZ=G4lrint(aTrack.GetDefinition()->GetPDGCharge()/eplus); 99 tA=targetNucleus.GetA_asInt(); 113 tA=targetNucleus.GetA_asInt(); 100 tZ=targetNucleus.GetZ_asInt(); 114 tZ=targetNucleus.GetZ_asInt(); 101 G4double timePrimary = aTrack.GetGlobalTime( 115 G4double timePrimary = aTrack.GetGlobalTime(); 102 G4LorentzVector mom(aTrack.Get4Momentum()); 116 G4LorentzVector mom(aTrack.Get4Momentum()); 103 //G4cout << "proj mom : " << mom << G4endl; 117 //G4cout << "proj mom : " << mom << G4endl; 104 G4LorentzRotation toBreit(mom.boostVector()) 118 G4LorentzRotation toBreit(mom.boostVector()); 105 119 106 G4bool swapped=SetLighterAsProjectile(mom, t 120 G4bool swapped=SetLighterAsProjectile(mom, toBreit); 107 //G4cout << "after swap, swapped? / mom " < 121 //G4cout << "after swap, swapped? / mom " << swapped << " / " << mom <<G4endl; 108 G4ReactionProductVector * result = 0; 122 G4ReactionProductVector * result = 0; 109 G4ReactionProductVector * cascaders=0; //new 123 G4ReactionProductVector * cascaders=0; //new G4ReactionProductVector; 110 // G4double m_nucl(0); // to check energ 124 // G4double m_nucl(0); // to check energy balance 111 125 112 126 113 // G4double m1=G4ParticleTable::GetPartic 127 // G4double m1=G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(pZ,pA); 114 // G4cout << "Entering the decision point 128 // G4cout << "Entering the decision point " 115 // << (mom.t()-mom.mag())/pA << " 129 // << (mom.t()-mom.mag())/pA << " " 116 // << pA<<" "<< pZ<<" " 130 // << pA<<" "<< pZ<<" " 117 // << tA<<" "<< tZ<<G4endl 131 // << tA<<" "<< tZ<<G4endl 118 // << " "<<mom.t()-mom.mag()<<" " 132 // << " "<<mom.t()-mom.mag()<<" " 119 // << mom.t()- m1<<G4endl; 133 // << mom.t()- m1<<G4endl; 120 if( (mom.t()-mom.mag())/pA < 50*MeV ) 134 if( (mom.t()-mom.mag())/pA < 50*MeV ) 121 { 135 { 122 // G4cout << "Using pre-compound only 136 // G4cout << "Using pre-compound only, E= "<<mom.t()-mom.mag()<<G4endl; 123 // m_nucl = mom.mag(); 137 // m_nucl = mom.mag(); 124 cascaders=FuseNucleiAndPrompound(mom); << 138 cascaders=FuseNucleiAndPrompound(mom); 125 if( !cascaders ) << 139 if( !cascaders ) 126 { << 140 { 127 << 141 128 // abort!! happens for too low e << 142 // abort!! happens for too low energy for nuclei to fuse 129 << 143 130 theResult.Clear(); << 144 theResult.Clear(); 131 theResult.SetStatusChange(isAliv << 145 theResult.SetStatusChange(isAlive); 132 theResult.SetEnergyChange(aTrack << 146 theResult.SetEnergyChange(aTrack.GetKineticEnergy()); 133 theResult.SetMomentumChange(aTra << 147 theResult.SetMomentumChange(aTrack.Get4Momentum().vect().unit()); 134 return &theResult; << 148 return &theResult; 135 } << 149 } 136 } 150 } 137 else 151 else 138 { 152 { 139 result=Interact(mom,toBreit); 153 result=Interact(mom,toBreit); 140 154 141 if(! result ) << 155 if(! result ) 142 { << 156 { 143 // abort!! << 157 // abort!! 144 << 158 145 G4cerr << "G4BinaryLightIonReaction << 159 G4cerr << "G4BinaryLightIonReaction no final state for: " << G4endl; 146 G4cerr << " Primary " << aTrack.Get << 160 G4cerr << " Primary " << aTrack.GetDefinition() 147 << ", (A,Z)=(" << aTrack.GetDefi << 161 << ", (A,Z)=(" << aTrack.GetDefinition()->GetBaryonNumber() 148 << "," << aTrack.GetDefinition() << 162 << "," << aTrack.GetDefinition()->GetPDGCharge()/eplus << ") " 149 << ", kinetic energy " << aTrack << 163 << ", kinetic energy " << aTrack.GetKineticEnergy() 150 << G4endl; << 164 << G4endl; 151 G4cerr << " Target nucleus (A,Z)=(" << 165 G4cerr << " Target nucleus (A,Z)=(" 152 << (swapped?pA:tA) << "," << 166 << (swapped?pA:tA) << "," 153 << (swapped?pZ:tZ) << ")" << << 167 << (swapped?pZ:tZ) << ")" << G4endl; 154 G4cerr << " if frequent, please sub << 168 G4cerr << " if frequent, please submit above information as bug report" 155 << G4endl << G4endl; << 169 << G4endl << G4endl; 156 << 157 theResult.Clear(); << 158 theResult.SetStatusChange(isAlive); << 159 theResult.SetEnergyChange(aTrack.Ge << 160 theResult.SetMomentumChange(aTrack. << 161 return &theResult; << 162 } << 163 << 164 // Calculate excitation energy, << 165 G4double theStatisticalExEnergy = GetProj << 166 << 167 << 168 pInitialState = mom; << 169 //G4cout << "BLIC: pInitialState from << 170 pInitialState.setT(pInitialState.getT() + << 171 G4ParticleTable::GetParticleTable()->GetIo << 172 //G4cout << "BLIC: target nucleus added << 173 << 174 delete target3dNucleus;target3dNucleus=0; << 175 delete projectile3dNucleus;projectile3dNu << 176 << 177 G4ReactionProductVector * spectators= new << 178 << 179 cascaders = new G4ReactionProductVector; << 180 170 181 G4LorentzVector pspectators=SortResult(re << 171 theResult.Clear(); 182 // this also sets spectatorA and << 172 theResult.SetStatusChange(isAlive); >> 173 theResult.SetEnergyChange(aTrack.GetKineticEnergy()); >> 174 theResult.SetMomentumChange(aTrack.Get4Momentum().vect().unit()); >> 175 return &theResult; >> 176 } 183 177 184 // pFinalState=std::accumulate(casca << 178 // Calculate excitation energy, >> 179 G4double theStatisticalExEnergy = GetProjectileExcitation(); 185 180 186 std::vector<G4ReactionProduct *>::iterato << 187 181 188 // G4cout << "pInitialState, pFin << 182 pInitialState = mom; >> 183 //G4cout << "BLIC: pInitialState from aTrack : " << pInitialState; >> 184 pInitialState.setT(pInitialState.getT() + >> 185 G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(tZ,tA)); >> 186 //G4cout << "BLIC: target nucleus added : " << pInitialState << G4endl; >> 187 >> 188 delete target3dNucleus;target3dNucleus=0; >> 189 delete projectile3dNucleus;projectile3dNucleus=0; >> 190 >> 191 G4ReactionProductVector * spectators= new G4ReactionProductVector; >> 192 >> 193 cascaders = new G4ReactionProductVector; >> 194 >> 195 G4LorentzVector pspectators=SortResult(result,spectators,cascaders); >> 196 // this also sets spectatorA and spectatorZ >> 197 >> 198 // pFinalState=std::accumulate(cascaders->begin(),cascaders->end(),pFinalState,ReactionProduct4Mom); >> 199 >> 200 std::vector<G4ReactionProduct *>::iterator iter; >> 201 >> 202 // G4cout << "pInitialState, pFinalState / pspectators"<< pInitialState << " / " << pFinalState << " / " << pspectators << G4endl; 189 // if ( spectA-spectatorA !=0 || spec 203 // if ( spectA-spectatorA !=0 || spectZ-spectatorZ !=0) 190 // { 204 // { 191 // G4cout << "spect Nucl != spect 205 // G4cout << "spect Nucl != spectators: nucl a,z; spect a,z" << 192 // spectatorA <<" "<< spectatorZ << 206 // spectatorA <<" "<< spectatorZ <<" ; " << spectA <<" "<< spectZ << G4endl; 193 // } 207 // } 194 delete result; << 208 delete result; 195 result=0; << 209 result=0; 196 G4LorentzVector momentum(pInitialState-pF << 210 G4LorentzVector momentum(pInitialState-pFinalState); 197 G4int loopcount(0); << 211 G4int loopcount(0); 198 //G4cout << "BLIC: momentum, pspectato << 212 //G4cout << "BLIC: momentum, pspectators : " << momentum << " / " << pspectators << G4endl; 199 while (std::abs(momentum.e()-pspectators. << 213 while (std::abs(momentum.e()-pspectators.e()) > 10*MeV) /* Loop checking, 31.08.2015, G.Folger */ 200 << 214 // see if on loopcount 201 { << 202 G4LorentzVector pCorrect(pInitialState- << 203 //G4cout << "BLIC:: BIC nonconservatio << 204 // Correct outgoing casacde particles.. << 205 G4bool EnergyIsCorrect=EnergyAndMomentu << 206 if ( ! EnergyIsCorrect && debug_G4Binar << 207 { << 208 G4cout << "Warning - G4BinaryLightIon << 209 } << 210 pFinalState=G4LorentzVector(0,0,0,0); << 211 for(iter=cascaders->begin(); iter!=casc << 212 { << 213 pFinalState += G4LorentzVector( (*ite << 214 } << 215 momentum=pInitialState-pFinalState; << 216 if (++loopcount > 10 ) << 217 { << 218 break; << 219 } << 220 } << 221 << 222 // Check if Energy/Momemtum is now ok, if << 223 if ( std::abs(momentum.e()-pspectators.e( << 224 { << 225 for (iter=spectators->begin();iter!=specta << 226 { 215 { 227 delete *iter; << 216 G4LorentzVector pCorrect(pInitialState-pspectators); >> 217 //G4cout << "BLIC:: BIC nonconservation? (pInitialState-pFinalState) / spectators :" << momentum << " / " << pspectators << "pCorrect "<< pCorrect<< G4endl; >> 218 // Correct outgoing casacde particles.... to have momentum of (initial state - spectators) >> 219 G4bool EnergyIsCorrect=EnergyAndMomentumCorrector(cascaders, pCorrect); >> 220 if ( ! EnergyIsCorrect && debug_G4BinaryLightIonReactionResults) >> 221 { >> 222 G4cout << "Warning - G4BinaryLightIonReaction E/P correction for cascaders failed" << G4endl; >> 223 } >> 224 pFinalState=G4LorentzVector(0,0,0,0); >> 225 for(iter=cascaders->begin(); iter!=cascaders->end(); iter++) >> 226 { >> 227 pFinalState += G4LorentzVector( (*iter)->GetMomentum(), (*iter)->GetTotalEnergy() ); >> 228 } >> 229 momentum=pInitialState-pFinalState; >> 230 if (++loopcount > 10 ) >> 231 { >> 232 if ( momentum.vect().mag() - momentum.e()> 10*keV ) >> 233 { >> 234 G4cerr << "G4BinaryLightIonReaction.cc: Cannot correct 4-momentum of cascade particles" << G4endl; >> 235 throw G4HadronicException(__FILE__, __LINE__, "G4BinaryCasacde::ApplyCollision()"); >> 236 } else { >> 237 break; >> 238 } >> 239 >> 240 } 228 } 241 } 229 delete spectators; << 242 230 for(iter=cascaders->begin(); iter!=cascad << 243 if (spectatorA > 0 ) 231 { << 244 { 232 delete *iter; << 245 // check spectator momentum 233 } << 246 if ( momentum.vect().mag() - momentum.e()< 10*keV ) 234 delete cascaders; << 247 { 235 << 248 // DeExciteSpectatorNucleus() also handles also case of A=1, Z=0,1 236 G4cout << "G4BinaryLightIonReaction.cc: m << 249 DeExciteSpectatorNucleus(spectators, cascaders, theStatisticalExEnergy, momentum); 237 << " initial - final " << momentum << 250 238 << " .. pInitialState/pFinalState/s << 251 } else { // momentum non-conservation --> fail 239 << pInitialState << G4endl << 252 for (iter=spectators->begin();iter!=spectators->end();iter++) 240 << pFinalState << G4endl << 253 { 241 << pspectators << G4endl << 254 delete *iter; 242 << " .. A,Z " << spectatorA <<" "<< << 255 } 243 G4cout << "G4BinaryLightIonReaction inval << 256 delete spectators; 244 G4cout << " Primary " << aTrack.GetDefini << 257 for(iter=cascaders->begin(); iter!=cascaders->end(); iter++) 245 << ", (A,Z)=(" << aTrack.GetDefin << 258 { 246 << "," << aTrack.GetDefinition()- << 259 delete *iter; 247 << ", kinetic energy " << aTrack. << 260 } 248 << G4endl; << 261 delete cascaders; 249 G4cout << " Target nucleus (A,Z)=(" << t << 262 250 << "," << targetNucleus.GetZ_ << 263 G4cout << "G4BinaryLightIonReaction.cc: mom check: " << momentum 251 G4cout << " if frequent, please submit ab << 264 << " 3.mag "<< momentum.vect().mag() << G4endl 252 << G4endl << G4endl; << 265 << " .. pInitialState/pFinalState/spectators " << pInitialState <<" " >> 266 << pFinalState << " " << pspectators << G4endl >> 267 << " .. A,Z " << spectatorA <<" "<< spectatorZ << G4endl; >> 268 G4cout << "G4BinaryLightIonReaction invalid final state for: " << G4endl; >> 269 G4cout << " Primary " << aTrack.GetDefinition() >> 270 << ", (A,Z)=(" << aTrack.GetDefinition()->GetBaryonNumber() >> 271 << "," << aTrack.GetDefinition()->GetPDGCharge()/eplus << ") " >> 272 << ", kinetic energy " << aTrack.GetKineticEnergy() >> 273 << G4endl; >> 274 G4cout << " Target nucleus (A,Z)=(" << targetNucleus.GetA_asInt() >> 275 << "," << targetNucleus.GetZ_asInt() << ")" << G4endl; >> 276 G4cout << " if frequent, please submit above information as bug report" >> 277 << G4endl << G4endl; 253 #ifdef debug_G4BinaryLightIonReaction 278 #ifdef debug_G4BinaryLightIonReaction 254 G4ExceptionDescription ed; 279 G4ExceptionDescription ed; 255 ed << "G4BinaryLightIonreaction: Ter 280 ed << "G4BinaryLightIonreaction: Terminate for above error" << G4endl; 256 G4Exception("G4BinaryLightIonreactio 281 G4Exception("G4BinaryLightIonreaction::ApplyYourSelf()", "BLIC001", FatalException, 257 ed); 282 ed); 258 283 259 #endif 284 #endif 260 theResult.Clear(); << 285 theResult.Clear(); 261 theResult.SetStatusChange(isAlive); << 286 theResult.SetStatusChange(isAlive); 262 theResult.SetEnergyChange(aTrack.GetKinet << 287 theResult.SetEnergyChange(aTrack.GetKineticEnergy()); 263 theResult.SetMomentumChange(aTrack.Get4Mo << 288 theResult.SetMomentumChange(aTrack.Get4Momentum().vect().unit()); 264 return &theResult; << 289 return &theResult; 265 << 290 } 266 } << 291 } else { // no spectators 267 if (spectatorA > 0 ) << 292 delete spectators; 268 { << 293 } 269 // DeExciteSpectatorNucleus() also << 270 DeExciteSpectatorNucleus(specta << 271 } else { // no spectators << 272 delete spectators; << 273 } << 274 } 294 } 275 // Rotate to lab 295 // Rotate to lab 276 G4LorentzRotation toZ; 296 G4LorentzRotation toZ; 277 toZ.rotateZ(-1*mom.phi()); 297 toZ.rotateZ(-1*mom.phi()); 278 toZ.rotateY(-1*mom.theta()); 298 toZ.rotateY(-1*mom.theta()); 279 G4LorentzRotation toLab(toZ.inverse()); 299 G4LorentzRotation toLab(toZ.inverse()); 280 300 281 // Fill the particle change, while rotating. 301 // Fill the particle change, while rotating. Boost from projectile breit-frame in case we swapped. 282 // theResult.Clear(); 302 // theResult.Clear(); 283 theResult.Clear(); 303 theResult.Clear(); 284 theResult.SetStatusChange(stopAndKill); 304 theResult.SetStatusChange(stopAndKill); 285 G4LorentzVector ptot(0); 305 G4LorentzVector ptot(0); 286 #ifdef debug_BLIR_result << 306 G4ReactionProductVector::iterator iter; 287 G4LorentzVector p_raw; << 307 #ifdef debug_BLIR_result 288 #endif << 308 G4LorentzVector p_raw; 289 //G4int i=0; << 309 #endif >> 310 //G4int i=0; 290 311 291 G4ReactionProductVector::iterator iter << 292 for(iter=cascaders->begin(); iter!=cascaders 312 for(iter=cascaders->begin(); iter!=cascaders->end(); iter++) 293 { 313 { 294 if((*iter)->GetNewlyAdded()) 314 if((*iter)->GetNewlyAdded()) 295 { 315 { 296 G4DynamicParticle * aNewDP = 316 G4DynamicParticle * aNewDP = 297 new G4DynamicParticle((*iter)->GetDe 317 new G4DynamicParticle((*iter)->GetDefinition(), 298 (*iter)->GetTotalEnergy(), 318 (*iter)->GetTotalEnergy(), 299 (*iter)->GetMomentum() ); 319 (*iter)->GetMomentum() ); 300 G4LorentzVector tmp = aNewDP->Get4Moment 320 G4LorentzVector tmp = aNewDP->Get4Momentum(); 301 #ifdef debug_BLIR_result 321 #ifdef debug_BLIR_result 302 p_raw+= tmp; 322 p_raw+= tmp; 303 #endif 323 #endif 304 if(swapped) 324 if(swapped) 305 { 325 { 306 tmp*=toBreit.inverse(); 326 tmp*=toBreit.inverse(); 307 tmp.setVect(-tmp.vect()); 327 tmp.setVect(-tmp.vect()); 308 } 328 } 309 tmp *= toLab; 329 tmp *= toLab; 310 aNewDP->Set4Momentum(tmp); 330 aNewDP->Set4Momentum(tmp); 311 G4HadSecondary aNew = G4HadSecondary(aNe 331 G4HadSecondary aNew = G4HadSecondary(aNewDP); 312 G4double time = 0; 332 G4double time = 0; //(*iter)->GetCreationTime(); 313 //if(time < 0.0) { time = 0.0; } 333 //if(time < 0.0) { time = 0.0; } 314 aNew.SetTime(timePrimary + time); 334 aNew.SetTime(timePrimary + time); 315 //aNew.SetCreatorModelID((*iter)-> << 335 aNew.SetCreatorModelType((*iter)->GetCreatorModel()); 316 aNew.SetCreatorModelID(theBLIR_ID) << 317 336 318 theResult.AddSecondary(aNew); 337 theResult.AddSecondary(aNew); 319 ptot += tmp; 338 ptot += tmp; 320 //G4cout << "BLIC: Secondary " < 339 //G4cout << "BLIC: Secondary " << aNew->GetDefinition()->GetParticleName() 321 // <<" "<< aNew->GetMomen 340 // <<" "<< aNew->GetMomentum()<<" "<< aNew->GetTotalEnergy() << G4endl; 322 } 341 } 323 delete *iter; 342 delete *iter; 324 } 343 } 325 delete cascaders; 344 delete cascaders; 326 345 327 #ifdef debug_BLIR_result 346 #ifdef debug_BLIR_result 328 //G4cout << "Result analysis, secondaries " 347 //G4cout << "Result analysis, secondaries " << theResult.GetNumberOfSecondaries() << G4endl; 329 //G4cout << "p_tot_raw " << p_raw << " sum p 348 //G4cout << "p_tot_raw " << p_raw << " sum p final " << ptot << G4endl; 330 G4double m_nucl= G4ParticleTable::GetPartic 349 G4double m_nucl= G4ParticleTable::GetParticleTable()->GetIonTable()-> 331 GetIonMass(targetNucleus.GetZ_asInt( 350 GetIonMass(targetNucleus.GetZ_asInt(),targetNucleus.GetA_asInt()); 332 // delete? tZ=targetNucleus.GetZ_asInt(); 351 // delete? tZ=targetNucleus.GetZ_asInt(); 333 352 334 //G4cout << "BLIC Energy conservation initia 353 //G4cout << "BLIC Energy conservation initial/primary/nucleus/final/delta(init-final) " 335 // << aTrack.GetTotalEnergy() + m_nuc 354 // << aTrack.GetTotalEnergy() + m_nucl <<" "<< aTrack.GetTotalEnergy() <<" "<< m_nucl <<" "<<ptot.e() 336 // <<" "<< aTrack.GetTotalEnergy() + m_ 355 // <<" "<< aTrack.GetTotalEnergy() + m_nucl - ptot.e() << G4endl; 337 G4cout << "BLIC momentum conservation " << a 356 G4cout << "BLIC momentum conservation " << aTrack.Get4Momentum()+ G4LorentzVector(m_nucl) 338 << " ptot " << ptot << " delta " << aTra 357 << " ptot " << ptot << " delta " << aTrack.Get4Momentum()+ G4LorentzVector(m_nucl) - ptot 339 << " 3mom.mag() " << (aTrack.Get4 358 << " 3mom.mag() " << (aTrack.Get4Momentum()+ G4LorentzVector(m_nucl) - ptot).vect().mag() << G4endl; 340 #endif 359 #endif 341 360 342 if(debug_G4BinaryLightIonReactionResults) G4 << 361 if(getenv("BLICDEBUG") ) G4cerr << " ######### Binary Light Ion Reaction number ends ######### " << G4endl; 343 362 344 return &theResult; 363 return &theResult; 345 } 364 } 346 365 347 //-------------------------------------------- 366 //-------------------------------------------------------------------------------- 348 367 349 //******************************************** 368 //**************************************************************************** 350 G4bool G4BinaryLightIonReaction::EnergyAndMome 369 G4bool G4BinaryLightIonReaction::EnergyAndMomentumCorrector( 351 G4ReactionProductVector* Output, G4Lorentz 370 G4ReactionProductVector* Output, G4LorentzVector& TotalCollisionMom) 352 //******************************************** 371 //**************************************************************************** 353 { 372 { 354 const int nAttemptScale = 2500; 373 const int nAttemptScale = 2500; 355 const double ErrLimit = 1.E-6; 374 const double ErrLimit = 1.E-6; 356 if (Output->empty()) 375 if (Output->empty()) 357 return TRUE; 376 return TRUE; 358 G4LorentzVector SumMom(0,0,0,0); 377 G4LorentzVector SumMom(0,0,0,0); 359 G4double SumMass = 0; 378 G4double SumMass = 0; 360 G4double TotalCollisionMass = TotalCo 379 G4double TotalCollisionMass = TotalCollisionMom.m(); 361 size_t i = 0; 380 size_t i = 0; 362 // Calculate sum hadron 4-momenta and summin 381 // Calculate sum hadron 4-momenta and summing hadron mass 363 for(i = 0; i < Output->size(); i++) 382 for(i = 0; i < Output->size(); i++) 364 { 383 { 365 SumMom += G4LorentzVector((*Output)[i]->G 384 SumMom += G4LorentzVector((*Output)[i]->GetMomentum(),(*Output)[i]->GetTotalEnergy()); 366 SumMass += (*Output)[i]->GetDefinition()-> 385 SumMass += (*Output)[i]->GetDefinition()->GetPDGMass(); 367 } 386 } 368 // G4cout << " E/P corrector, SumMass, Sum 387 // G4cout << " E/P corrector, SumMass, SumMom.m2, TotalMass " 369 // << SumMass <<" "<< SumMom.m2() << 388 // << SumMass <<" "<< SumMom.m2() <<" "<<TotalCollisionMass<< G4endl; 370 if (SumMass > TotalCollisionMass) return FAL 389 if (SumMass > TotalCollisionMass) return FALSE; 371 SumMass = SumMom.m2(); 390 SumMass = SumMom.m2(); 372 if (SumMass < 0) return FALSE; 391 if (SumMass < 0) return FALSE; 373 SumMass = std::sqrt(SumMass); 392 SumMass = std::sqrt(SumMass); 374 393 375 // Compute c.m.s. hadron velocity and boost 394 // Compute c.m.s. hadron velocity and boost KTV to hadron c.m.s. 376 G4ThreeVector Beta = -SumMom.boostVector(); 395 G4ThreeVector Beta = -SumMom.boostVector(); 377 //G4cout << " == pre boost 2 "<< SumMo 396 //G4cout << " == pre boost 2 "<< SumMom.e()<< " "<< SumMom.mag()<<" "<< Beta <<G4endl; 378 //--old Output->Boost(Beta); 397 //--old Output->Boost(Beta); 379 for(i = 0; i < Output->size(); i++) 398 for(i = 0; i < Output->size(); i++) 380 { 399 { 381 G4LorentzVector mom = G4LorentzVector((*Ou 400 G4LorentzVector mom = G4LorentzVector((*Output)[i]->GetMomentum(),(*Output)[i]->GetTotalEnergy()); 382 mom *= Beta; 401 mom *= Beta; 383 (*Output)[i]->SetMomentum(mom.vect()); 402 (*Output)[i]->SetMomentum(mom.vect()); 384 (*Output)[i]->SetTotalEnergy(mom.e()); 403 (*Output)[i]->SetTotalEnergy(mom.e()); 385 } 404 } 386 405 387 // Scale total c.m.s. hadron energy (hadron 406 // Scale total c.m.s. hadron energy (hadron system mass). 388 // It should be equal interaction mass 407 // It should be equal interaction mass 389 G4double Scale = 0,OldScale=0; 408 G4double Scale = 0,OldScale=0; 390 G4double factor = 1.; 409 G4double factor = 1.; 391 G4int cAttempt = 0; 410 G4int cAttempt = 0; 392 G4double Sum = 0; 411 G4double Sum = 0; 393 G4bool success = false; 412 G4bool success = false; 394 for(cAttempt = 0; cAttempt < nAttemptScale; 413 for(cAttempt = 0; cAttempt < nAttemptScale; cAttempt++) 395 { 414 { 396 Sum = 0; 415 Sum = 0; 397 for(i = 0; i < Output->size(); i++) 416 for(i = 0; i < Output->size(); i++) 398 { 417 { 399 G4LorentzVector HadronMom = G4LorentzVec 418 G4LorentzVector HadronMom = G4LorentzVector((*Output)[i]->GetMomentum(),(*Output)[i]->GetTotalEnergy()); 400 HadronMom.setVect(HadronMom.vect()+ fact 419 HadronMom.setVect(HadronMom.vect()+ factor*Scale*HadronMom.vect()); 401 G4double E = std::sqrt(HadronMom.vect(). 420 G4double E = std::sqrt(HadronMom.vect().mag2() + sqr((*Output)[i]->GetDefinition()->GetPDGMass())); 402 HadronMom.setE(E); 421 HadronMom.setE(E); 403 (*Output)[i]->SetMomentum(HadronMom.vect 422 (*Output)[i]->SetMomentum(HadronMom.vect()); 404 (*Output)[i]->SetTotalEnergy(HadronMom.e 423 (*Output)[i]->SetTotalEnergy(HadronMom.e()); 405 Sum += E; 424 Sum += E; 406 } 425 } 407 OldScale=Scale; 426 OldScale=Scale; 408 Scale = TotalCollisionMass/Sum - 1; 427 Scale = TotalCollisionMass/Sum - 1; 409 // G4cout << "E/P corr - " << cAttempt << 428 // G4cout << "E/P corr - " << cAttempt << " " << Scale << G4endl; 410 if (std::abs(Scale) <= ErrLimit 429 if (std::abs(Scale) <= ErrLimit 411 || OldScale == Scale) // protect ' 430 || OldScale == Scale) // protect 'frozen' situation and divide by 0 in calculating new factor below 412 { 431 { 413 if (debug_G4BinaryLightIonReactionResult 432 if (debug_G4BinaryLightIonReactionResults) G4cout << "E/p corrector: " << cAttempt << G4endl; 414 success = true; 433 success = true; 415 break; 434 break; 416 } 435 } 417 if ( cAttempt > 10 ) 436 if ( cAttempt > 10 ) 418 { 437 { 419 // G4cout << " speed it up? " << 438 // G4cout << " speed it up? " << std::abs(OldScale/(OldScale-Scale)) << G4endl; 420 factor=std::max(1.,G4Log(std::abs(OldSca 439 factor=std::max(1.,G4Log(std::abs(OldScale/(OldScale-Scale)))); 421 // G4cout << " ? factor ? " << factor 440 // G4cout << " ? factor ? " << factor << G4endl; 422 } 441 } 423 } 442 } 424 443 425 if( (!success) && debug_G4BinaryLightIonRea 444 if( (!success) && debug_G4BinaryLightIonReactionResults) 426 { 445 { 427 G4cout << "G4G4BinaryLightIonReaction::Ene 446 G4cout << "G4G4BinaryLightIonReaction::EnergyAndMomentumCorrector - Warning"<<G4endl; 428 G4cout << " Scale not unity at end of it 447 G4cout << " Scale not unity at end of iteration loop: "<<TotalCollisionMass<<" "<<Sum<<" "<<Scale<<G4endl; 429 G4cout << " Increase number of attempts 448 G4cout << " Increase number of attempts or increase ERRLIMIT"<<G4endl; 430 } 449 } 431 450 432 // Compute c.m.s. interaction velocity and K 451 // Compute c.m.s. interaction velocity and KTV back boost 433 Beta = TotalCollisionMom.boostVector(); 452 Beta = TotalCollisionMom.boostVector(); 434 //--old Output->Boost(Beta); 453 //--old Output->Boost(Beta); 435 for(i = 0; i < Output->size(); i++) 454 for(i = 0; i < Output->size(); i++) 436 { 455 { 437 G4LorentzVector mom = G4LorentzVector((*Ou 456 G4LorentzVector mom = G4LorentzVector((*Output)[i]->GetMomentum(),(*Output)[i]->GetTotalEnergy()); 438 mom *= Beta; 457 mom *= Beta; 439 (*Output)[i]->SetMomentum(mom.vect()); 458 (*Output)[i]->SetMomentum(mom.vect()); 440 (*Output)[i]->SetTotalEnergy(mom.e()); 459 (*Output)[i]->SetTotalEnergy(mom.e()); 441 } 460 } 442 return TRUE; 461 return TRUE; 443 } 462 } 444 G4bool G4BinaryLightIonReaction::SetLighterAsP 463 G4bool G4BinaryLightIonReaction::SetLighterAsProjectile(G4LorentzVector & mom,const G4LorentzRotation & toBreit) 445 { 464 { 446 G4bool swapped = false; 465 G4bool swapped = false; 447 if(tA<pA) 466 if(tA<pA) 448 { 467 { 449 swapped = true; 468 swapped = true; 450 G4int tmp(0); 469 G4int tmp(0); 451 tmp = tA; tA=pA; pA=tmp; 470 tmp = tA; tA=pA; pA=tmp; 452 tmp = tZ; tZ=pZ; pZ=tmp; 471 tmp = tZ; tZ=pZ; pZ=tmp; 453 G4double m1=G4ParticleTable::GetParticle 472 G4double m1=G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(pZ,pA); 454 G4LorentzVector it(m1, G4ThreeVector(0,0 473 G4LorentzVector it(m1, G4ThreeVector(0,0,0)); 455 mom = toBreit*it; 474 mom = toBreit*it; 456 } 475 } 457 return swapped; 476 return swapped; 458 } 477 } 459 G4ReactionProductVector * G4BinaryLightIonReac 478 G4ReactionProductVector * G4BinaryLightIonReaction::FuseNucleiAndPrompound(const G4LorentzVector & mom) 460 { 479 { 461 // Check if kinematically nuclei can fuse. 480 // Check if kinematically nuclei can fuse. 462 G4double mFused=G4ParticleTable::GetParticl 481 G4double mFused=G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(pZ+tZ,pA+tA); 463 G4double mTarget=G4ParticleTable::GetPartic 482 G4double mTarget=G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(tZ,tA); 464 G4LorentzVector pCompound(mom.e()+mTarget,m 483 G4LorentzVector pCompound(mom.e()+mTarget,mom.vect()); 465 G4double m2Compound=pCompound.m2(); 484 G4double m2Compound=pCompound.m2(); 466 if (m2Compound < sqr(mFused) ) { 485 if (m2Compound < sqr(mFused) ) { 467 //G4cout << "G4BLIC: projectile p, mTarge 486 //G4cout << "G4BLIC: projectile p, mTarget, mFused, mCompound, delta: " <<mom << " " << mTarget << " " << mFused 468 // << " " << sqrt(m2Compound)<< " " 487 // << " " << sqrt(m2Compound)<< " " << sqrt(m2Compound) - mFused << G4endl; 469 return 0; 488 return 0; 470 } 489 } 471 490 472 G4Fragment aPreFrag; 491 G4Fragment aPreFrag; 473 aPreFrag.SetZandA_asInt(pZ+tZ, pA+tA); 492 aPreFrag.SetZandA_asInt(pZ+tZ, pA+tA); 474 aPreFrag.SetNumberOfParticles(pA); 493 aPreFrag.SetNumberOfParticles(pA); 475 aPreFrag.SetNumberOfCharged(pZ); 494 aPreFrag.SetNumberOfCharged(pZ); 476 aPreFrag.SetNumberOfHoles(0); 495 aPreFrag.SetNumberOfHoles(0); 477 //GF FIXME: whyusing plop in z direction? t 496 //GF FIXME: whyusing plop in z direction? this will not conserve momentum? 478 //G4ThreeVector plop(0.,0., mom.vect().mag( 497 //G4ThreeVector plop(0.,0., mom.vect().mag()); 479 //G4LorentzVector aL(mom.t()+mTarget, plop) 498 //G4LorentzVector aL(mom.t()+mTarget, plop); 480 G4LorentzVector aL(mom.t()+mTarget,mom.vect 499 G4LorentzVector aL(mom.t()+mTarget,mom.vect()); 481 aPreFrag.SetMomentum(aL); 500 aPreFrag.SetMomentum(aL); 482 501 483 502 484 //G4cout << "Fragment INFO "<< pA+tA 503 //G4cout << "Fragment INFO "<< pA+tA <<" "<<pZ+tZ<<" " 485 // << aL <<" "<<G4endl << aPreF 504 // << aL <<" "<<G4endl << aPreFrag << G4endl; 486 G4ReactionProductVector * cascaders = thePr 505 G4ReactionProductVector * cascaders = theProjectileFragmentation->DeExcite(aPreFrag); 487 //G4double tSum = 0; 506 //G4double tSum = 0; 488 for(size_t count = 0; count<cascaders->size 507 for(size_t count = 0; count<cascaders->size(); count++) 489 { 508 { 490 cascaders->operator[](count)->SetNewlyAd 509 cascaders->operator[](count)->SetNewlyAdded(true); 491 //tSum += cascaders->operator[](count)-> 510 //tSum += cascaders->operator[](count)->GetKineticEnergy(); 492 } 511 } 493 // G4cout << "Exiting pre-compound on 512 // G4cout << "Exiting pre-compound only, E= "<<tSum<<G4endl; 494 return cascaders; 513 return cascaders; 495 } 514 } 496 G4ReactionProductVector * G4BinaryLightIonReac 515 G4ReactionProductVector * G4BinaryLightIonReaction::Interact(G4LorentzVector & mom, const G4LorentzRotation & toBreit) 497 { 516 { 498 G4ReactionProductVector * result = 0; 517 G4ReactionProductVector * result = 0; 499 G4double projectileMass(0); 518 G4double projectileMass(0); 500 G4LorentzVector it; 519 G4LorentzVector it; 501 520 502 G4int tryCount(0); 521 G4int tryCount(0); 503 do 522 do 504 { 523 { 505 ++tryCount; 524 ++tryCount; 506 projectile3dNucleus = new G4Fancy3DNu 525 projectile3dNucleus = new G4Fancy3DNucleus; 507 projectile3dNucleus->Init(pA, pZ); 526 projectile3dNucleus->Init(pA, pZ); 508 projectile3dNucleus->CenterNucleons() 527 projectile3dNucleus->CenterNucleons(); 509 projectileMass=G4ParticleTable::GetPa 528 projectileMass=G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass( 510 projectile3dNucleus->GetCharge( 529 projectile3dNucleus->GetCharge(),projectile3dNucleus->GetMassNumber()); 511 it=toBreit * G4LorentzVector(projecti 530 it=toBreit * G4LorentzVector(projectileMass,G4ThreeVector(0,0,0)); 512 531 513 target3dNucleus = new G4Fancy3DNucleu 532 target3dNucleus = new G4Fancy3DNucleus; 514 target3dNucleus->Init(tA, tZ); 533 target3dNucleus->Init(tA, tZ); 515 G4double impactMax = target3dNucleus- 534 G4double impactMax = target3dNucleus->GetOuterRadius()+projectile3dNucleus->GetOuterRadius(); 516 // G4cout << "out radius - nuc 535 // G4cout << "out radius - nucleus - projectile " << target3dNucleus->GetOuterRadius()/fermi << " - " << projectile3dNucleus->GetOuterRadius()/fermi << G4endl; 517 G4double aX=(2.*G4UniformRand()-1.)*i 536 G4double aX=(2.*G4UniformRand()-1.)*impactMax; 518 G4double aY=(2.*G4UniformRand()-1.)*i 537 G4double aY=(2.*G4UniformRand()-1.)*impactMax; 519 G4ThreeVector pos(aX, aY, -2.*impactM 538 G4ThreeVector pos(aX, aY, -2.*impactMax-5.*fermi); 520 539 521 G4KineticTrackVector * initalState = 540 G4KineticTrackVector * initalState = new G4KineticTrackVector; 522 projectile3dNucleus->StartLoop(); 541 projectile3dNucleus->StartLoop(); 523 G4Nucleon * aNuc; 542 G4Nucleon * aNuc; 524 G4LorentzVector tmpV(0,0,0,0); 543 G4LorentzVector tmpV(0,0,0,0); 525 #ifdef debug_BLIR_finalstate 544 #ifdef debug_BLIR_finalstate 526 G4LorentzVector pinitial; 545 G4LorentzVector pinitial; 527 #endif 546 #endif 528 G4LorentzVector nucleonMom(1./pA*mom) 547 G4LorentzVector nucleonMom(1./pA*mom); 529 nucleonMom.setZ(nucleonMom.vect().mag 548 nucleonMom.setZ(nucleonMom.vect().mag()); 530 nucleonMom.setX(0); 549 nucleonMom.setX(0); 531 nucleonMom.setY(0); 550 nucleonMom.setY(0); 532 theFermi.Init(pA,pZ); 551 theFermi.Init(pA,pZ); 533 while( (aNuc=projectile3dNucleus->Get 552 while( (aNuc=projectile3dNucleus->GetNextNucleon()) ) /* Loop checking, 31.08.2015, G.Folger */ 534 { 553 { 535 G4LorentzVector p4 = aNuc->GetMome 554 G4LorentzVector p4 = aNuc->GetMomentum(); 536 tmpV+=p4; 555 tmpV+=p4; 537 G4ThreeVector nucleonPosition(aNuc 556 G4ThreeVector nucleonPosition(aNuc->GetPosition()); 538 G4double density=(projectile3dNucl 557 G4double density=(projectile3dNucleus->GetNuclearDensity())->GetDensity(nucleonPosition); 539 nucleonPosition += pos; 558 nucleonPosition += pos; 540 G4KineticTrack * it1 = new G4Kinet 559 G4KineticTrack * it1 = new G4KineticTrack(aNuc, nucleonPosition, nucleonMom ); 541 it1->SetState(G4KineticTrack::outs 560 it1->SetState(G4KineticTrack::outside); 542 G4double pfermi= theFermi.GetFermi 561 G4double pfermi= theFermi.GetFermiMomentum(density); 543 G4double mass = aNuc->GetDefinitio 562 G4double mass = aNuc->GetDefinition()->GetPDGMass(); 544 G4double Efermi= std::sqrt( sqr(ma 563 G4double Efermi= std::sqrt( sqr(mass) + sqr(pfermi)) - mass; 545 it1->SetProjectilePotential(-Eferm 564 it1->SetProjectilePotential(-Efermi); 546 initalState->push_back(it1); 565 initalState->push_back(it1); 547 #ifdef debug_BLIR_finalstate 566 #ifdef debug_BLIR_finalstate 548 pinitial += it1->Get4Momentum() 567 pinitial += it1->Get4Momentum(); 549 #endif 568 #endif 550 } 569 } 551 570 552 result=theModel->Propagate(initalStat 571 result=theModel->Propagate(initalState, target3dNucleus); 553 #ifdef debug_BLIR_finalstate 572 #ifdef debug_BLIR_finalstate 554 if( result && result->size()>0) << 573 if( result && result->size()>0) 555 { << 574 { 556 G4cout << " Cascade result " << G4endl; << 575 G4LorentzVector presult; 557 G4LorentzVector presult; << 576 G4ReactionProductVector::iterator iter; 558 G4ReactionProductVector::iterator << 577 G4ReactionProduct xp; 559 G4ReactionProduct xp; << 578 for (iter=result->begin(); iter !=result->end(); ++iter) 560 for (iter=result->begin(); iter ! << 579 { 561 { << 580 presult += G4LorentzVector((*iter)->GetMomentum(),(*iter)->GetTotalEnergy()); 562 presult += G4LorentzVector((*ite << 581 } 563 G4cout << (*iter)->GetDefinition()->Ge << 582 564 << "("<< (*iter)->GetMomentum().x()<<" << 583 G4cout << "BLIC check result : initial " << pinitial << " mass tgt " << target3dNucleus->GetMass() 565 << (*iter)->GetMomentum().y()<<"," << 584 << " final " << presult 566 << (*iter)->GetMomentum().z()<<";" << 585 << " IF - FF " << pinitial +G4LorentzVector(target3dNucleus->GetMass()) - presult << G4endl; 567 << (*iter)->GetTotalEnergy() <<")"< << 568 } << 569 << 570 G4cout << "BLIC check result : in << 571 << " final " << presult << 572 << " IF - FF " << pinitial +G << 573 586 574 } << 587 } 575 #endif 588 #endif 576 if( result && result->size()==0) 589 if( result && result->size()==0) 577 { 590 { 578 delete result; 591 delete result; 579 result=0; 592 result=0; 580 } 593 } 581 if ( ! result ) 594 if ( ! result ) 582 { 595 { 583 delete target3dNucleus; 596 delete target3dNucleus; 584 delete projectile3dNucleus; 597 delete projectile3dNucleus; 585 } 598 } 586 599 587 // std::for_each(initalState->begin() 600 // std::for_each(initalState->begin(), initalState->end(), Delete<G4KineticTrack>()); 588 // delete initalState; 601 // delete initalState; 589 602 590 } while (! result && tryCount< 150); / 603 } while (! result && tryCount< 150); /* Loop checking, 31.08.2015, G.Folger */ 591 return result; 604 return result; 592 } 605 } 593 G4double G4BinaryLightIonReaction::GetProjecti 606 G4double G4BinaryLightIonReaction::GetProjectileExcitation() 594 { 607 { 595 608 596 G4Nucleon * aNuc; 609 G4Nucleon * aNuc; 597 // the projectileNucleus excitation ener 610 // the projectileNucleus excitation energy estimate... 598 G4double theStatisticalExEnergy = 0; 611 G4double theStatisticalExEnergy = 0; 599 projectile3dNucleus->StartLoop(); 612 projectile3dNucleus->StartLoop(); 600 while( (aNuc=projectile3dNucleus->GetNex 613 while( (aNuc=projectile3dNucleus->GetNextNucleon()) ) /* Loop checking, 31.08.2015, G.Folger */ 601 { 614 { 602 //G4cout << " Nucleon : " << a 615 //G4cout << " Nucleon : " << aNuc->GetDefinition()->GetParticleName() <<" "<< aNuc->AreYouHit() <<" "<<aNuc->GetMomentum()<<G4endl; 603 if(aNuc->AreYouHit()) { 616 if(aNuc->AreYouHit()) { 604 G4ThreeVector aPosition(aNuc->GetP 617 G4ThreeVector aPosition(aNuc->GetPosition()); 605 G4double localDensity = projectile 618 G4double localDensity = projectile3dNucleus->GetNuclearDensity()->GetDensity(aPosition); 606 G4double localPfermi = theFermi.Ge 619 G4double localPfermi = theFermi.GetFermiMomentum(localDensity); 607 G4double nucMass = aNuc->GetDefini 620 G4double nucMass = aNuc->GetDefinition()->GetPDGMass(); 608 G4double localFermiEnergy = std::s 621 G4double localFermiEnergy = std::sqrt(nucMass*nucMass + localPfermi*localPfermi) - nucMass; 609 G4double deltaE = localFermiEnergy 622 G4double deltaE = localFermiEnergy - (aNuc->GetMomentum().t()-aNuc->GetMomentum().mag()); 610 theStatisticalExEnergy += deltaE; 623 theStatisticalExEnergy += deltaE; 611 } 624 } 612 } 625 } 613 return theStatisticalExEnergy; 626 return theStatisticalExEnergy; 614 } 627 } 615 628 616 G4LorentzVector G4BinaryLightIonReaction::Sort 629 G4LorentzVector G4BinaryLightIonReaction::SortResult(G4ReactionProductVector * result, G4ReactionProductVector * spectators,G4ReactionProductVector * cascaders) 617 { 630 { 618 unsigned int i(0); 631 unsigned int i(0); 619 spectatorA=spectatorZ=0; 632 spectatorA=spectatorZ=0; 620 G4LorentzVector pspectators(0,0,0,0); 633 G4LorentzVector pspectators(0,0,0,0); 621 pFinalState=G4LorentzVector(0,0,0,0); 634 pFinalState=G4LorentzVector(0,0,0,0); 622 for(i=0; i<result->size(); i++) 635 for(i=0; i<result->size(); i++) 623 { 636 { 624 if( (*result)[i]->GetNewlyAdded() ) 637 if( (*result)[i]->GetNewlyAdded() ) 625 { 638 { 626 pFinalState += G4LorentzVector( (*res 639 pFinalState += G4LorentzVector( (*result)[i]->GetMomentum(), (*result)[i]->GetTotalEnergy() ); 627 cascaders->push_back((*result)[i]); 640 cascaders->push_back((*result)[i]); 628 } 641 } 629 else { 642 else { 630 // G4cout <<" spectator ... 643 // G4cout <<" spectator ... "; 631 pspectators += G4LorentzVector( (*res 644 pspectators += G4LorentzVector( (*result)[i]->GetMomentum(), (*result)[i]->GetTotalEnergy() ); 632 spectators->push_back((*result)[i]); 645 spectators->push_back((*result)[i]); 633 spectatorA++; 646 spectatorA++; 634 spectatorZ+= G4lrint((*result)[i]->Ge 647 spectatorZ+= G4lrint((*result)[i]->GetDefinition()->GetPDGCharge()/eplus); 635 } 648 } 636 649 637 // G4cout << (*result)[i]<< " " 650 // G4cout << (*result)[i]<< " " 638 // << (*result)[i]->GetDefinition 651 // << (*result)[i]->GetDefinition()->GetParticleName() << " " 639 // << (*result)[i]->GetMomentum() 652 // << (*result)[i]->GetMomentum()<< " " 640 // << (*result)[i]->GetTotalEnerg 653 // << (*result)[i]->GetTotalEnergy() << G4endl; 641 } 654 } 642 //G4cout << "pFinalState / pspectators, 655 //G4cout << "pFinalState / pspectators, (A,Z), p " << pFinalState << " / " << spectators->size() 643 // << " (" << spectatorA << ", "<< sp 656 // << " (" << spectatorA << ", "<< spectatorZ << "), 4-mom: " << pspectators << G4endl; 644 657 645 return pspectators; 658 return pspectators; 646 } 659 } 647 660 648 void G4BinaryLightIonReaction::DeExciteSpectat 661 void G4BinaryLightIonReaction::DeExciteSpectatorNucleus(G4ReactionProductVector * spectators, G4ReactionProductVector * cascaders, 649 662 G4double theStatisticalExEnergy, G4LorentzVector & pSpectators) 650 { 663 { 651 // call precompound model 664 // call precompound model 652 G4ReactionProductVector * proFrag = 0; 665 G4ReactionProductVector * proFrag = 0; 653 G4LorentzVector pFragment(0.,0.,0.,0.); 666 G4LorentzVector pFragment(0.,0.,0.,0.); 654 // G4cout << " == pre boost 1 "<< mome 667 // G4cout << " == pre boost 1 "<< momentum.e()<< " "<< momentum.mag()<<G4endl; 655 G4LorentzRotation boost_fragments; 668 G4LorentzRotation boost_fragments; 656 // G4cout << " == post boost 1 "<< mom 669 // G4cout << " == post boost 1 "<< momentum.e()<< " "<< momentum.mag()<<G4endl; 657 // G4LorentzRotation boost_spectator_mom 670 // G4LorentzRotation boost_spectator_mom(-momentum.boostVector()); 658 // G4cout << "- momentum " << boost_spe 671 // G4cout << "- momentum " << boost_spectator_mom * momentum << G4endl; 659 G4LorentzVector pFragments(0,0,0,0); 672 G4LorentzVector pFragments(0,0,0,0); 660 673 661 if(spectatorZ>0 && spectatorA>1) 674 if(spectatorZ>0 && spectatorA>1) 662 { 675 { 663 // Make the fragment 676 // Make the fragment 664 G4Fragment aProRes; 677 G4Fragment aProRes; 665 aProRes.SetZandA_asInt(spectatorZ, spect 678 aProRes.SetZandA_asInt(spectatorZ, spectatorA); 666 aProRes.SetNumberOfParticles(0); 679 aProRes.SetNumberOfParticles(0); 667 aProRes.SetNumberOfCharged(0); 680 aProRes.SetNumberOfCharged(0); 668 aProRes.SetNumberOfHoles(pA-spectatorA); 681 aProRes.SetNumberOfHoles(pA-spectatorA); 669 G4double mFragment=G4ParticleTable::GetP 682 G4double mFragment=G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(spectatorZ,spectatorA); 670 pFragment=G4LorentzVector(0,0,0,mFragmen 683 pFragment=G4LorentzVector(0,0,0,mFragment+std::max(0.,theStatisticalExEnergy) ); 671 aProRes.SetMomentum(pFragment); 684 aProRes.SetMomentum(pFragment); 672 685 673 proFrag = theHandler->BreakItUp(aProRes) 686 proFrag = theHandler->BreakItUp(aProRes); 674 687 675 boost_fragments = G4LorentzRotation(pSpe 688 boost_fragments = G4LorentzRotation(pSpectators.boostVector()); 676 689 677 // G4cout << " Fragment a,z, Mass Fr 690 // G4cout << " Fragment a,z, Mass Fragment, mass spect-mom, exitationE " 678 // << spectatorA <<" "<< spectator 691 // << spectatorA <<" "<< spectatorZ <<" "<< mFragment <<" " 679 // << momentum.mag() <<" "<< momen 692 // << momentum.mag() <<" "<< momentum.mag() - mFragment 680 // << " "<<theStatisticalExEnergy 693 // << " "<<theStatisticalExEnergy 681 // << " "<< boost_fragments*pFragm 694 // << " "<< boost_fragments*pFragment<< G4endl; 682 G4ReactionProductVector::iterator ispect 695 G4ReactionProductVector::iterator ispectator; 683 for (ispectator=spectators->begin();ispe 696 for (ispectator=spectators->begin();ispectator!=spectators->end();ispectator++) 684 { 697 { 685 delete *ispectator; 698 delete *ispectator; 686 } 699 } 687 } 700 } 688 else if(spectatorA!=0) 701 else if(spectatorA!=0) 689 { 702 { 690 G4ReactionProductVector::iterator ispecta 703 G4ReactionProductVector::iterator ispectator; 691 for (ispectator=spectators->begin();ispec 704 for (ispectator=spectators->begin();ispectator!=spectators->end();ispectator++) 692 { 705 { 693 (*ispectator)->SetNewlyAdded(true); 706 (*ispectator)->SetNewlyAdded(true); 694 cascaders->push_back(*ispectator); 707 cascaders->push_back(*ispectator); 695 pFinalState+=G4LorentzVector((*ispect 708 pFinalState+=G4LorentzVector((*ispectator)->GetMomentum(),(*ispectator)->GetTotalEnergy()); 696 //G4cout << "BLIC: spectator 709 //G4cout << "BLIC: spectatorA>0, Z=0 from spectator " 697 // << (*ispectator)->GetDefi 710 // << (*ispectator)->GetDefinition()->GetParticleName() << " " 698 // << (*ispectator)->GetMome 711 // << (*ispectator)->GetMomentum()<< " " 699 // << (*ispectator)->GetTota 712 // << (*ispectator)->GetTotalEnergy() << G4endl; 700 } 713 } 701 714 702 } 715 } 703 // / if (spectators) 716 // / if (spectators) 704 delete spectators; 717 delete spectators; 705 718 706 // collect the evaporation part and boost t 719 // collect the evaporation part and boost to spectator frame 707 G4ReactionProductVector::iterator ii; 720 G4ReactionProductVector::iterator ii; 708 if(proFrag) 721 if(proFrag) 709 { 722 { 710 for(ii=proFrag->begin(); ii!=proFrag->en 723 for(ii=proFrag->begin(); ii!=proFrag->end(); ii++) 711 { 724 { 712 (*ii)->SetNewlyAdded(true); 725 (*ii)->SetNewlyAdded(true); 713 G4LorentzVector tmp((*ii)->GetMomentu 726 G4LorentzVector tmp((*ii)->GetMomentum(),(*ii)->GetTotalEnergy()); 714 tmp *= boost_fragments; 727 tmp *= boost_fragments; 715 (*ii)->SetMomentum(tmp.vect()); 728 (*ii)->SetMomentum(tmp.vect()); 716 (*ii)->SetTotalEnergy(tmp.e()); 729 (*ii)->SetTotalEnergy(tmp.e()); 717 // result->push_back(*ii); 730 // result->push_back(*ii); 718 pFragments += tmp; 731 pFragments += tmp; 719 } 732 } 720 } 733 } 721 734 722 // G4cout << "Fragmented p, momentum, de 735 // G4cout << "Fragmented p, momentum, delta " << pFragments <<" "<<momentum 723 // <<" "<< pFragments-momentum < 736 // <<" "<< pFragments-momentum << G4endl; 724 737 725 // correct p/E of Cascade secondaries 738 // correct p/E of Cascade secondaries 726 G4LorentzVector pCas=pInitialState - pFragm 739 G4LorentzVector pCas=pInitialState - pFragments; 727 740 728 //G4cout <<"BLIC: Going to correct from 741 //G4cout <<"BLIC: Going to correct from " << pFinalState << " to " << pCas << G4endl; 729 // the creation of excited fragment did vi 742 // the creation of excited fragment did violate E/p, so correct cascaders to get overall conservation. 730 G4bool EnergyIsCorrect=EnergyAndMomentumCor 743 G4bool EnergyIsCorrect=EnergyAndMomentumCorrector(cascaders, pCas); 731 if ( ! EnergyIsCorrect && debug_G4BinaryLig 744 if ( ! EnergyIsCorrect && debug_G4BinaryLightIonReactionResults) 732 { 745 { 733 G4cout << "G4BinaryLightIonReaction E/P 746 G4cout << "G4BinaryLightIonReaction E/P correction for nucleus failed, will try to correct overall" << G4endl; 734 } 747 } 735 748 736 // Add deexcitation secondaries 749 // Add deexcitation secondaries 737 if(proFrag) 750 if(proFrag) 738 { 751 { 739 for(ii=proFrag->begin(); ii!=proFrag->en 752 for(ii=proFrag->begin(); ii!=proFrag->end(); ii++) 740 { 753 { 741 cascaders->push_back(*ii); 754 cascaders->push_back(*ii); 742 } 755 } 743 delete proFrag; 756 delete proFrag; 744 } 757 } 745 //G4cout << "EnergyIsCorrect? " << Energ 758 //G4cout << "EnergyIsCorrect? " << EnergyIsCorrect << G4endl; 746 if ( ! EnergyIsCorrect ) 759 if ( ! EnergyIsCorrect ) 747 { 760 { 748 // G4cout <<" ! EnergyIsCorrect " << 761 // G4cout <<" ! EnergyIsCorrect " << pFinalState << " to " << pInitialState << G4endl; 749 if (! EnergyAndMomentumCorrector(cascade 762 if (! EnergyAndMomentumCorrector(cascaders,pInitialState)) 750 { 763 { 751 if(debug_G4BinaryLightIonReactionResu 764 if(debug_G4BinaryLightIonReactionResults) 752 G4cout << "G4BinaryLightIonReactio 765 G4cout << "G4BinaryLightIonReaction E/P corrections failed" << G4endl; 753 } 766 } 754 } 767 } 755 768 756 } 769 } 757 770 758 771