[819] | 1 | // |
---|
| 2 | // ******************************************************************** |
---|
| 3 | // * License and Disclaimer * |
---|
| 4 | // * * |
---|
| 5 | // * The Geant4 software is copyright of the Copyright Holders of * |
---|
| 6 | // * the Geant4 Collaboration. It is provided under the terms and * |
---|
| 7 | // * conditions of the Geant4 Software License, included in the file * |
---|
| 8 | // * LICENSE and available at http://cern.ch/geant4/license . These * |
---|
| 9 | // * include a list of copyright holders. * |
---|
| 10 | // * * |
---|
| 11 | // * Neither the authors of this software system, nor their employing * |
---|
| 12 | // * institutes,nor the agencies providing financial support for this * |
---|
| 13 | // * work make any representation or warranty, express or implied, * |
---|
| 14 | // * regarding this software system or assume any liability for its * |
---|
| 15 | // * use. Please see the license in the file LICENSE and URL above * |
---|
| 16 | // * for the full disclaimer and the limitation of liability. * |
---|
| 17 | // * * |
---|
| 18 | // * This code implementation is the result of the scientific and * |
---|
| 19 | // * technical work of the GEANT4 collaboration. * |
---|
| 20 | // * By using, copying, modifying or distributing the software (or * |
---|
| 21 | // * any work based on the software) you agree to acknowledge its * |
---|
| 22 | // * use in resulting scientific publications, and indicate your * |
---|
| 23 | // * acceptance of all terms of the Geant4 Software license. * |
---|
| 24 | // ******************************************************************** |
---|
| 25 | // |
---|
| 26 | #include "G4BinaryLightIonReaction.hh" |
---|
| 27 | #include "G4LorentzVector.hh" |
---|
| 28 | #include "G4LorentzRotation.hh" |
---|
| 29 | #include <algorithm> |
---|
| 30 | #include "G4ReactionProductVector.hh" |
---|
| 31 | #include <vector> |
---|
| 32 | #include "G4ping.hh" |
---|
| 33 | #include "G4Delete.hh" |
---|
| 34 | #include "G4Neutron.hh" |
---|
| 35 | #include "G4VNuclearDensity.hh" |
---|
| 36 | #include "G4FermiMomentum.hh" |
---|
| 37 | #include "G4HadTmpUtil.hh" |
---|
| 38 | #include <cmath> |
---|
| 39 | |
---|
[1196] | 40 | G4BinaryLightIonReaction::G4BinaryLightIonReaction() |
---|
| 41 | : G4HadronicInteraction("Binary Cascade"), theModel() , |
---|
| 42 | theHandler(0) , theProjectileFragmentation(0) |
---|
| 43 | { |
---|
| 44 | theHandler= new G4ExcitationHandler; |
---|
| 45 | SetPrecompound(new G4PreCompoundModel(theHandler)); |
---|
| 46 | } |
---|
[819] | 47 | |
---|
[1196] | 48 | G4HadFinalState *G4BinaryLightIonReaction:: |
---|
[819] | 49 | ApplyYourself(const G4HadProjectile &aTrack, G4Nucleus & targetNucleus ) |
---|
[1196] | 50 | { |
---|
[819] | 51 | static G4int eventcounter=0; |
---|
| 52 | eventcounter++; |
---|
| 53 | if(getenv("BLICDEBUG") ) G4cerr << " ######### Binary Light Ion Reaction number starts ######### "<<eventcounter<<G4endl; |
---|
| 54 | G4ping debug("debug_G4BinaryLightIonReaction"); |
---|
[1340] | 55 | G4int a1=aTrack.GetDefinition()->GetBaryonNumber(); |
---|
| 56 | G4int z1=G4lrint(aTrack.GetDefinition()->GetPDGCharge()); |
---|
| 57 | G4int a2=targetNucleus.GetA_asInt(); |
---|
| 58 | G4int z2=targetNucleus.GetZ_asInt(); |
---|
[819] | 59 | debug.push_back(a1); |
---|
| 60 | debug.push_back(z1); |
---|
| 61 | debug.push_back(a2); |
---|
| 62 | debug.push_back(z2); |
---|
| 63 | // debug.push_back(m2); |
---|
| 64 | G4LorentzVector mom(aTrack.Get4Momentum()); |
---|
| 65 | debug.push_back(mom); |
---|
| 66 | debug.dump(); |
---|
| 67 | G4LorentzRotation toBreit(mom.boostVector()); |
---|
| 68 | |
---|
| 69 | G4bool swapped = false; |
---|
| 70 | if(a2<a1) |
---|
| 71 | { |
---|
| 72 | swapped = true; |
---|
[1340] | 73 | G4int tmp(0); |
---|
[819] | 74 | tmp = a2; a2=a1; a1=tmp; |
---|
| 75 | tmp = z2; z2=z1; z1=tmp; |
---|
[1340] | 76 | G4double m1=G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(z1,a1); |
---|
[819] | 77 | G4LorentzVector it(m1, G4ThreeVector(0,0,0)); |
---|
| 78 | mom = toBreit*it; |
---|
| 79 | } |
---|
| 80 | |
---|
| 81 | G4ReactionProductVector * result = NULL; |
---|
| 82 | G4ReactionProductVector * cascaders= new G4ReactionProductVector; |
---|
| 83 | G4double m_nucl(0); // to check energy balance |
---|
| 84 | |
---|
| 85 | |
---|
[1340] | 86 | // G4double m1=G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(z1,a1); |
---|
[819] | 87 | // G4cout << "Entering the decision point " |
---|
| 88 | // << (mom.t()-mom.mag())/a1 << " " |
---|
| 89 | // << a1<<" "<< z1<<" " |
---|
| 90 | // << a2<<" "<< z2<<G4endl |
---|
| 91 | // << " "<<mom.t()-mom.mag()<<" " |
---|
| 92 | // << mom.t()- m1<<G4endl; |
---|
| 93 | if( (mom.t()-mom.mag())/a1 < 50*MeV ) |
---|
| 94 | { |
---|
| 95 | // G4cout << "Using pre-compound only, E= "<<mom.t()-mom.mag()<<G4endl; |
---|
| 96 | // m_nucl = mom.mag(); |
---|
| 97 | delete cascaders; |
---|
| 98 | G4Fragment aPreFrag; |
---|
| 99 | aPreFrag.SetA(a1+a2); |
---|
| 100 | aPreFrag.SetZ(z1+z2); |
---|
[1340] | 101 | aPreFrag.SetNumberOfParticles(a1); |
---|
| 102 | aPreFrag.SetNumberOfCharged(z1); |
---|
[819] | 103 | aPreFrag.SetNumberOfHoles(0); |
---|
| 104 | G4ThreeVector plop(0.,0., mom.vect().mag()); |
---|
[1340] | 105 | G4double m2=G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(z2,a2); |
---|
[819] | 106 | m_nucl=m2; |
---|
| 107 | G4LorentzVector aL(mom.t()+m2, plop); |
---|
| 108 | aPreFrag.SetMomentum(aL); |
---|
| 109 | G4ParticleDefinition * preFragDef; |
---|
| 110 | preFragDef = G4ParticleTable::GetParticleTable() |
---|
[1340] | 111 | ->FindIon(z1+z2,a1+a2,0,z1+z2); |
---|
[819] | 112 | aPreFrag.SetParticleDefinition(preFragDef); |
---|
| 113 | |
---|
| 114 | // G4cout << "Fragment INFO "<< a1+a2 <<" "<<z1+z2<<" " |
---|
| 115 | // << aL <<" "<<preFragDef->GetParticleName()<<G4endl; |
---|
[1196] | 116 | cascaders = theProjectileFragmentation->DeExcite(aPreFrag); |
---|
[819] | 117 | G4double tSum = 0; |
---|
| 118 | for(size_t count = 0; count<cascaders->size(); count++) |
---|
| 119 | { |
---|
| 120 | cascaders->operator[](count)->SetNewlyAdded(true); |
---|
| 121 | tSum += cascaders->operator[](count)->GetKineticEnergy(); |
---|
| 122 | } |
---|
| 123 | // G4cout << "Exiting pre-compound only, E= "<<tSum<<G4endl; |
---|
| 124 | } |
---|
| 125 | else |
---|
| 126 | { |
---|
| 127 | |
---|
| 128 | |
---|
| 129 | G4V3DNucleus * fancyNucleus = NULL; |
---|
| 130 | G4Fancy3DNucleus * projectile = NULL; |
---|
| 131 | G4double m1(0) ,m2(0); |
---|
| 132 | G4LorentzVector it; |
---|
| 133 | |
---|
| 134 | G4FermiMomentum theFermi; |
---|
| 135 | G4int tryCount(0); |
---|
| 136 | while(!result) |
---|
| 137 | { |
---|
| 138 | projectile = new G4Fancy3DNucleus; |
---|
| 139 | projectile->Init(a1, z1); |
---|
| 140 | projectile->CenterNucleons(); |
---|
| 141 | m1=G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass( |
---|
| 142 | projectile->GetCharge(),projectile->GetMassNumber()); |
---|
| 143 | it=toBreit * G4LorentzVector(m1,G4ThreeVector(0,0,0)); |
---|
| 144 | fancyNucleus = new G4Fancy3DNucleus; |
---|
| 145 | fancyNucleus->Init(a2, z2); |
---|
| 146 | m2=G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass( |
---|
| 147 | fancyNucleus->GetCharge(),fancyNucleus->GetMassNumber()); |
---|
| 148 | m_nucl = ( swapped ) ? m1 : m2; |
---|
| 149 | // G4cout << " mass table, nucleus, delta : " << m2 <<" "<< fancyNucleus->GetMass() |
---|
| 150 | // <<" "<<m2-fancyNucleus->GetMass() << G4endl; |
---|
| 151 | G4double impactMax = fancyNucleus->GetOuterRadius()+projectile->GetOuterRadius(); |
---|
| 152 | // G4cout << "out radius - nucleus - projectile " << fancyNucleus->GetOuterRadius()/fermi << " - " << projectile->GetOuterRadius()/fermi << G4endl; |
---|
| 153 | G4double aX=(2.*G4UniformRand()-1.)*impactMax; |
---|
| 154 | G4double aY=(2.*G4UniformRand()-1.)*impactMax; |
---|
| 155 | G4ThreeVector pos(aX, aY, -2.*impactMax-5.*fermi); |
---|
| 156 | debug.push_back("Impact parameter"); |
---|
| 157 | debug.push_back(aX); |
---|
| 158 | debug.push_back(aY); |
---|
| 159 | debug.push_back(-2.*impactMax); |
---|
| 160 | debug.dump(); |
---|
| 161 | |
---|
| 162 | G4KineticTrackVector * initalState = new G4KineticTrackVector; |
---|
| 163 | projectile->StartLoop(); |
---|
| 164 | G4Nucleon * aNuc; |
---|
| 165 | G4LorentzVector tmpV(0,0,0,0); |
---|
| 166 | G4LorentzVector nucleonMom(1./a1*mom); |
---|
| 167 | nucleonMom.setZ(nucleonMom.vect().mag()); |
---|
| 168 | nucleonMom.setX(0); |
---|
| 169 | nucleonMom.setY(0); |
---|
| 170 | debug.push_back(" projectile nucleon momentum"); |
---|
| 171 | debug.push_back(nucleonMom); |
---|
| 172 | debug.dump(); |
---|
| 173 | theFermi.Init(a1,z1); |
---|
| 174 | while( (aNuc=projectile->GetNextNucleon()) ) |
---|
| 175 | { |
---|
| 176 | G4LorentzVector p4 = aNuc->GetMomentum(); |
---|
| 177 | tmpV+=p4; |
---|
| 178 | G4ThreeVector nucleonPosition(aNuc->GetPosition()); |
---|
| 179 | G4double density=(projectile->GetNuclearDensity())->GetDensity(nucleonPosition); |
---|
| 180 | nucleonPosition += pos; |
---|
| 181 | G4KineticTrack * it = new G4KineticTrack(aNuc, nucleonPosition, nucleonMom ); |
---|
| 182 | it->SetState(G4KineticTrack::outside); |
---|
| 183 | G4double pfermi= theFermi.GetFermiMomentum(density); |
---|
| 184 | G4double mass = aNuc->GetDefinition()->GetPDGMass(); |
---|
| 185 | G4double Efermi= std::sqrt( sqr(mass) + sqr(pfermi)) - mass; |
---|
| 186 | it->SetProjectilePotential(-Efermi); |
---|
| 187 | initalState->push_back(it); |
---|
| 188 | } |
---|
| 189 | debug.push_back(" Sum of proj. nucleon momentum"); |
---|
| 190 | debug.push_back(tmpV); |
---|
| 191 | debug.dump(); |
---|
| 192 | |
---|
| 193 | result=theModel.Propagate(initalState, fancyNucleus); |
---|
| 194 | debug.push_back("################# Result size"); |
---|
| 195 | if (result) { |
---|
| 196 | debug.push_back(result->size()); |
---|
| 197 | } else { |
---|
| 198 | debug.push_back(" -none-"); |
---|
| 199 | } |
---|
| 200 | debug.dump(); |
---|
| 201 | // std::for_each(initalState->begin(), initalState->end(), Delete<G4KineticTrack>()); |
---|
| 202 | // delete initalState; |
---|
| 203 | |
---|
| 204 | if(! result || result->size()==0) |
---|
| 205 | { |
---|
| 206 | if (result) {delete result; result=0;} |
---|
| 207 | delete fancyNucleus; |
---|
| 208 | delete projectile; |
---|
| 209 | if (++tryCount > 200) |
---|
| 210 | { |
---|
| 211 | // abort!! |
---|
| 212 | |
---|
| 213 | G4cerr << "G4BinaryLightIonReaction no final state for: " << G4endl; |
---|
| 214 | G4cerr << " Primary " << aTrack.GetDefinition() |
---|
| 215 | << ", (A,Z)=(" << aTrack.GetDefinition()->GetBaryonNumber() |
---|
| 216 | << "," << aTrack.GetDefinition()->GetPDGCharge() << ") " |
---|
| 217 | << ", kinetic energy " << aTrack.GetKineticEnergy() |
---|
| 218 | << G4endl; |
---|
[1340] | 219 | G4cerr << " Target nucleus (A,Z)=(" << targetNucleus.GetA_asInt() |
---|
| 220 | << "," << targetNucleus.GetZ_asInt() << G4endl; |
---|
[819] | 221 | G4cerr << " if frequent, please submit above information as bug report" |
---|
| 222 | << G4endl << G4endl; |
---|
| 223 | |
---|
| 224 | theResult.Clear(); |
---|
| 225 | theResult.SetStatusChange(isAlive); |
---|
| 226 | theResult.SetEnergyChange(aTrack.GetKineticEnergy()); |
---|
| 227 | theResult.SetMomentumChange(aTrack.Get4Momentum().vect().unit()); |
---|
| 228 | return &theResult; |
---|
| 229 | |
---|
| 230 | } |
---|
| 231 | } |
---|
| 232 | else |
---|
| 233 | { |
---|
| 234 | break; |
---|
| 235 | } |
---|
| 236 | } |
---|
| 237 | debug.push_back(" Attempts to create final state"); |
---|
| 238 | debug.push_back(tryCount); |
---|
| 239 | debug.dump(); |
---|
| 240 | debug.push_back("################# Through the loop ? "); debug.dump(); |
---|
| 241 | //inverse transformation in case we swapped. |
---|
| 242 | G4int resA(0), resZ(0); |
---|
| 243 | G4Nucleon * aNuc; |
---|
| 244 | // fancyNucleus->StartLoop(); |
---|
| 245 | // while( (aNuc=fancyNucleus->GetNextNucleon()) ) |
---|
| 246 | // { |
---|
| 247 | // G4cout << " tgt Nucleon : " << aNuc->GetDefinition()->GetParticleName() <<" "<< aNuc->AreYouHit() <<" "<<aNuc->GetMomentum()<<G4endl; |
---|
| 248 | // } |
---|
| 249 | G4ReactionProductVector * spectators= new G4ReactionProductVector; |
---|
| 250 | debug.push_back("getting at the hits"); debug.dump(); |
---|
| 251 | // the projectile excitation energy estimate... |
---|
| 252 | G4double theStatisticalExEnergy = 0; |
---|
| 253 | projectile->StartLoop(); |
---|
| 254 | while( (aNuc=projectile->GetNextNucleon()) ) |
---|
| 255 | { |
---|
| 256 | // G4cout << " Nucleon : " << aNuc->GetDefinition()->GetParticleName() <<" "<< aNuc->AreYouHit() <<" "<<aNuc->GetMomentum()<<G4endl; |
---|
| 257 | debug.push_back("getting the hits"); debug.dump(); |
---|
| 258 | if(!aNuc->AreYouHit()) |
---|
| 259 | { |
---|
| 260 | resA++; |
---|
| 261 | resZ+=G4lrint(aNuc->GetDefinition()->GetPDGCharge()); |
---|
| 262 | } |
---|
| 263 | else |
---|
| 264 | { |
---|
| 265 | debug.push_back(" ##### a hit ##### "); debug.dump(); |
---|
| 266 | G4ThreeVector aPosition(aNuc->GetPosition()); |
---|
| 267 | G4double localDensity = projectile->GetNuclearDensity()->GetDensity(aPosition); |
---|
| 268 | G4double localPfermi = theFermi.GetFermiMomentum(localDensity); |
---|
| 269 | G4double nucMass = aNuc->GetDefinition()->GetPDGMass(); |
---|
| 270 | G4double localFermiEnergy = std::sqrt(nucMass*nucMass + localPfermi*localPfermi) - nucMass; |
---|
| 271 | G4double deltaE = localFermiEnergy - (aNuc->GetMomentum().t()-aNuc->GetMomentum().mag()); |
---|
| 272 | theStatisticalExEnergy += deltaE; |
---|
| 273 | } |
---|
| 274 | debug.push_back("collected a hit"); |
---|
| 275 | debug.push_back(aNuc->GetMomentum()); |
---|
| 276 | debug.dump(); |
---|
| 277 | } |
---|
| 278 | delete fancyNucleus; |
---|
| 279 | delete projectile; |
---|
| 280 | G4ping debug("debug_G4BinaryLightIonReaction_1"); |
---|
| 281 | debug.push_back("have the hits. A,Z, excitE"); |
---|
| 282 | debug.push_back(resA); |
---|
| 283 | debug.push_back(resZ); |
---|
| 284 | debug.push_back(theStatisticalExEnergy); |
---|
| 285 | debug.dump(); |
---|
| 286 | // Calculate excitation energy |
---|
| 287 | G4LorentzVector iState = mom; |
---|
| 288 | iState.setT(iState.getT()+m2); |
---|
| 289 | |
---|
| 290 | G4LorentzVector fState(0,0,0,0); |
---|
| 291 | G4LorentzVector pspectators(0,0,0,0); |
---|
| 292 | unsigned int i(0); |
---|
| 293 | // G4int spectA(0),spectZ(0); |
---|
| 294 | for(i=0; i<result->size(); i++) |
---|
| 295 | { |
---|
| 296 | if( (*result)[i]->GetNewlyAdded() ) |
---|
| 297 | { |
---|
| 298 | fState += G4LorentzVector( (*result)[i]->GetMomentum(), (*result)[i]->GetTotalEnergy() ); |
---|
| 299 | cascaders->push_back((*result)[i]); |
---|
| 300 | // G4cout <<" secondary ... "; |
---|
| 301 | debug.push_back("secondary "); |
---|
| 302 | debug.push_back((*result)[i]->GetDefinition()->GetParticleName()); |
---|
| 303 | debug.push_back(G4LorentzVector((*result)[i]->GetMomentum(),(*result)[i]->GetTotalEnergy())); |
---|
| 304 | debug.dump(); |
---|
| 305 | } |
---|
| 306 | else { |
---|
| 307 | // G4cout <<" spectator ... "; |
---|
| 308 | pspectators += G4LorentzVector( (*result)[i]->GetMomentum(), (*result)[i]->GetTotalEnergy() ); |
---|
| 309 | spectators->push_back((*result)[i]); |
---|
| 310 | debug.push_back("spectator "); |
---|
| 311 | debug.push_back((*result)[i]->GetDefinition()->GetParticleName()); |
---|
| 312 | debug.push_back(G4LorentzVector((*result)[i]->GetMomentum(),(*result)[i]->GetTotalEnergy())); |
---|
| 313 | debug.dump(); |
---|
| 314 | // spectA++; |
---|
| 315 | // spectZ+= G4lrint((*result)[i]->GetDefinition()->GetPDGCharge()); |
---|
| 316 | } |
---|
| 317 | |
---|
| 318 | // G4cout << (*result)[i]<< " " |
---|
| 319 | // << (*result)[i]->GetDefinition()->GetParticleName() << " " |
---|
| 320 | // << (*result)[i]->GetMomentum()<< " " |
---|
| 321 | // << (*result)[i]->GetTotalEnergy() << G4endl; |
---|
| 322 | } |
---|
| 323 | // if ( spectA-resA !=0 || spectZ-resZ !=0) |
---|
| 324 | // { |
---|
| 325 | // G4cout << "spect Nucl != spectators: nucl a,z; spect a,z" << |
---|
| 326 | // resA <<" "<< resZ <<" ; " << spectA <<" "<< spectZ << G4endl; |
---|
| 327 | // } |
---|
| 328 | delete result; |
---|
| 329 | debug.push_back(" iState - (fState+pspectators) "); |
---|
| 330 | debug.push_back(iState-fState-pspectators); |
---|
| 331 | debug.dump(); |
---|
| 332 | G4LorentzVector momentum(iState-fState); |
---|
| 333 | G4int loopcount(0); |
---|
| 334 | while (std::abs(momentum.e()-pspectators.e()) > 10*MeV) |
---|
| 335 | { |
---|
| 336 | debug.push_back("the momentum balance"); |
---|
| 337 | debug.push_back(iState); |
---|
| 338 | debug.push_back(fState); |
---|
| 339 | debug.push_back(momentum-pspectators); |
---|
| 340 | debug.push_back(momentum); |
---|
| 341 | debug.dump(); |
---|
| 342 | G4LorentzVector pCorrect(iState-pspectators); |
---|
| 343 | G4bool EnergyIsCorrect=EnergyAndMomentumCorrector(cascaders, pCorrect); |
---|
| 344 | if ( ! EnergyIsCorrect && getenv("debug_G4BinaryLightIonReactionResults")) |
---|
| 345 | { |
---|
| 346 | G4cout << "Warning - G4BinaryLightIonReaction E/P correction for cascaders failed" << G4endl; |
---|
| 347 | } |
---|
| 348 | fState=G4LorentzVector(); |
---|
| 349 | for(i=0; i<cascaders->size(); i++) |
---|
| 350 | { |
---|
| 351 | fState += G4LorentzVector( (*cascaders)[i]->GetMomentum(), (*cascaders)[i]->GetTotalEnergy() ); |
---|
| 352 | } |
---|
| 353 | momentum=iState-fState; |
---|
| 354 | debug.push_back("the momentum balance after correction"); |
---|
| 355 | debug.push_back(iState); |
---|
| 356 | debug.push_back(fState); |
---|
| 357 | debug.push_back(momentum-pspectators); |
---|
| 358 | debug.push_back(momentum); |
---|
| 359 | debug.dump(); |
---|
| 360 | if (++loopcount > 10 ) |
---|
| 361 | { |
---|
| 362 | if ( momentum.vect().mag() > momentum.e() ) |
---|
| 363 | { |
---|
| 364 | G4cerr << "G4BinaryLightIonReaction.cc: Cannot correct 4-momentum of cascade particles" << G4endl; |
---|
| 365 | throw G4HadronicException(__FILE__, __LINE__, "G4BinaryCasacde::ApplyCollision()"); |
---|
| 366 | } else { |
---|
| 367 | break; |
---|
| 368 | } |
---|
| 369 | |
---|
| 370 | } |
---|
| 371 | } |
---|
| 372 | |
---|
| 373 | |
---|
| 374 | |
---|
| 375 | |
---|
| 376 | // call precompound model |
---|
| 377 | G4ReactionProductVector * proFrag = NULL; |
---|
| 378 | G4LorentzVector pFragment; |
---|
| 379 | // G4cout << " == pre boost 1 "<< momentum.e()<< " "<< momentum.mag()<<G4endl; |
---|
| 380 | G4LorentzRotation boost_fragments; |
---|
| 381 | // G4cout << " == post boost 1 "<< momentum.e()<< " "<< momentum.mag()<<G4endl; |
---|
| 382 | // G4LorentzRotation boost_spectator_mom(-momentum.boostVector()); |
---|
| 383 | // G4cout << "- momentum " << boost_spectator_mom * momentum << G4endl; |
---|
| 384 | G4LorentzVector pFragments(0); |
---|
| 385 | if(resZ>0 && resA>1) |
---|
| 386 | { |
---|
| 387 | // Make the fragment |
---|
| 388 | G4Fragment aProRes; |
---|
| 389 | aProRes.SetA(resA); |
---|
| 390 | aProRes.SetZ(resZ); |
---|
| 391 | aProRes.SetNumberOfParticles(0); |
---|
| 392 | aProRes.SetNumberOfCharged(0); |
---|
[1340] | 393 | aProRes.SetNumberOfHoles(a1-resA); |
---|
[819] | 394 | G4double mFragment=G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(resZ,resA); |
---|
| 395 | G4LorentzVector pFragment(0,0,0,mFragment+std::max(0.,theStatisticalExEnergy) ); |
---|
| 396 | aProRes.SetMomentum(pFragment); |
---|
| 397 | G4ParticleDefinition * resDef; |
---|
| 398 | resDef = G4ParticleTable::GetParticleTable()->FindIon(resZ,resA,0,resZ); |
---|
| 399 | aProRes.SetParticleDefinition(resDef); |
---|
[1196] | 400 | proFrag = theHandler->BreakItUp(aProRes); |
---|
[819] | 401 | if ( momentum.vect().mag() > momentum.e() ) |
---|
| 402 | { |
---|
| 403 | G4cerr << "mom check: " << momentum |
---|
| 404 | << " 3.mag "<< momentum.vect().mag() << G4endl |
---|
| 405 | << " .. iState/fState/spectators " << iState <<" " |
---|
| 406 | << fState << " " << pspectators << G4endl |
---|
| 407 | << " .. A,Z " << resA <<" "<< resZ << G4endl; |
---|
| 408 | G4cerr << "G4BinaryLightIonReaction no final state for: " << G4endl; |
---|
| 409 | G4cerr << " Primary " << aTrack.GetDefinition() |
---|
| 410 | << ", (A,Z)=(" << aTrack.GetDefinition()->GetBaryonNumber() |
---|
| 411 | << "," << aTrack.GetDefinition()->GetPDGCharge() << ") " |
---|
| 412 | << ", kinetic energy " << aTrack.GetKineticEnergy() |
---|
| 413 | << G4endl; |
---|
[1340] | 414 | G4cerr << " Target nucleus (A,Z)=(" << targetNucleus.GetA_asInt() |
---|
| 415 | << "," << targetNucleus.GetZ_asInt() << G4endl; |
---|
[819] | 416 | G4cerr << " if frequent, please submit above information as bug report" |
---|
| 417 | << G4endl << G4endl; |
---|
| 418 | |
---|
| 419 | theResult.Clear(); |
---|
| 420 | theResult.SetStatusChange(isAlive); |
---|
| 421 | theResult.SetEnergyChange(aTrack.GetKineticEnergy()); |
---|
| 422 | theResult.SetMomentumChange(aTrack.Get4Momentum().vect().unit()); |
---|
| 423 | return &theResult; |
---|
| 424 | } |
---|
| 425 | |
---|
| 426 | G4LorentzRotation boost_fragments_here(momentum.boostVector()); |
---|
| 427 | boost_fragments = boost_fragments_here; |
---|
| 428 | |
---|
| 429 | // G4cout << " Fragment a,z, Mass Fragment, mass spect-mom, exitationE " |
---|
| 430 | // << resA <<" "<< resZ <<" "<< mFragment <<" " |
---|
| 431 | // << momentum.mag() <<" "<< momentum.mag() - mFragment |
---|
| 432 | // << " "<<theStatisticalExEnergy |
---|
| 433 | // << " "<< boost_fragments*pFragment<< G4endl; |
---|
| 434 | G4ReactionProductVector::iterator ispectator; |
---|
| 435 | for (ispectator=spectators->begin();ispectator!=spectators->end();ispectator++) |
---|
| 436 | { |
---|
| 437 | delete *ispectator; |
---|
| 438 | } |
---|
| 439 | } |
---|
| 440 | else if(resA!=0) |
---|
| 441 | { |
---|
| 442 | G4ReactionProductVector::iterator ispectator; |
---|
| 443 | for (ispectator=spectators->begin();ispectator!=spectators->end();ispectator++) |
---|
| 444 | { |
---|
| 445 | (*ispectator)->SetNewlyAdded(true); |
---|
| 446 | cascaders->push_back(*ispectator); |
---|
| 447 | pFragments+=G4LorentzVector((*ispectator)->GetMomentum(),(*ispectator)->GetTotalEnergy()); |
---|
| 448 | // G4cout << "from spectator " |
---|
| 449 | // << (*ispectator)->GetDefinition()->GetParticleName() << " " |
---|
| 450 | // << (*ispectator)->GetMomentum()<< " " |
---|
| 451 | // << (*ispectator)->GetTotalEnergy() << G4endl; |
---|
| 452 | } |
---|
| 453 | } |
---|
| 454 | if (spectators) delete spectators; |
---|
| 455 | |
---|
| 456 | // collect the evaporation part |
---|
| 457 | debug.push_back("the nucleon count balance"); |
---|
| 458 | debug.push_back(resA); |
---|
| 459 | debug.push_back(resZ); |
---|
| 460 | if(proFrag) debug.push_back(proFrag->size()); |
---|
| 461 | debug.dump(); |
---|
| 462 | G4ReactionProductVector::iterator ii; |
---|
| 463 | if(proFrag) for(ii=proFrag->begin(); ii!=proFrag->end(); ii++) |
---|
| 464 | { |
---|
| 465 | (*ii)->SetNewlyAdded(true); |
---|
| 466 | G4LorentzVector tmp((*ii)->GetMomentum(),(*ii)->GetTotalEnergy()); |
---|
| 467 | tmp *= boost_fragments; |
---|
| 468 | (*ii)->SetMomentum(tmp.vect()); |
---|
| 469 | (*ii)->SetTotalEnergy(tmp.e()); |
---|
| 470 | // result->push_back(*ii); |
---|
| 471 | pFragments += tmp; |
---|
| 472 | } |
---|
| 473 | |
---|
| 474 | // G4cout << "Fragmented p, momentum, delta " << pFragments <<" "<<momentum |
---|
| 475 | // <<" "<< pFragments-momentum << G4endl; |
---|
| 476 | debug.push_back("################# done with evaporation"); debug.dump(); |
---|
| 477 | |
---|
| 478 | // correct p/E of Cascade secondaries |
---|
| 479 | G4LorentzVector pCas=iState - pFragments; |
---|
| 480 | // G4cout <<" Going to correct from " << fState << " to " << pCas << G4endl; |
---|
| 481 | G4bool EnergyIsCorrect=EnergyAndMomentumCorrector(cascaders, pCas); |
---|
| 482 | if ( ! EnergyIsCorrect ) |
---|
| 483 | { |
---|
| 484 | if(getenv("debug_G4BinaryLightIonReactionResults")) |
---|
| 485 | G4cout << "G4BinaryLightIonReaction E/P correction for nucleus failed, will try to correct overall" << G4endl; |
---|
| 486 | } |
---|
| 487 | |
---|
| 488 | // Add deexcitation secondaries |
---|
| 489 | if(proFrag) for(ii=proFrag->begin(); ii!=proFrag->end(); ii++) |
---|
| 490 | { |
---|
| 491 | cascaders->push_back(*ii); |
---|
| 492 | } |
---|
| 493 | if (proFrag) delete proFrag; |
---|
| 494 | |
---|
| 495 | if ( ! EnergyIsCorrect ) |
---|
| 496 | { |
---|
| 497 | if (! EnergyAndMomentumCorrector(cascaders,iState)) |
---|
| 498 | { |
---|
| 499 | if(getenv("debug_G4BinaryLightIonReactionResults")) |
---|
| 500 | G4cout << "G4BinaryLightIonReaction E/P corrections failed" << G4endl; |
---|
| 501 | } |
---|
| 502 | } |
---|
| 503 | |
---|
| 504 | } |
---|
| 505 | |
---|
| 506 | // Rotate to lab |
---|
| 507 | G4LorentzRotation toZ; |
---|
| 508 | toZ.rotateZ(-1*mom.phi()); |
---|
| 509 | toZ.rotateY(-1*mom.theta()); |
---|
| 510 | G4LorentzRotation toLab(toZ.inverse()); |
---|
| 511 | |
---|
| 512 | // Fill the particle change, while rotating. Boost from projectile breit-frame in case we swapped. |
---|
| 513 | // theResult.Clear(); |
---|
| 514 | theResult.Clear(); |
---|
| 515 | theResult.SetStatusChange(stopAndKill); |
---|
| 516 | G4double Etot(0); |
---|
| 517 | size_t i=0; |
---|
| 518 | for(i=0; i<cascaders->size(); i++) |
---|
| 519 | { |
---|
| 520 | if((*cascaders)[i]->GetNewlyAdded()) |
---|
| 521 | { |
---|
| 522 | G4DynamicParticle * aNew = |
---|
| 523 | new G4DynamicParticle((*cascaders)[i]->GetDefinition(), |
---|
| 524 | (*cascaders)[i]->GetTotalEnergy(), |
---|
| 525 | (*cascaders)[i]->GetMomentum() ); |
---|
| 526 | G4LorentzVector tmp = aNew->Get4Momentum(); |
---|
| 527 | if(swapped) |
---|
| 528 | { |
---|
| 529 | tmp*=toBreit.inverse(); |
---|
| 530 | tmp.setVect(-tmp.vect()); |
---|
| 531 | } |
---|
| 532 | tmp *= toLab; |
---|
| 533 | aNew->Set4Momentum(tmp); |
---|
| 534 | theResult.AddSecondary(aNew); |
---|
| 535 | Etot += tmp.e(); |
---|
| 536 | // G4cout << "LIBIC: Secondary " << aNew->GetDefinition()->GetParticleName() |
---|
| 537 | // <<" "<< aNew->GetMomentum() |
---|
| 538 | // <<" "<< aNew->GetTotalEnergy() |
---|
| 539 | // << G4endl; |
---|
| 540 | } |
---|
| 541 | delete (*cascaders)[i]; |
---|
| 542 | } |
---|
| 543 | if(cascaders) delete cascaders; |
---|
| 544 | |
---|
| 545 | G4ping debug1("debug_G4BinaryLightIonReactionResults"); |
---|
| 546 | debug1.push_back("Result analysis, secondaries"); |
---|
| 547 | debug1.push_back(theResult.GetNumberOfSecondaries()); |
---|
| 548 | debug1.dump(); |
---|
| 549 | debug1.push_back(" Energy conservation initial/final/delta(init-final) "); |
---|
| 550 | debug1.push_back(aTrack.GetTotalEnergy() + m_nucl); |
---|
| 551 | debug1.push_back(aTrack.GetTotalEnergy()); |
---|
| 552 | debug1.push_back(m_nucl); |
---|
| 553 | debug1.push_back(Etot); |
---|
| 554 | debug1.push_back(aTrack.GetTotalEnergy() + m_nucl - Etot); |
---|
| 555 | debug1.dump(); |
---|
| 556 | |
---|
| 557 | if(getenv("BLICDEBUG") ) G4cerr << " ######### Binary Light Ion Reaction number ends ######### "<<eventcounter<<G4endl; |
---|
| 558 | |
---|
| 559 | return &theResult; |
---|
| 560 | } |
---|
| 561 | |
---|
| 562 | //**************************************************************************** |
---|
| 563 | G4bool G4BinaryLightIonReaction::EnergyAndMomentumCorrector( |
---|
| 564 | G4ReactionProductVector* Output, G4LorentzVector& TotalCollisionMom) |
---|
| 565 | //**************************************************************************** |
---|
| 566 | { |
---|
| 567 | const int nAttemptScale = 2500; |
---|
| 568 | const double ErrLimit = 1.E-6; |
---|
| 569 | if (Output->empty()) |
---|
| 570 | return TRUE; |
---|
| 571 | G4LorentzVector SumMom(0); |
---|
| 572 | G4double SumMass = 0; |
---|
| 573 | G4double TotalCollisionMass = TotalCollisionMom.m(); |
---|
| 574 | size_t i = 0; |
---|
| 575 | // Calculate sum hadron 4-momenta and summing hadron mass |
---|
| 576 | for(i = 0; i < Output->size(); i++) |
---|
| 577 | { |
---|
| 578 | SumMom += G4LorentzVector((*Output)[i]->GetMomentum(),(*Output)[i]->GetTotalEnergy()); |
---|
| 579 | SumMass += (*Output)[i]->GetDefinition()->GetPDGMass(); |
---|
| 580 | } |
---|
| 581 | // G4cout << " E/P corrector, SumMass, SumMom.m2, TotalMass " |
---|
| 582 | // << SumMass <<" "<< SumMom.m2() <<" "<<TotalCollisionMass<< G4endl; |
---|
| 583 | if (SumMass > TotalCollisionMass) return FALSE; |
---|
| 584 | SumMass = SumMom.m2(); |
---|
| 585 | if (SumMass < 0) return FALSE; |
---|
| 586 | SumMass = std::sqrt(SumMass); |
---|
| 587 | |
---|
| 588 | // Compute c.m.s. hadron velocity and boost KTV to hadron c.m.s. |
---|
| 589 | G4ThreeVector Beta = -SumMom.boostVector(); |
---|
| 590 | // G4cout << " == pre boost 2 "<< SumMom.e()<< " "<< SumMom.mag()<<" "<< Beta <<G4endl; |
---|
| 591 | //--old Output->Boost(Beta); |
---|
| 592 | for(i = 0; i < Output->size(); i++) |
---|
| 593 | { |
---|
| 594 | G4LorentzVector mom = G4LorentzVector((*Output)[i]->GetMomentum(),(*Output)[i]->GetTotalEnergy()); |
---|
| 595 | mom *= Beta; |
---|
| 596 | (*Output)[i]->SetMomentum(mom.vect()); |
---|
| 597 | (*Output)[i]->SetTotalEnergy(mom.e()); |
---|
| 598 | } |
---|
| 599 | |
---|
| 600 | // Scale total c.m.s. hadron energy (hadron system mass). |
---|
| 601 | // It should be equal interaction mass |
---|
| 602 | G4double Scale = 0,OldScale=0; |
---|
| 603 | G4double factor = 1.; |
---|
| 604 | G4int cAttempt = 0; |
---|
| 605 | G4double Sum = 0; |
---|
| 606 | G4bool success = false; |
---|
| 607 | for(cAttempt = 0; cAttempt < nAttemptScale; cAttempt++) |
---|
| 608 | { |
---|
| 609 | Sum = 0; |
---|
| 610 | for(i = 0; i < Output->size(); i++) |
---|
| 611 | { |
---|
| 612 | G4LorentzVector HadronMom = G4LorentzVector((*Output)[i]->GetMomentum(),(*Output)[i]->GetTotalEnergy()); |
---|
| 613 | HadronMom.setVect(HadronMom.vect()+ factor*Scale*HadronMom.vect()); |
---|
| 614 | G4double E = std::sqrt(HadronMom.vect().mag2() + sqr((*Output)[i]->GetDefinition()->GetPDGMass())); |
---|
| 615 | HadronMom.setE(E); |
---|
| 616 | (*Output)[i]->SetMomentum(HadronMom.vect()); |
---|
| 617 | (*Output)[i]->SetTotalEnergy(HadronMom.e()); |
---|
| 618 | Sum += E; |
---|
| 619 | } |
---|
| 620 | OldScale=Scale; |
---|
| 621 | Scale = TotalCollisionMass/Sum - 1; |
---|
| 622 | // G4cout << "E/P corr - " << cAttempt << " " << Scale << G4endl; |
---|
| 623 | if (std::abs(Scale) <= ErrLimit |
---|
| 624 | || OldScale == Scale) // protect 'frozen' situation and divide by 0 in calculating new factor below |
---|
| 625 | { |
---|
| 626 | if (getenv("debug_G4BinaryLightIonReactionResults")) G4cout << "E/p corrector: " << cAttempt << G4endl; |
---|
| 627 | success = true; |
---|
| 628 | break; |
---|
| 629 | } |
---|
| 630 | if ( cAttempt > 10 ) |
---|
| 631 | { |
---|
| 632 | // G4cout << " speed it up? " << std::abs(OldScale/(OldScale-Scale)) << G4endl; |
---|
| 633 | factor=std::max(1.,std::log(std::abs(OldScale/(OldScale-Scale)))); |
---|
| 634 | // G4cout << " ? factor ? " << factor << G4endl; |
---|
| 635 | } |
---|
| 636 | } |
---|
| 637 | |
---|
| 638 | if( (!success) && getenv("debug_G4BinaryLightIonReactionResults")) |
---|
| 639 | { |
---|
| 640 | G4cout << "G4G4BinaryLightIonReaction::EnergyAndMomentumCorrector - Warning"<<G4endl; |
---|
| 641 | G4cout << " Scale not unity at end of iteration loop: "<<TotalCollisionMass<<" "<<Sum<<" "<<Scale<<G4endl; |
---|
| 642 | G4cout << " Increase number of attempts or increase ERRLIMIT"<<G4endl; |
---|
| 643 | } |
---|
| 644 | |
---|
| 645 | // Compute c.m.s. interaction velocity and KTV back boost |
---|
| 646 | Beta = TotalCollisionMom.boostVector(); |
---|
| 647 | //--old Output->Boost(Beta); |
---|
| 648 | for(i = 0; i < Output->size(); i++) |
---|
| 649 | { |
---|
| 650 | G4LorentzVector mom = G4LorentzVector((*Output)[i]->GetMomentum(),(*Output)[i]->GetTotalEnergy()); |
---|
| 651 | mom *= Beta; |
---|
| 652 | (*Output)[i]->SetMomentum(mom.vect()); |
---|
| 653 | (*Output)[i]->SetTotalEnergy(mom.e()); |
---|
| 654 | } |
---|
| 655 | return TRUE; |
---|
| 656 | } |
---|