Ignore:
Timestamp:
Nov 5, 2010, 3:45:55 PM (14 years ago)
Author:
garnier
Message:

update ti head

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/source/processes/hadronic/models/cascade/cascade/src/G4BigBanger.cc

    r1337 r1340  
    2424// ********************************************************************
    2525//
    26 // $Id: G4BigBanger.cc,v 1.30 2010/06/25 09:43:58 gunter Exp $
    27 // Geant4 tag: $Name: geant4-09-04-beta-01 $
     26// $Id: G4BigBanger.cc,v 1.40 2010/09/28 20:15:00 mkelsey Exp $
     27// Geant4 tag: $Name: hadr-casc-V09-03-85 $
    2828//
    2929// 20100114  M. Kelsey -- Remove G4CascadeMomentum, use G4LorentzVector directly
    … …  
    3434// 20100413  M. Kelsey -- Pass G4CollisionOutput by ref to ::collide()
    3535// 20100517  M. Kelsey -- Inherit from common base class, clean up code
     36// 20100628  M. Kelsey -- Use old "bindingEnergy" fn as wrapper, add balance
     37//              checking after bang.
     38// 20100630  M. Kelsey -- Just do simple boost for target, instead of using
     39//              G4LorentzConverter with dummy bullet.
     40// 20100701  M. Kelsey -- Re-throw momentum list, not just angles!
     41// 20100714  M. Kelsey -- Move conservation checking to base class
     42// 20100726  M. Kelsey -- Move std::vector<> buffer to .hh file
     43// 20100923  M. Kelsey -- Migrate to integer A and Z
    3644
    3745#include "G4BigBanger.hh"
    … …  
    4149#include "G4InuclSpecialFunctions.hh"
    4250#include "G4ParticleLargerEkin.hh"
    43 #include "G4LorentzConvertor.hh"
    44 #include "G4HadTmpUtil.hh"
    45 #include "G4NucleiProperties.hh"
    4651#include <algorithm>
    4752
    … …  
    5055typedef std::vector<G4InuclElementaryParticle>::iterator particleIterator;
    5156
    52 G4BigBanger::G4BigBanger() : G4VCascadeCollider("G4BigBanger") {}
     57G4BigBanger::G4BigBanger() : G4CascadeColliderBase("G4BigBanger") {}
    5358
    5459void
    … …  
    5661                     G4CollisionOutput& output) {
    5762
    58   if (verboseLevel > 3) {
    59     G4cout << " >>> G4BigBanger::collide" << G4endl;
    60   }
     63  if (verboseLevel) G4cout << " >>> G4BigBanger::collide" << G4endl;
    6164
    6265  // primitive explosion model A -> nucleons to prevent too exotic evaporation
    63 
    64   const G4double small_ekin = 1.0e-6;
    65 
    66   G4LorentzVector totscm;
    67   G4LorentzVector totlab;
    6866
    6967  G4InuclNuclei* nuclei_target = dynamic_cast<G4InuclNuclei*>(target);
    … …  
    7371  }
    7472
    75   G4double A = nuclei_target->getA();
    76   G4double Z = nuclei_target->getZ();
     73  G4int A = nuclei_target->getA();
     74  G4int Z = nuclei_target->getZ();
     75
    7776  G4LorentzVector PEX = nuclei_target->getMomentum();
    7877  G4double EEXS = nuclei_target->getExitationEnergy();
    79   G4InuclElementaryParticle dummy(small_ekin, 1);
    80   G4LorentzConvertor toTheNucleiSystemRestFrame;
    81  
    82   toTheNucleiSystemRestFrame.setBullet(dummy);
    83   toTheNucleiSystemRestFrame.setTarget(nuclei_target);
    84   toTheNucleiSystemRestFrame.toTheTargetRestFrame();
    85  
    86   G4double etot = (EEXS - G4NucleiProperties::GetBindingEnergy(G4lrint(A), G4lrint(Z) ) ) * MeV/GeV;  // To Bertini units
     78
     79  G4ThreeVector toTheLabFrame = PEX.boostVector();      // From rest to lab
     80
     81  // This "should" be difference between E-target and sum of m(nucleons)
     82  G4double etot = (EEXS - bindingEnergy(A,Z)) * MeV/GeV;  // To Bertini units
    8783  if (etot < 0.0) etot = 0.0;
    8884 
    … …  
    9086    G4cout << " BigBanger: target " << G4endl;
    9187    nuclei_target->printParticle();
    92     G4cout << " BigBanger: a " << A << " z " << Z << " eexs " << EEXS << " etot " <<
    93       etot << " nm " << nuclei_target->getMass() << G4endl;
    94   }
    95  
     88    G4cout << " etot " << etot << G4endl;
     89  }
     90
     91  if (verboseLevel > 3) {
     92    G4LorentzVector PEXrest = PEX;
     93    PEXrest.boost(-toTheLabFrame);
     94    G4cout << " target rest frame: px " << PEXrest.px() << " py "
     95           << PEXrest.py() << " pz " << PEXrest.pz() << " E " << PEXrest.e()
     96           << G4endl;
     97  }
     98
    9699  generateBangInSCM(etot, A, Z);
    97100 
    … …  
    102105  }
    103106
    104   if(!particles.empty()) { // convert back to Lab
    105     particleIterator ipart;
    106    
    107     for(ipart = particles.begin(); ipart != particles.end(); ipart++) {
    108       if (verboseLevel > 2) {
    109         totscm += ipart->getMomentum();
    110       }
    111       G4LorentzVector mom =
    112         toTheNucleiSystemRestFrame.backToTheLab(ipart->getMomentum());
    113       ipart->setMomentum(mom);
    114      
    115       if (verboseLevel > 2) {
    116         mom = ipart->getMomentum();
    117         totlab += mom;
    118       }
    119     }
    120 
    121     std::sort(particles.begin(), particles.end(), G4ParticleLargerEkin());
    122 
    123     if (verboseLevel > 2) {
    124       G4cout << " In SCM: total outgoing momentum " << G4endl
    125              << " E " << totscm.e() << " px " << totscm.x()
    126              << " py " << totscm.y() << " pz " << totscm.z() << G4endl;
    127       G4cout << " In Lab: mom cons " << G4endl
    128              << " E " << PEX.e() + 0.001 * EEXS - totlab.e()
    129              << " px " << PEX.x() - totlab.x()
    130              << " py " << PEX.y() - totlab.y()
    131              << " pz " << PEX.z() - totlab.z() << G4endl;
    132     }
     107  if (particles.empty()) {      // No bang!  Don't know why...
     108    G4cerr << " >>> G4BigBanger unable to process fragment "
     109           << nuclei_target->getDefinition()->GetParticleName() << G4endl;
     110
     111    // FIXME:  This will violate baryon number, momentum, energy, etc.
     112    return;
     113  }
     114
     115  // convert back to Lab
     116  G4LorentzVector totscm;
     117  G4LorentzVector totlab;
     118
     119  if (verboseLevel > 2) G4cout << " BigBanger: boosting to lab" << G4endl;
     120
     121  particleIterator ipart;
     122  for(ipart = particles.begin(); ipart != particles.end(); ipart++) {
     123    G4LorentzVector mom = ipart->getMomentum();
     124    if (verboseLevel > 2) totscm += mom;
     125
     126    mom.boost(toTheLabFrame);
     127    if (verboseLevel > 2) totlab += mom;
     128
     129    ipart->setMomentum(mom);
     130    if (verboseLevel > 2) ipart->printParticle();
     131  }
     132 
     133  std::sort(particles.begin(), particles.end(), G4ParticleLargerEkin());
     134
     135  validateOutput(0, target, particles);         // Checks <vector> directly
     136 
     137  if (verboseLevel > 2) {
     138    G4cout << " In SCM: total outgoing momentum " << G4endl
     139           << " E " << totscm.e() << " px " << totscm.x()
     140           << " py " << totscm.y() << " pz " << totscm.z() << G4endl;
     141    G4cout << " In Lab: mom cons " << G4endl
     142           << " E " << PEX.e() - totlab.e()     // PEX now includes EEXS
     143           << " px " << PEX.x() - totlab.x()
     144           << " py " << PEX.y() - totlab.y()
     145           << " pz " << PEX.z() - totlab.z() << G4endl;
    133146  }
    134147
    135148  output.addOutgoingParticles(particles);
    136 
    137   return;
    138149}                   
    139150
    140 void G4BigBanger::generateBangInSCM(G4double etot, G4double a, G4double z) {
     151void G4BigBanger::generateBangInSCM(G4double etot, G4int a, G4int z) {
    141152  if (verboseLevel > 3) {
    142153    G4cout << " >>> G4BigBanger::generateBangInSCM" << G4endl;
    143154  }
    144155
    145   // Proton and neutron masses
    146   const G4double mp = G4InuclElementaryParticle::getParticleMass(1);
    147   const G4double mn = G4InuclElementaryParticle::getParticleMass(2);
    148 
    149156  const G4double ang_cut = 0.9999;
    150157  const G4int itry_max = 1000;
    151158 
    152   G4int ia = int(a + 0.1);
    153   G4int iz = int(z + 0.1);
    154 
    155   if (verboseLevel > 2) {
    156     G4cout << " ia " << ia << " iz " << iz << G4endl;
     159  if (verboseLevel > 2) {
     160    G4cout << " a " << a << " z " << z << G4endl;
    157161  }
    158162
    159163  particles.clear();    // Reset output vector before filling
    160164 
    161   if(ia == 1) {
    162     // abnormal situation
    163     G4double m = iz > 0 ? mp : mn;
    164     G4double pmod = std::sqrt((etot + 2.0 * m) * etot);
    165     G4LorentzVector mom = generateWithRandomAngles(pmod, m);
    166 
    167     G4int knd = iz > 0 ? 1 : 2;
    168 
    169     particles.push_back(G4InuclElementaryParticle(mom, knd, 8)); // modelId included
    170 
     165  if (a == 1) {         // Special -- bare nucleon doesn't really "explode"
     166    G4int knd = (z>0) ? 1 : 2;
     167    particles.push_back(G4InuclElementaryParticle(knd)); // zero momentum
    171168    return;
    172   }; 
     169  }
    173170     
    174   generateMomentumModules(etot, a, z);
     171  // NOTE:  If distribution fails, need to regenerate magnitudes and angles!
     172  //*** generateMomentumModules(etot, a, z);
     173
     174  scm_momentums.reserve(a);
     175  G4LorentzVector tot_mom;
     176
    175177  G4bool bad = true;
    176178  G4int itry = 0;
    177 
    178179  while(bad && itry < itry_max) {
    179180    itry++;
    180     std::vector<G4LorentzVector> scm_momentums;
    181     G4LorentzVector tot_mom;
    182 
    183     if(ia == 2) {
    184       // FIXME:  This isn't actually a correct four-vector, wrong mass!
     181    scm_momentums.clear();
     182
     183    generateMomentumModules(etot, a, z);
     184    if (a == 2) {
     185      // This is only a three-vector, not a four-vector
    185186      G4LorentzVector mom = generateWithRandomAngles(momModules[0]);
    186 
    187       tot_mom += mom;           
    188 
    189187      scm_momentums.push_back(mom);
    190 
    191       G4LorentzVector mom1 = -mom;
    192 
    193       scm_momentums.push_back(mom1); 
     188      scm_momentums.push_back(-mom);    // Only safe since three-vector!
    194189      bad = false;
    195190    } else {
    196       for(G4int i = 0; i < ia - 2; i++) {
    197         // FIXME:  This isn't actually a correct four-vector, wrong mass!
     191      tot_mom *= 0.;            // Easy way to reset accumulator
     192
     193      for(G4int i = 0; i < a-2; i++) {          // All but last two are thrown
     194      // This is only a three-vector, not a four-vector
    198195        G4LorentzVector mom = generateWithRandomAngles(momModules[i]);
    199 
     196        scm_momentums.push_back(mom);
    200197        tot_mom += mom;         
    201 
    202         scm_momentums.push_back(mom);
    203198      };
    204199
    205200      //                handle last two
    206201      G4double tot_mod = tot_mom.rho();
    207       G4double ct = -0.5 * (tot_mod * tot_mod + momModules[ia - 2] * momModules[ia - 2] -
    208                             momModules[ia - 1] * momModules[ia - 1]) / tot_mod / momModules[ia - 2];
    209 
    210       if (verboseLevel > 2) {
    211         G4cout << " ct last " << ct << G4endl;
    212       }
     202      G4double ct = -0.5*(tot_mod*tot_mod + momModules[a-2]*momModules[a-2]
     203                          - momModules[a-1]*momModules[a-1]) / tot_mod
     204        / momModules[a-2];
     205
     206      if (verboseLevel > 2) G4cout << " ct last " << ct << G4endl;
    213207 
    214208      if(std::fabs(ct) < ang_cut) {
    215         // FIXME:  These are not actually four-vectors, just three-momenta
    216         G4LorentzVector mom2 = generateWithFixedTheta(ct, momModules[ia - 2]);
    217         //       rotate to the normal system
     209        // This is only a three-vector, not a four-vector
     210        G4LorentzVector mom2 = generateWithFixedTheta(ct, momModules[a - 2]);
     211
     212        // rotate to the normal system
    218213        G4LorentzVector apr = tot_mom/tot_mod;
    219214        G4double a_tr = std::sqrt(apr.x()*apr.x() + apr.y()*apr.y());
    … …  
    224219
    225220        scm_momentums.push_back(mom);
    226         //               and the last one
    227         G4LorentzVector mom1= - mom - tot_mom;
     221
     222        // and the last one (again, not actually a four-vector!)
     223        G4LorentzVector mom1 = -mom - tot_mom;
     224
    228225        scm_momentums.push_back(mom1); 
    229226        bad = false;
    230       };
    231     };   
    232     if(!bad) {
    233       for(G4int i = 0; i < ia; i++) {
    234         G4int knd = i < iz ? 1 : 2;
    235 
    236         particles.push_back(G4InuclElementaryParticle(scm_momentums[i], knd, 8));
    237       };
     227      } // if (abs(ct) < ang_cut)
     228    }   // (a > 2)
     229  }     // while (bad && itry<itry_max)
     230
     231  if (!bad) {
     232    for(G4int i = 0; i < a; i++) {
     233      G4int knd = i < z ? 1 : 2;
     234      particles.push_back(G4InuclElementaryParticle(scm_momentums[i], knd, 8));
    238235    };
    239   }; 
    240   if (verboseLevel > 2) {
    241     if(itry == itry_max) G4cout << " BigBanger -> can not generate bang " << G4endl;
     236  };
     237
     238  if (verboseLevel > 2) {
     239    if (itry == itry_max) G4cout << " BigBanger -> can not generate bang " << G4endl;
    242240  }
    243241
    … …  
    245243}
    246244           
    247 void G4BigBanger::generateMomentumModules(G4double etot, G4double a, G4double z) {
     245void G4BigBanger::generateMomentumModules(G4double etot, G4int a, G4int z) {
    248246  if (verboseLevel > 3) {
    249247    G4cout << " >>> G4BigBanger::generateMomentumModules" << G4endl;
    … …  
    256254  momModules.clear();           // Reset buffer for filling
    257255
    258   G4int ia = G4int(a + 0.1);
    259   G4int iz = G4int(z + 0.1);
    260 
    261256  G4double xtot = 0.0;
    262   G4double promax = maxProbability(a);
    263  
    264   G4int i;
    265   for(i = 0; i < ia; i++) {
    266     G4double x = generateX(ia, a, promax);
    267 
    268     if (verboseLevel > 2) {
    269       G4cout << " i " << i << " x " << x << G4endl;
     257
     258  if (a > 2) {                  // For "large" nuclei, energy is distributed
     259    G4double promax = maxProbability(a);
     260   
     261    for(G4int i = 0; i < a; i++) {
     262      G4double x = generateX(a, promax);
     263     
     264      if (verboseLevel > 2) {
     265        G4cout << " i " << i << " x " << x << G4endl;
     266      }
     267      momModules.push_back(x);
     268      xtot += x;
    270269    }
    271     momModules.push_back(x);
    272     xtot += x;
    273   };
    274   for(i = 0; i < ia; i++) {
    275     G4double m = i < iz ? mp : mn;
     270  } else {                      // Two-body case is special, must be 50%
     271    xtot = 1.;
     272    momModules.push_back(0.5);
     273    momModules.push_back(0.5);
     274  }
     275
     276  for(G4int i = 0; i < a; i++) {
     277    G4double m = i < z ? mp : mn;
    276278
    277279    momModules[i] = momModules[i] * etot / xtot;
    … …  
    286288}
    287289
    288 G4double G4BigBanger::xProbability(G4double x, G4int ia) const {
    289 
    290 
    291   if (verboseLevel > 3) {
    292     G4cout << " >>> G4BigBanger::xProbability" << G4endl;
    293   }
    294 
    295   G4int ihalf = ia / 2;
     290G4double G4BigBanger::xProbability(G4double x, G4int a) const {
     291  if (verboseLevel > 3) G4cout << " >>> G4BigBanger::xProbability" << G4endl;
     292
    296293  G4double ekpr = 0.0;
    297294
    … …  
    299296    ekpr = x * x;
    300297
    301     if(2 * ihalf == ia) { // even A
    302       ekpr *= std::sqrt(1.0 - x) * std::pow((1.0 - x), G4int(G4double(3 * ia - 6) / 2.0));
     298    if (a%2 == 0) { // even A
     299      ekpr *= std::sqrt(1.0 - x) * std::pow((1.0 - x), (3*a-6)/2);
    303300    }
    304301    else {
    305       ekpr *= std::pow((1.0 - x), G4int(G4double(3 * ia - 5) / 2.0));
     302      ekpr *= std::pow((1.0 - x), (3*a-5)/2);
    306303    };
    307304  };
    … …  
    310307}
    311308
    312 G4double G4BigBanger::maxProbability(G4double a) const {
    313 
     309G4double G4BigBanger::maxProbability(G4int a) const {
    314310  if (verboseLevel > 3) {
    315311    G4cout << " >>> G4BigBanger::maxProbability" << G4endl;
    316312  }
    317313
    318   return xProbability(1.0 / (a - 1.0) / 1.5, G4int(a + 0.1));
    319 }
    320 
    321 G4double G4BigBanger::generateX(G4int ia,
    322                                 G4double a,
    323                                 G4double promax) const {
    324 
    325   if (verboseLevel > 3) {
    326     G4cout << " >>> G4BigBanger::generateX" << G4endl;
    327   }
     314  return xProbability(2./3./(a-1.0), a);
     315}
     316
     317G4double G4BigBanger::generateX(G4int a, G4double promax) const {
     318  if (verboseLevel > 3) G4cout << " >>> G4BigBanger::generateX" << G4endl;
    328319
    329320  const G4int itry_max = 1000;
    … …  
    335326    x = inuclRndm();
    336327
    337     if(xProbability(x, ia) >= promax * inuclRndm()) return x;
     328    if(xProbability(x, a) >= promax * inuclRndm()) return x;
    338329  };
    339330  if (verboseLevel > 2) {
Note: See TracChangeset for help on using the changeset viewer.