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