Ignore:
Timestamp:
Jun 18, 2010, 11:42:07 AM (14 years ago)
Author:
garnier
Message:

update geant4-09-04-beta-cand-01 interfaces-V09-03-09 vis-V09-03-08

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r962 r1315  
    2323// * acceptance of all terms of the Geant4 Software license.          *
    2424// ********************************************************************
     25// $Id: G4LorentzConvertor.cc,v 1.23 2010/05/21 17:56:34 mkelsey Exp $
     26// Geant4 tag: $Name: geant4-09-04-beta-cand-01 $
    2527//
     28// 20100108  Michael Kelsey -- Use G4LorentzVector internally
     29// 20100112  M. Kelsey -- Remove G4CascadeMomentum, use G4LorentzVector directly
     30// 20100308  M. Kelsey -- Bug fix in toTheTargetRestFrame: scm_momentum should
     31//              be data member, not local.
     32// 20100409  M. Kelsey -- Protect std::sqrt(ga) against round-off negatives
     33// 20100519  M. Kelsey -- Add interfaces to pass G4InuclParticles directly
     34
    2635#include "G4LorentzConvertor.hh"
     36#include "G4ThreeVector.hh"
    2737#include "G4HadronicException.hh"
     38#include "G4InuclParticle.hh"
     39
     40
     41const G4double G4LorentzConvertor::small = 1.0e-10;
    2842
    2943G4LorentzConvertor::G4LorentzConvertor()
    30   : verboseLevel(2), degenerated(false) {
     44  : verboseLevel(0), degenerated(false) {
    3145
    3246  if (verboseLevel > 3) {
    … …  
    3549}
    3650
     51// Extract four-vectors from input particles
     52void G4LorentzConvertor::setBullet(const G4InuclParticle* bullet) {
     53  setBullet(bullet->getMomentum());
     54}
     55
     56void G4LorentzConvertor::setTarget(const G4InuclParticle* target) {
     57  setTarget(target->getMomentum());
     58}
     59
     60// Boost bullet and target four-vectors into destired frame
     61
    3762void G4LorentzConvertor::toTheCenterOfMass() {
    38    
    3963  if (verboseLevel > 3) {
    4064    G4cout << " >>> G4LorentzConvertor::toTheCenterOfMass" << G4endl;
    4165  }
    4266
    43   const G4double small = 1.0e-10;
    44 
    45   v2 = 0.0;
    46 
    47   G4double pv = 0.0;
    48 
    49    G4double e_sum = target_mom[0] + bullet_mom[0];
    50 
    51   velocity.resize(4);
    52   G4int i(0);
    53   for(i = 1; i < 4; i++) {
    54     velocity[i] = (target_mom[i] + bullet_mom[i]) / e_sum;
    55     v2 += velocity[i] * velocity[i];
    56     pv += target_mom[i] * velocity[i];
    57   };
    58    
    59   gamma = 1.0 / std::sqrt(std::fabs(1.0 - v2));
    60   ecm_tot = e_sum / gamma;
    61 
    62   G4double pa = 0.0;
    63 
    64   G4double pb = 0.0;
    65 
    66   G4double xx = pv * (gamma - 1.0) / v2 - target_mom[0] * gamma;
    67 
    68   for(i = 1; i < 4; i++) {
    69     scm_momentum[i] = -target_mom[i] - velocity[i] * xx;
    70 
    71     if (verboseLevel > 3) {
    72       G4cout << " i " << i << " pscm(i) " << scm_momentum[i] << G4endl;
    73     }
    74 
    75     pa += scm_momentum[i] * scm_momentum[i];
    76     pb += scm_momentum[i] * velocity[i];
    77   };
    78   ga = v2 - pb * pb / pa;
    79   if(ga < small) {
    80     ga = small;
    81     degenerated = true;
    82 
    83     if (verboseLevel > 3) {
    84       G4cout << " degenerated case " << G4endl;
    85     }
    86 
    87   } else {
    88     ga = std::sqrt(ga);
    89   };
    90 
    91   if (verboseLevel > 3) {
    92     G4cout << " ga " << ga << " v2 " << v2 << " pb " << pb <<
    93       " pb * pb / pa " << pb * pb / pa << " pv " << pv << G4endl;
    94   }
    95 
    96   pscm = std::sqrt(pa);
     67  G4LorentzVector cm4v = target_mom + bullet_mom;
     68  velocity = cm4v.boostVector();
     69
     70  // "SCM" is reverse target momentum in the CM frame
     71  scm_momentum = target_mom;
     72  scm_momentum.boost(-velocity);
     73  scm_momentum.setVect(-scm_momentum.vect());
     74
     75  if (verboseLevel > 3)
     76    G4cout << " i 1 pscm(i) " << scm_momentum.x() << G4endl
     77           << " i 2 pscm(i) " << scm_momentum.y() << G4endl
     78           << " i 3 pscm(i) " << scm_momentum.z() << G4endl;
     79
     80  // Compute kinematic quantities for rotate() functions
     81  v2 = velocity.mag2();
     82  gamma = cm4v.e()/cm4v.m();
     83
     84  ecm_tot = cm4v.m();
     85
     86  G4double pscm = scm_momentum.rho();
     87  G4double pa   = scm_momentum.vect().mag2();
     88  G4double pb   = scm_momentum.vect().dot(velocity);
     89
     90  G4double gasq = v2-pb*pb/pa;
     91  ga = (gasq > small*small) ? std::sqrt(gasq) : 0.;
     92
    9793  gb = pb / pscm;
    9894  gbpp = gb / pscm;
    9995  gapp = ga * pscm;
    100 }
    101 
    102 G4CascadeMomentum G4LorentzConvertor::rotate(const G4CascadeMomentum& mom) const {
    103 
    104   if (verboseLevel > 3) {
    105     G4cout << " >>> G4LorentzConvertor::rotate(G4CascadeMomentum)" << G4endl;
    106   }
    107 
    108   G4CascadeMomentum mom_rot;
    109 
    110   if (verboseLevel > 3) {
    111     G4cout << " ga " << ga << " gbpp " << gbpp << " gapp " << gapp << G4endl; 
    112     G4cout << " gegenerated " << degenerated << G4endl;
    113     G4cout << " before rotation: px " << mom[1] << " py " << mom[2] <<
    114       " pz " << mom[3] << G4endl;
    115   }
    116 
    117   if(degenerated) {
    118     mom_rot = mom;
    119   } else {
    120     mom_rot[1] = mom[1] * (velocity[1] - gbpp * scm_momentum[1]) / ga +
    121       mom[2] * (scm_momentum[2] * velocity[3] - scm_momentum[3] * velocity[2]) / gapp +
    122       mom[3] * scm_momentum[1] / pscm;
    123     mom_rot[2] = mom[1] * (velocity[2] - gbpp * scm_momentum[2]) / ga +
    124       mom[2] * (scm_momentum[3] * velocity[1] - scm_momentum[1] * velocity[3]) / gapp +
    125       mom[3] * scm_momentum[2] / pscm;
    126     mom_rot[3] = mom[1] * (velocity[3] - gbpp * scm_momentum[3]) / ga +
    127       mom[2] * (scm_momentum[1] * velocity[2] - scm_momentum[2] * velocity[1]) / gapp +
    128       mom[3] * scm_momentum[3] / pscm;
     96
     97  degenerated = (ga < small);
     98  if (degenerated && verboseLevel > 3)
     99    G4cout << " degenerated case " << G4endl;
     100
     101  if (verboseLevel > 3) {
     102    G4double pv = target_mom.vect().dot(velocity);
     103    G4cout << " ga " << ga << " v2 " << v2 << " pb " << pb
     104           << " pb * pb / pa " << pb * pb / pa << " pv " << pv << G4endl;
     105  }
     106}
     107
     108void G4LorentzConvertor::toTheTargetRestFrame() {
     109  if (verboseLevel > 3) {
     110    G4cout << " >>> G4LorentzConvertor::toTheTargetRestFrame" << G4endl;
     111  }
     112
     113  velocity = target_mom.boostVector();
     114
     115  // "SCM" is bullet momentum in the target's frame
     116  scm_momentum = bullet_mom;
     117  scm_momentum.boost(-velocity);
     118
     119  if (verboseLevel > 3)
     120    G4cout << " rf: i 1 pscm(i) " << scm_momentum.x() << G4endl
     121           << " rf: i 2 pscm(i) " << scm_momentum.y() << G4endl
     122           << " rf: i 3 pscm(i) " << scm_momentum.z() << G4endl;
     123
     124  // Compute kinematic quantities for rotate() functions
     125  v2 = velocity.mag2();
     126  gamma = target_mom.e() / target_mom.m();
     127
     128  G4double pscm = scm_momentum.rho();
     129  G4double pa   = scm_momentum.vect().mag2();
     130  G4double pb   = velocity.dot(scm_momentum.vect());
     131
     132  G4double gasq = v2-pb*pb/pa;
     133  ga = (gasq > small*small) ? std::sqrt(gasq) : 0.;
     134
     135  gb = pb/pscm;
     136  gbpp = gb/pscm;
     137  gapp = ga*pscm;
     138
     139  degenerated = (ga < small);
     140  if (degenerated && verboseLevel > 3)
     141    G4cout << " degenerated case " << G4endl;
     142}
     143
     144G4LorentzVector
     145G4LorentzConvertor::backToTheLab(const G4LorentzVector& mom) const {
     146  if (verboseLevel > 3) {
     147    G4cout << " >>> G4LorentzConvertor::backToTheLab" << G4endl
     148           << " at rest: px " << mom.x() << " py " << mom.y() << " pz "
     149           << mom.z() << " e " << mom.e() << G4endl
     150           << " v2 " << v2 << G4endl;
     151  }
     152
     153  G4LorentzVector mom1 = mom;
     154  if (v2 > small) mom1.boost(velocity);
     155
     156  if (verboseLevel > 3)
     157    G4cout << " at lab: px " << mom1.x() << " py " << mom1.y() << " pz "
     158           << mom1.z() << G4endl;
     159
     160  return mom1;
     161}
     162
     163
     164// Bullet kinematics in target rest frame (LAB frame, usually)
     165
     166G4double G4LorentzConvertor::getKinEnergyInTheTRS() const {
     167  G4double pv = bullet_mom.vect().dot(target_mom.vect());
     168 
     169   G4double ekin_trf =
     170    (target_mom.e() * bullet_mom.e() - pv) / target_mom.m()
     171    - bullet_mom.m();
     172 
     173  return ekin_trf;
     174}
     175
     176G4double G4LorentzConvertor::getTRSMomentum() const {
     177  G4LorentzVector bmom = bullet_mom;
     178  bmom.boost(-target_mom.boostVector());
     179  return bmom.rho();
     180}
     181
     182G4LorentzVector G4LorentzConvertor::rotate(const G4LorentzVector& mom) const {
     183  if (verboseLevel > 3) {
     184    G4cout << " >>> G4LorentzConvertor::rotate(G4LorentzVector)" << G4endl
     185           << " ga " << ga << " gbpp " << gbpp << " gapp " << gapp << G4endl
     186           << " degenerated " << degenerated << G4endl
     187           << " before rotation: px " << mom.x() << " py " << mom.y()
     188           << " pz " << mom.z() << G4endl;
     189  }
     190
     191  G4LorentzVector mom_rot = mom;
     192  if (!degenerated) {
     193    G4ThreeVector vscm = velocity - gbpp*scm_momentum.vect();
     194    G4ThreeVector vxcm = scm_momentum.vect().cross(velocity);
     195
     196    mom_rot.setVect(mom.x()*vscm/ga + mom.y()*vxcm/gapp +
     197                    mom.z()*scm_momentum.vect().unit() );
    129198  };
    130199
    131200  if (verboseLevel > 3) {
    132     G4cout << " after rotation: px " << mom_rot[1] << " py " << mom_rot[2] <<
    133       " pz " << mom_rot[3] << G4endl;
     201    G4cout << " after rotation: px " << mom_rot.x() << " py " << mom_rot.y()
     202           << " pz " << mom_rot.z() << G4endl;
    134203  }
    135204
    … …  
    137206}
    138207
    139 G4CascadeMomentum G4LorentzConvertor::rotate(const G4CascadeMomentum& mom1,
    140                                             const G4CascadeMomentum& mom) const {
    141 
    142   if (verboseLevel > 3) {
    143     G4cout << " >>> G4LorentzConvertor::rotate(G4CascadeMomentum,G4CascadeMomentum)" << G4endl;
    144   }
    145 
    146   const G4double small = 1.0e-10;
    147 
    148   G4CascadeMomentum mom_rot;
    149 
    150   G4double pp = 0.0;
    151 
    152   G4double pv = 0.0;
    153 
    154   for(G4int i = 0; i < 4; i++) {
    155     pp += mom1[i] * mom1[i];
    156     pv += mom1[i] * velocity[i];
     208G4LorentzVector G4LorentzConvertor::rotate(const G4LorentzVector& mom1,
     209                                           const G4LorentzVector& mom) const {
     210  if (verboseLevel > 3) {
     211    G4cout << " >>> G4LorentzConvertor::rotate(G4LorentzVector,G4LorentzVector)"
     212           << G4endl
     213           << " before rotation: px " << mom.x() << " py " << mom.y()
     214           << " pz " << mom.z() << G4endl;
     215  }
     216
     217  G4double pp = mom1.vect().mag2();
     218  G4double pv = mom1.vect().dot(velocity);
     219
     220  G4double ga1 = v2 - pv * pv / pp;
     221  if (verboseLevel > 3) {
     222    G4cout << " ga1 " << ga1 << " small? " << (ga1 <= small) << G4endl;
     223  }
     224
     225  G4LorentzVector mom_rot = mom;
     226
     227  if (ga1 > small) {
     228    ga1 = std::sqrt(ga1);
     229
     230    G4double gb1 = pv / pp;
     231
     232    pp = std::sqrt(pp);
     233
     234    G4double ga1pp = ga1 * pp;
     235
     236    if (verboseLevel > 3) {
     237      G4cout << " gb1 " << gb1 << " ga1pp " << ga1pp << G4endl;
     238    }
     239
     240    G4ThreeVector vmom1 = velocity - gb1*mom1;
     241    G4ThreeVector vxm1  = mom1.vect().cross(velocity);
     242
     243    mom_rot.setVect(mom.x()*vmom1/ga1 + mom.y()*vxm1/ga1pp +
     244                    mom.z()*mom1.vect().unit() );
    157245  };
    158246
    159   G4double ga1 = v2 - pv * pv / pp;
    160 
    161   if(ga1 < small) {
    162     mom_rot = mom;
    163   } else { 
    164     ga1 = std::sqrt(ga1);
    165 
    166     G4double gb1 = pv / pp;
    167 
    168     pp = std::sqrt(pp);
    169 
    170     G4double ga1pp = ga1 * pp;
    171 
    172     mom_rot[1] = mom[1] * (velocity[1] - gb1 * mom1[1]) / ga1 +
    173       mom[2] * (mom1[2] * velocity[3] - mom1[3] * velocity[2]) / ga1pp +
    174       mom[3] * mom1[1] / pp;
    175     mom_rot[2] = mom[1] * (velocity[2] - gb1 * mom1[2]) / ga1 +
    176       mom[2] * (mom1[3] * velocity[1] - mom1[1] * velocity[3]) / ga1pp +
    177       mom[3] * mom1[2] / pp;
    178     mom_rot[3] = mom[1] * (velocity[3] - gb1 * mom1[3]) / ga1 +
    179       mom[2] * (mom1[1] * velocity[2] - mom1[2] * velocity[1]) / ga1pp +
    180       mom[3] * mom1[3] / pp;
    181   };
     247  if (verboseLevel > 3) {
     248    G4cout << " after rotation: px " << mom_rot.x() << " py " << mom_rot.y()
     249           << " pz " << mom_rot.z() << G4endl;
     250  }
    182251
    183252  return mom_rot;
    184253}
    185254
    186 void G4LorentzConvertor::toTheTargetRestFrame() {
    187    
    188   if (verboseLevel > 3) {
    189     G4cout << " >>> G4LorentzConvertor::toTheTargetRestFrame" << G4endl;
    190   }
    191 
    192   const G4double small = 1.0e-10;
    193 
    194   gamma = target_mom[0] / target_mass;
    195   v2 = 0.0;
    196 
    197   G4double pv = 0.0;
    198 
    199   //  G4double e_sum = target_mom[0] + bullet_mom[0];
    200 
    201   velocity.resize(4);
    202   G4int i(0);
    203   for(i = 1; i < 4; i++) {
    204     velocity[i] = target_mom[i] / target_mom[0];
    205     v2 += velocity[i] * velocity[i];
    206     pv += bullet_mom[i] * velocity[i];
    207   };
    208 
    209   G4double pa = 0.0;
    210 
    211   G4double pb = 0.0;
    212 
    213   G4double xx = 0.0;
    214 
    215   if(v2 > small) xx = pv * (gamma - 1.0) / v2 - bullet_mom[0] * gamma;
    216   for(i = 1; i < 4; i++) {
    217     scm_momentum[i] = bullet_mom[i] + velocity[i] * xx;
    218 
    219     if (verboseLevel > 3) {
    220       G4cout << " rf: i " << i << " pscm(i) " << scm_momentum[i] << G4endl;
    221     }
    222     pa += scm_momentum[i] * scm_momentum[i];
    223     pb += scm_momentum[i] * velocity[i];
    224   };
    225 
    226   ga = v2 - pb * pb / pa;
    227   if(ga < small) {
    228     ga = small;
    229     degenerated = true;
    230   } else {
    231     ga = std::sqrt(ga);
    232   }; 
    233   pscm = std::sqrt(pa);
    234   plab = pscm;
    235   gb = pb / pscm;
    236   gbpp = gb / pscm;
    237   gapp = ga * pscm;   
    238 }
    239 
    240 G4CascadeMomentum G4LorentzConvertor::backToTheLab(const G4CascadeMomentum& mom) const {
    241 
    242   if (verboseLevel > 3) {
    243     G4cout << " >>> G4LorentzConvertor::backToTheLab" << G4endl;
    244   }
    245 
    246   const G4double small = 1.0e-10;
    247 
    248   if (verboseLevel > 3) {
    249     G4cout << " at rest: px " << mom[1] << " py " << mom[2] << " pz " << mom[3] <<
    250       " e " << mom[0] << G4endl;
    251     G4cout << " v2 " << v2 << G4endl;   
    252   }
    253 
    254   G4CascadeMomentum mom1;
    255 
    256   if(v2 < small) {
    257     mom1 = mom;
    258   } else {
    259     G4double pv = 0.0;
    260 
    261     G4int i(0);
    262     for(i = 1; i < 4; i++) pv += mom[i] * velocity[i];
    263 
    264     G4double xx = pv * (gamma - 1.0) / v2 + mom[0] * gamma;
    265 
    266     for(i = 1; i < 4; i++) mom1[i] = mom[i] + velocity[i] * xx;
    267   };
    268 
    269   if (verboseLevel > 3) {
    270     G4cout << " at lab: px " << mom1[1] << " py " << mom1[2] << " pz " << mom1[3] << G4endl;
    271   }
    272 
    273   return mom1;
    274 }
    275 
    276255G4bool G4LorentzConvertor::reflectionNeeded() const {
    277 
    278   if (verboseLevel > 3) {
    279     G4cout << " >>> G4LorentzConvertor::reflectionNeeded" << G4endl;
    280   }
    281 
    282   const G4double small = 1.0e-10;
    283 
    284   if(v2 < small) {
    285     return false;
    286   }  else {   
    287     if(degenerated) return (scm_momentum[3] < 0.0);
    288     else
    289     {
    290       throw G4HadronicException(__FILE__, __LINE__, "G4LorentzConvertor::reflectionNeeded - return value undefined");
    291       return false;
    292     }
    293   };
    294 }
    295 
    296 
    297 
    298 
    299 
    300 
    301 
     256  if (verboseLevel > 3) {
     257    G4cout << " >>> G4LorentzConvertor::reflectionNeeded: v2 = " << v2
     258           << " SCM z = " << scm_momentum.z() << " degenerated? "
     259           << degenerated << G4endl;
     260  }
     261
     262  if (v2 < small && !degenerated)
     263    throw G4HadronicException(__FILE__, __LINE__, "G4LorentzConvertor::reflectionNeeded - return value undefined");
     264
     265  return (v2>=small && (!degenerated || scm_momentum.z() < 0.0));
     266}
     267
     268
     269
     270
     271
     272
     273
Note: See TracChangeset for help on using the changeset viewer.