G4LorentzConvertor.cc

Go to the documentation of this file.
00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00026 // $Id$
00027 //
00028 // 20100108  Michael Kelsey -- Use G4LorentzVector internally
00029 // 20100112  M. Kelsey -- Remove G4CascadeMomentum, use G4LorentzVector directly
00030 // 20100308  M. Kelsey -- Bug fix in toTheTargetRestFrame: scm_momentum should
00031 //              be data member, not local.
00032 // 20100409  M. Kelsey -- Protect std::sqrt(ga) against round-off negatives
00033 // 20100519  M. Kelsey -- Add interfaces to pass G4InuclParticles directly
00034 // 20100617  M. Kelsey -- Add more diagnostic messages with multiple levels
00035 // 20100713  M. Kelsey -- All diagnostic messages should be verbose > 1
00036 // 20110525  M. Kelsey -- Add some diagnostic messages in both ::rotate()
00037 // 20110602  M. Kelsey -- Simplify some kinematics, dropping intermediate calcs
00038 
00039 #include "G4LorentzConvertor.hh"
00040 #include "G4ThreeVector.hh"
00041 #include "G4HadronicException.hh"
00042 #include "G4InuclParticle.hh"
00043 
00044 
00045 const G4double G4LorentzConvertor::small = 1.0e-10;
00046 
00047 G4LorentzConvertor::G4LorentzConvertor() 
00048   : verboseLevel(0), v2(0.), ecm_tot(0.), valong(0.), degenerated(false) {}
00049 
00050 G4LorentzConvertor::
00051 G4LorentzConvertor(const G4LorentzVector& bmom, G4double bmass, 
00052                    const G4LorentzVector& tmom, G4double tmass) 
00053   : verboseLevel(0), v2(0.), ecm_tot(0.), valong(0.), degenerated(false) {
00054   setBullet(bmom, bmass);
00055   setTarget(tmom, tmass);
00056 }
00057 
00058 G4LorentzConvertor::
00059 G4LorentzConvertor(const G4InuclParticle* bullet, 
00060                    const G4InuclParticle* target) 
00061   : verboseLevel(0), v2(0.), ecm_tot(0.), valong(0.), degenerated(false) {
00062   setBullet(bullet);
00063   setTarget(target);
00064 }
00065 
00066 // Extract four-vectors from input particles
00067 void G4LorentzConvertor::setBullet(const G4InuclParticle* bullet) {
00068   setBullet(bullet->getMomentum());
00069 }
00070 
00071 void G4LorentzConvertor::setTarget(const G4InuclParticle* target) {
00072   setTarget(target->getMomentum());
00073 }
00074 
00075 
00076 // Boost bullet and target four-vectors into desired frame
00077 
00078 void G4LorentzConvertor::toTheCenterOfMass() {
00079   if (verboseLevel > 2)
00080     G4cout << " >>> G4LorentzConvertor::toTheCenterOfMass" << G4endl;
00081 
00082   velocity = (target_mom+bullet_mom).boostVector();
00083   if (verboseLevel > 3) G4cout << " boost " << velocity << G4endl;
00084 
00085   // "SCM" is reverse target momentum in the CM frame
00086   scm_momentum = target_mom;
00087   scm_momentum.boost(-velocity);
00088   scm_momentum.setVect(-scm_momentum.vect());
00089 
00090   if (verboseLevel > 3) G4cout << " pscm " << scm_momentum.vect() << G4endl;
00091 
00092   fillKinematics();
00093 }
00094 
00095 void G4LorentzConvertor::toTheTargetRestFrame() {
00096   if (verboseLevel > 2)
00097     G4cout << " >>> G4LorentzConvertor::toTheTargetRestFrame" << G4endl;
00098 
00099   velocity = target_mom.boostVector();
00100   if (verboseLevel > 3) G4cout << " boost " << velocity << G4endl;
00101 
00102   // "SCM" is bullet momentum in the target's frame
00103   scm_momentum = bullet_mom;
00104   scm_momentum.boost(-velocity);
00105 
00106   if (verboseLevel > 3) G4cout << " pseudo-pscm " << scm_momentum.vect() << G4endl;
00107 
00108   fillKinematics();
00109 }
00110 
00111 // Compute kinematic quantities for rotate() functions
00112 
00113 void G4LorentzConvertor::fillKinematics() {
00114   ecm_tot = (target_mom+bullet_mom).m();
00115 
00116   scm_direction = scm_momentum.vect().unit();
00117   valong = velocity.dot(scm_direction);
00118 
00119   v2 = velocity.mag2();
00120 
00121   G4double pvsq = v2 - valong*valong;           // velocity perp to scm_momentum
00122   if (verboseLevel > 3) G4cout << " pvsq " << pvsq << G4endl;
00123 
00124   degenerated = (pvsq < small);
00125   if (degenerated && verboseLevel > 2) 
00126     G4cout << " degenerated case (already along Z) " << G4endl; 
00127 
00128   if (verboseLevel > 3) {
00129     G4cout << " v2 " << v2 << " valong " << valong
00130            << " valong*valong " << valong*valong << G4endl;
00131   }
00132 }
00133 
00134 G4LorentzVector 
00135 G4LorentzConvertor::backToTheLab(const G4LorentzVector& mom) const {
00136   if (verboseLevel > 2)
00137     G4cout << " >>> G4LorentzConvertor::backToTheLab" << G4endl;
00138 
00139   if (verboseLevel > 3)
00140     G4cout << " at rest: px " << mom.x() << " py " << mom.y() << " pz "
00141            << mom.z() << " e " << mom.e() << G4endl
00142            << " v2 " << v2 << G4endl;
00143 
00144   G4LorentzVector mom1 = mom;
00145   if (v2 > small) mom1.boost(velocity);
00146 
00147   if (verboseLevel > 3)
00148     G4cout << " at lab: px " << mom1.x() << " py " << mom1.y() << " pz "
00149            << mom1.z() << G4endl;
00150 
00151   return mom1;
00152 }
00153 
00154 
00155 // Bullet kinematics in target rest frame (LAB frame, usually)
00156 
00157 G4double G4LorentzConvertor::getKinEnergyInTheTRS() const {
00158   if (verboseLevel > 2)
00159     G4cout << " >>> G4LorentzConvertor::getKinEnergyInTheTRS" << G4endl;
00160 
00161   G4LorentzVector bmom = bullet_mom;
00162   bmom.boost(-target_mom.boostVector());
00163   return bmom.e()-bmom.m();
00164 }
00165 
00166 G4double G4LorentzConvertor::getTRSMomentum() const {
00167   if (verboseLevel > 2)
00168     G4cout << " >>> G4LorentzConvertor::getTRSMomentum" << G4endl;
00169 
00170   G4LorentzVector bmom = bullet_mom;
00171   bmom.boost(-target_mom.boostVector());
00172   return bmom.rho();
00173 }
00174 
00175 G4LorentzVector G4LorentzConvertor::rotate(const G4LorentzVector& mom) const {
00176   if (verboseLevel > 2)
00177     G4cout << " >>> G4LorentzConvertor::rotate(G4LorentzVector)" << G4endl;
00178 
00179   if (verboseLevel > 3) {
00180     G4cout << " valong " << valong << " degenerated " << degenerated << G4endl
00181            << " before rotation: px " << mom.x() << " py " << mom.y()
00182            << " pz " << mom.z() << G4endl;
00183   }
00184 
00185   G4LorentzVector mom_rot = mom;
00186   if (!degenerated) {
00187     if (verboseLevel > 2)
00188       G4cout << " rotating to align with reference z axis " << G4endl;
00189 
00190     G4ThreeVector vscm = velocity - valong*scm_direction;
00191     G4ThreeVector vxcm = scm_direction.cross(velocity);
00192 
00193     if (vscm.mag() > small && vxcm.mag() > small) {     // Double check
00194       if (verboseLevel > 3) {
00195         G4cout << " reference z axis " << scm_direction
00196                << " vscm " << vscm << " vxcm " << vxcm << G4endl;
00197       }
00198       
00199       mom_rot.setVect(mom.x()*vscm.unit() + mom.y()*vxcm.unit() +
00200                       mom.z()*scm_direction);
00201     } else {
00202       if (verboseLevel) 
00203         G4cerr << ">>> G4LorentzVector::rotate zero with !degenerated" << G4endl;
00204     }
00205   }
00206 
00207   if (verboseLevel > 3) {
00208     G4cout << " after rotation: px " << mom_rot.x() << " py " << mom_rot.y()
00209            << " pz " << mom_rot.z() << G4endl;
00210   }
00211 
00212   return mom_rot;
00213 }
00214 
00215 G4LorentzVector G4LorentzConvertor::rotate(const G4LorentzVector& mom1, 
00216                                            const G4LorentzVector& mom) const {
00217   if (verboseLevel > 2)
00218     G4cout << " >>> G4LorentzConvertor::rotate(G4LorentzVector,G4LorentzVector)"
00219            << G4endl;
00220 
00221   if (verboseLevel > 3) {
00222     G4cout << " before rotation: px " << mom.x() << " py " << mom.y()
00223            << " pz " << mom.z() << G4endl;
00224   }
00225 
00226   G4ThreeVector mom1_dir = mom1.vect().unit();
00227   G4double pv = velocity.dot(mom1_dir);
00228 
00229   G4double vperp = v2 - pv*pv;          // velocity perpendicular to mom1
00230   if (verboseLevel > 3) {
00231     G4cout << " vperp " << vperp << " small? " << (vperp <= small) << G4endl;
00232   }
00233 
00234   G4LorentzVector mom_rot = mom;
00235 
00236   if (vperp > small) {
00237     if (verboseLevel > 2)
00238       G4cout << " rotating to align with first z axis " << G4endl;
00239 
00240     G4ThreeVector vmom1 = velocity - pv*mom1_dir;
00241     G4ThreeVector vxm1  = mom1_dir.cross(velocity);
00242 
00243     if (vmom1.mag() > small && vxm1.mag() > small) {    // Double check
00244       if (verboseLevel > 3) {
00245         G4cout << " first z axis " << mom1_dir << G4endl
00246                << " vmom1 " << vmom1 << " vxm1 " << vxm1 << G4endl;
00247       }
00248       
00249       mom_rot.setVect(mom.x()*vmom1.unit() + mom.y()*vxm1.unit() +
00250                       mom.z()*mom1_dir );
00251     } else {
00252       if (verboseLevel) 
00253         G4cerr << ">>> G4LorentzVector::rotate zero with !degenerated" << G4endl;
00254     }
00255   }
00256 
00257   if (verboseLevel > 3) {
00258     G4cout << " after rotation: px " << mom_rot.x() << " py " << mom_rot.y()
00259            << " pz " << mom_rot.z() << G4endl;
00260   }
00261 
00262   return mom_rot;
00263 }
00264 
00265 G4bool G4LorentzConvertor::reflectionNeeded() const {
00266   if (verboseLevel > 2)
00267     G4cout << " >>> G4LorentzConvertor::reflectionNeeded (query)" << G4endl;
00268 
00269   if (verboseLevel > 3) {
00270     G4cout << " v2 = " << v2 << " SCM z = " << scm_momentum.z()
00271            << " degenerated? " << degenerated << G4endl;
00272   }
00273 
00274   if (v2 < small && !degenerated) 
00275     throw G4HadronicException(__FILE__, __LINE__, "G4LorentzConvertor::reflectionNeeded - return value undefined");
00276 
00277   if (verboseLevel > 2) {
00278     G4cout << " reflection across XY is"
00279            << ((v2>=small && (!degenerated || scm_momentum.z()<0.0))?"":" NOT")
00280            << " needed" << G4endl;
00281   }
00282 
00283   return (v2>=small && (!degenerated || scm_momentum.z()<0.0));
00284 }
00285 
00286 
00287 // Reporting functions for diagnostics
00288 
00289 void G4LorentzConvertor::printBullet() const {
00290   G4cout << " G4LC bullet: px " << bullet_mom.px() << " py " << bullet_mom.py()
00291          << " pz " << bullet_mom.pz() << " e " << bullet_mom.e()
00292          << " mass " << bullet_mom.m() << G4endl;
00293   }
00294 
00295 void G4LorentzConvertor::printTarget() const {
00296   G4cout << " G4LC target: px " << target_mom.px() << " py " << target_mom.py()
00297          << " pz " << target_mom.pz() << " e " << target_mom.e()
00298          << " mass " << target_mom.m() << G4endl;
00299 }
00300 

Generated on Mon May 27 17:48:49 2013 for Geant4 by  doxygen 1.4.7