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 "G4Absorber.hh" 26 #include "G4Absorber.hh" 27 #include "G4KineticTrack.hh" 27 #include "G4KineticTrack.hh" 28 #include "G4PionPlus.hh" 28 #include "G4PionPlus.hh" 29 #include "G4PionMinus.hh" 29 #include "G4PionMinus.hh" 30 #include "G4PionZero.hh" 30 #include "G4PionZero.hh" 31 #include "G4Proton.hh" 31 #include "G4Proton.hh" 32 #include "G4Neutron.hh" 32 #include "G4Neutron.hh" 33 33 34 #include "G4PhysicalConstants.hh" << 35 #include "G4SystemOfUnits.hh" << 36 #include "G4LorentzRotation.hh" 34 #include "G4LorentzRotation.hh" 37 35 38 G4Absorber::G4Absorber(G4double cutOnP) 36 G4Absorber::G4Absorber(G4double cutOnP) 39 { 37 { 40 theCutOnP = cutOnP; 38 theCutOnP = cutOnP; 41 theAbsorbers = new G4KineticTrackVector; 39 theAbsorbers = new G4KineticTrackVector; 42 theProducts = new G4KineticTrackVector; 40 theProducts = new G4KineticTrackVector; 43 } 41 } 44 42 45 43 46 G4Absorber::~G4Absorber() 44 G4Absorber::~G4Absorber() 47 { 45 { 48 delete theAbsorbers; 46 delete theAbsorbers; 49 delete theProducts; 47 delete theProducts; 50 } 48 } 51 49 52 50 53 bool G4Absorber::WillBeAbsorbed(const G4Kineti 51 bool G4Absorber::WillBeAbsorbed(const G4KineticTrack & kt) 54 { 52 { 55 // FixMe: actually only for pions 53 // FixMe: actually only for pions 56 // if(kt.Get4Momentum().vect().mag() < theCut 54 // if(kt.Get4Momentum().vect().mag() < theCutOnP) 57 // Cut on kinetic Energy... 55 // Cut on kinetic Energy... 58 if (kt.Get4Momentum().e() - kt.GetActualMass 56 if (kt.Get4Momentum().e() - kt.GetActualMass() < theCutOnP) 59 { 57 { 60 if(kt.GetDefinition() == G4PionPlus::Pio 58 if(kt.GetDefinition() == G4PionPlus::PionPlus() || 61 kt.GetDefinition() == G4PionZero::PionZero( 59 kt.GetDefinition() == G4PionZero::PionZero() || 62 kt.GetDefinition() == G4PionMinus::PionMinu 60 kt.GetDefinition() == G4PionMinus::PionMinus()) 63 { 61 { 64 return true; 62 return true; 65 } 63 } 66 } 64 } 67 return false; 65 return false; 68 } 66 } 69 67 70 68 71 69 72 G4bool G4Absorber::Absorb(G4KineticTrack & kt, 70 G4bool G4Absorber::Absorb(G4KineticTrack & kt, G4KineticTrackVector & tgt) 73 { 71 { 74 if(!FindAbsorbers(kt, tgt)) 72 if(!FindAbsorbers(kt, tgt)) 75 return false; 73 return false; 76 return FindProducts(kt); 74 return FindProducts(kt); 77 } 75 } 78 76 79 77 80 G4bool G4Absorber::FindAbsorbers(G4KineticTrac 78 G4bool G4Absorber::FindAbsorbers(G4KineticTrack & kt, 81 G4KineticTrackVector & tgt) 79 G4KineticTrackVector & tgt) 82 { 80 { 83 // Find a closest ( in space) pair of Nucleon << 84 // pi+ can be absorbed on np or nn resultin << 85 // pi- can be absorbed on np or pp resultin << 86 << 87 // @GF: FindAbsorbers is unused, logic is seri << 88 << 89 G4KineticTrack * kt1 = NULL; 81 G4KineticTrack * kt1 = NULL; 90 G4KineticTrack * kt2 = NULL; 82 G4KineticTrack * kt2 = NULL; 91 G4double dist1 = DBL_MAX; // dist to close << 83 G4double dist1 = DBL_MAX; 92 G4double dist2 = DBL_MAX; // dist to next << 84 G4double dist2 = DBL_MAX; 93 G4double charge1 = 0; 85 G4double charge1 = 0; 94 // G4double charge2 = 0; // charge2 is only << 86 G4double charge2 = 0; 95 G4double charge0 = kt.GetDefinition()->GetPD 87 G4double charge0 = kt.GetDefinition()->GetPDGCharge(); 96 G4ThreeVector pos = kt.GetPosition(); 88 G4ThreeVector pos = kt.GetPosition(); 97 89 98 std::vector<G4KineticTrack *>::iterator iter 90 std::vector<G4KineticTrack *>::iterator iter; 99 for(iter = tgt.begin(); iter != tgt.end(); + 91 for(iter = tgt.begin(); iter != tgt.end(); ++iter) 100 { 92 { 101 G4KineticTrack * curr = *iter; 93 G4KineticTrack * curr = *iter; 102 G4double dist = (pos-curr->GetPosition()). 94 G4double dist = (pos-curr->GetPosition()).mag(); 103 if(dist >= dist2) 95 if(dist >= dist2) 104 continue; 96 continue; 105 if(dist < dist1) 97 if(dist < dist1) 106 { 98 { 107 if(dist1 == DBL_MAX) // accept 1st as a << 99 if(dist1 == DBL_MAX) // accept the candidate 108 { 100 { 109 kt1 = curr; 101 kt1 = curr; 110 charge1 = kt1->GetDefinition()->GetPDGCharge 102 charge1 = kt1->GetDefinition()->GetPDGCharge(); 111 dist1 = dist; 103 dist1 = dist; 112 continue; 104 continue; 113 } 105 } 114 if(dist2 == DBL_MAX) // accept the candi << 106 if(dist2 == DBL_MAX) // accept the candidate put kt1 in kt2 115 { // @GF: should'nt we check if co << 107 { 116 kt2 = kt1; 108 kt2 = kt1; 117 // charge2 = charge1; << 109 charge2 = charge1; 118 dist2 = dist1; 110 dist2 = dist1; 119 kt1 = curr; 111 kt1 = curr; 120 charge1 = kt1->GetDefinition()->GetPDGCharge 112 charge1 = kt1->GetDefinition()->GetPDGCharge(); 121 dist1 = dist; 113 dist1 = dist; 122 continue; 114 continue; 123 } 115 } 124 // test the compatibility with charge conserva << 116 // test the compatibility with charge conservation 125 G4double charge = curr->GetDefinition()- 117 G4double charge = curr->GetDefinition()->GetPDGCharge(); 126 if((charge0+charge1+charge < 0.) || //te << 118 if((charge0+charge1+charge < 0.) || 127 (charge0+charge1+charge) > 2*eplus) 119 (charge0+charge1+charge) > 2*eplus) 128 { // incompatible: change kt1 with curr << 120 { // incomatible: change kt1 with curr. 129 kt1 = curr; 121 kt1 = curr; 130 charge1 = charge; 122 charge1 = charge; 131 dist1 = dist; 123 dist1 = dist; 132 } 124 } 133 else 125 else 134 { // compatible: change kt1 with curr an 126 { // compatible: change kt1 with curr and kt2 with kt1 135 kt2 = kt1; 127 kt2 = kt1; 136 // charge2 = charge1; << 128 charge2 = charge1; 137 dist2 = dist1; 129 dist2 = dist1; 138 kt1 = curr; 130 kt1 = curr; 139 charge1 = charge; 131 charge1 = charge; 140 dist1 = dist; 132 dist1 = dist; 141 } 133 } 142 continue; 134 continue; 143 } 135 } 144 // here if dist1 < dist < dist2 136 // here if dist1 < dist < dist2 145 if(dist2 == DBL_MAX) // accept the candida 137 if(dist2 == DBL_MAX) // accept the candidate 146 { 138 { 147 kt2 = curr; 139 kt2 = curr; 148 // charge2 = kt2->GetDefinition()->GetPDG << 140 charge2 = kt2->GetDefinition()->GetPDGCharge(); 149 dist2 = dist; 141 dist2 = dist; 150 continue; 142 continue; 151 } 143 } 152 // test the compatibility with charge conserva 144 // test the compatibility with charge conservation 153 G4double charge = curr->GetDefinition()->G 145 G4double charge = curr->GetDefinition()->GetPDGCharge(); 154 if((charge0+charge1+charge < 0.) || 146 if((charge0+charge1+charge < 0.) || 155 (charge0+charge1+charge) > 2*eplus) 147 (charge0+charge1+charge) > 2*eplus) 156 continue; // incomatible: do nothing 148 continue; // incomatible: do nothing 157 // compatible: change kt2 with curr 149 // compatible: change kt2 with curr 158 kt2 = curr; 150 kt2 = curr; 159 // charge2 = charge; << 151 charge2 = charge; 160 dist2 = dist; 152 dist2 = dist; 161 } 153 } 162 154 163 theAbsorbers->clear(); // do not delete trac 155 theAbsorbers->clear(); // do not delete tracks in theAbsorbers vector! 164 if((kt1 == NULL) || (kt2 == NULL)) 156 if((kt1 == NULL) || (kt2 == NULL)) 165 return false; 157 return false; 166 158 167 theAbsorbers->push_back(kt1); 159 theAbsorbers->push_back(kt1); 168 theAbsorbers->push_back(kt2); 160 theAbsorbers->push_back(kt2); 169 return true; 161 return true; 170 } 162 } 171 163 172 164 173 165 174 G4bool G4Absorber::FindProducts(G4KineticTrack 166 G4bool G4Absorber::FindProducts(G4KineticTrack & kt) 175 { 167 { 176 // Choose the products type 168 // Choose the products type 177 const G4ParticleDefinition * prod1; << 169 G4ParticleDefinition * prod1; 178 const G4ParticleDefinition * prod2; << 170 G4ParticleDefinition * prod2; 179 G4KineticTrack * abs1 = (*theAbsorbers)[0]; 171 G4KineticTrack * abs1 = (*theAbsorbers)[0]; 180 G4KineticTrack * abs2 = (*theAbsorbers)[1]; 172 G4KineticTrack * abs2 = (*theAbsorbers)[1]; 181 173 182 G4double charge = kt.GetDefinition()->GetPDG 174 G4double charge = kt.GetDefinition()->GetPDGCharge(); 183 if(charge == eplus) 175 if(charge == eplus) 184 { // a neutron become proton 176 { // a neutron become proton 185 prod1 = G4Proton::Proton(); 177 prod1 = G4Proton::Proton(); 186 if(abs1->GetDefinition() == G4Neutron::Neu 178 if(abs1->GetDefinition() == G4Neutron::Neutron()) 187 prod2 = abs2->GetDefinition(); 179 prod2 = abs2->GetDefinition(); 188 else 180 else 189 prod2 = G4Proton::Proton(); 181 prod2 = G4Proton::Proton(); 190 } 182 } 191 else if(charge == -eplus) 183 else if(charge == -eplus) 192 { // a proton become neutron 184 { // a proton become neutron 193 prod1 = G4Neutron::Neutron(); 185 prod1 = G4Neutron::Neutron(); 194 if(abs1->GetDefinition() == G4Proton::Prot 186 if(abs1->GetDefinition() == G4Proton::Proton()) 195 prod2 = abs2->GetDefinition(); 187 prod2 = abs2->GetDefinition(); 196 else 188 else 197 prod2 = G4Neutron::Neutron(); 189 prod2 = G4Neutron::Neutron(); 198 } 190 } 199 else // charge = 0: leave particle types un 191 else // charge = 0: leave particle types unchenged 200 { 192 { 201 prod1 = abs1->GetDefinition(); 193 prod1 = abs1->GetDefinition(); 202 prod2 = abs2->GetDefinition(); 194 prod2 = abs2->GetDefinition(); 203 } 195 } 204 196 205 // Translate to the CMS frame 197 // Translate to the CMS frame 206 G4LorentzVector momLab = kt.Get4Momentum()+a 198 G4LorentzVector momLab = kt.Get4Momentum()+abs1->Get4Momentum()+ 207 abs2->Get4Momentum(); 199 abs2->Get4Momentum(); 208 G4LorentzRotation toCMSFrame((-1)*momLab.boo 200 G4LorentzRotation toCMSFrame((-1)*momLab.boostVector()); 209 G4LorentzRotation toLabFrame(momLab.boostVec 201 G4LorentzRotation toLabFrame(momLab.boostVector()); 210 G4LorentzVector momCMS = toCMSFrame*momLab; 202 G4LorentzVector momCMS = toCMSFrame*momLab; 211 203 212 // Evaluate the final momentum of products 204 // Evaluate the final momentum of products 213 G4double ms1 = prod1->GetPDGMass(); << 205 G4double m1 = prod1->GetPDGMass(); 214 G4double ms2 = prod2->GetPDGMass(); << 206 G4double m2 = prod2->GetPDGMass(); 215 G4double e0 = momCMS.e(); 207 G4double e0 = momCMS.e(); 216 G4double squareP = (e0*e0*e0*e0-2*e0*e0*(ms1 << 208 G4double squareP = (e0*e0*e0*e0-2*e0*e0*(m1*m1+m2*m2)+ 217 (ms2*ms2-ms1*ms1)*(ms2*ms2-ms1*ms1))/(4*e0 << 209 (m2*m2-m1*m1)*(m2*m2-m1*m1))/(4*e0*e0); 218 // if(squareP < 0) // should never happen 210 // if(squareP < 0) // should never happen 219 // squareP = 0; 211 // squareP = 0; 220 G4ThreeVector mom1CMS = GetRandomDirection() 212 G4ThreeVector mom1CMS = GetRandomDirection(); 221 mom1CMS = std::sqrt(squareP)*mom1CMS; 213 mom1CMS = std::sqrt(squareP)*mom1CMS; 222 G4LorentzVector final4Mom1CMS(mom1CMS, std:: << 214 G4LorentzVector final4Mom1CMS(mom1CMS, std::sqrt(squareP+m1*m1)); 223 G4LorentzVector final4Mom2CMS((-1)*mom1CMS, << 215 G4LorentzVector final4Mom2CMS((-1)*mom1CMS, std::sqrt(squareP+m2*m2)); 224 216 225 // Go back to the lab frame 217 // Go back to the lab frame 226 G4LorentzVector mom1 = toLabFrame*final4Mom1 218 G4LorentzVector mom1 = toLabFrame*final4Mom1CMS; 227 G4LorentzVector mom2 = toLabFrame*final4Mom2 219 G4LorentzVector mom2 = toLabFrame*final4Mom2CMS; 228 220 229 // ------ debug 221 // ------ debug 230 /* 222 /* 231 G4LorentzVector temp = mom1+mom2; 223 G4LorentzVector temp = mom1+mom2; 232 224 233 cout << (1/MeV)*momLab.x() << " " << (1/MeV) 225 cout << (1/MeV)*momLab.x() << " " << (1/MeV)*momLab.y() << " " 234 << (1/MeV)*momLab.z() << " " << (1/MeV) 226 << (1/MeV)*momLab.z() << " " << (1/MeV)*momLab.t() << " " 235 << (1/MeV)*momLab.vect().mag() << " " < 227 << (1/MeV)*momLab.vect().mag() << " " << (1/MeV)*momLab.mag() << " " 236 << (1/MeV)*temp.x() << " " << (1/MeV)*t 228 << (1/MeV)*temp.x() << " " << (1/MeV)*temp.y() << " " 237 << (1/MeV)*temp.z() << " " << (1/MeV)*t 229 << (1/MeV)*temp.z() << " " << (1/MeV)*temp.t() << " " 238 << (1/MeV)*temp.vect().mag() << " " << 230 << (1/MeV)*temp.vect().mag() << " " << (1/MeV)*temp.mag() << " " 239 << (1/MeV)*std::sqrt(squareP) << endl; 231 << (1/MeV)*std::sqrt(squareP) << endl; 240 232 241 */ 233 */ 242 // ------ end debug 234 // ------ end debug 243 235 244 // Build two new kinetic tracks and add to pro 236 // Build two new kinetic tracks and add to products 245 G4KineticTrack * kt1 = new G4KineticTrack(pr 237 G4KineticTrack * kt1 = new G4KineticTrack(prod1, 0., abs1->GetPosition(), 246 mom1); 238 mom1); 247 G4KineticTrack * kt2 = new G4KineticTrack(pr 239 G4KineticTrack * kt2 = new G4KineticTrack(prod2, 0., abs2->GetPosition(), 248 mom2); 240 mom2); 249 // ------ debug 241 // ------ debug 250 /* 242 /* 251 G4LorentzVector initialMom1 = abs1->Get4Mome 243 G4LorentzVector initialMom1 = abs1->Get4Momentum(); 252 G4LorentzVector initialMom2 = abs2->Get4Mome 244 G4LorentzVector initialMom2 = abs2->Get4Momentum(); 253 G4LorentzVector pion4MomCMS = toCMSFrame*kt. 245 G4LorentzVector pion4MomCMS = toCMSFrame*kt.Get4Momentum(); 254 cout << (1/MeV)*initialMom1.x() << " " << (1 246 cout << (1/MeV)*initialMom1.x() << " " << (1/MeV)*initialMom1.y() << " " 255 << (1/MeV)*initialMom1.z() << " " << (1 247 << (1/MeV)*initialMom1.z() << " " << (1/MeV)*initialMom1.e() << " " 256 << (1/MeV)*initialMom1.vect().mag() << 248 << (1/MeV)*initialMom1.vect().mag() << " " 257 << (1/MeV)*initialMom2.x() << " " << (1 249 << (1/MeV)*initialMom2.x() << " " << (1/MeV)*initialMom2.y() << " " 258 << (1/MeV)*initialMom2.z() << " " << (1 250 << (1/MeV)*initialMom2.z() << " " << (1/MeV)*initialMom2.e() << " " 259 << (1/MeV)*initialMom2.vect().mag() << 251 << (1/MeV)*initialMom2.vect().mag() << " " 260 << (1/MeV)*mom1.x() << " " << (1/MeV)*m 252 << (1/MeV)*mom1.x() << " " << (1/MeV)*mom1.y() << " " 261 << (1/MeV)*mom1.z() << " " << (1/MeV)*m 253 << (1/MeV)*mom1.z() << " " << (1/MeV)*mom1.e() << " " 262 << (1/MeV)*mom1.vect().mag() << " " 254 << (1/MeV)*mom1.vect().mag() << " " 263 << (1/MeV)*mom2.x() << " " << (1/MeV)*m 255 << (1/MeV)*mom2.x() << " " << (1/MeV)*mom2.y() << " " 264 << (1/MeV)*mom2.z() << " " << (1/MeV)*m 256 << (1/MeV)*mom2.z() << " " << (1/MeV)*mom2.e() << " " 265 << (1/MeV)*mom2.vect().mag() << " " 257 << (1/MeV)*mom2.vect().mag() << " " 266 << (1/MeV)*pion4MomCMS.x() << " " << (1 258 << (1/MeV)*pion4MomCMS.x() << " " << (1/MeV)*pion4MomCMS.y() << " " 267 << (1/MeV)*pion4MomCMS.z() << " " << (1 259 << (1/MeV)*pion4MomCMS.z() << " " << (1/MeV)*pion4MomCMS.e() << " " 268 << (1/MeV)*pion4MomCMS.vect().mag() << 260 << (1/MeV)*pion4MomCMS.vect().mag() << " " 269 << (1/MeV)*final4Mom1CMS.x() << " " << 261 << (1/MeV)*final4Mom1CMS.x() << " " << (1/MeV)*final4Mom1CMS.y() << " " 270 << (1/MeV)*final4Mom1CMS.z() << " " << 262 << (1/MeV)*final4Mom1CMS.z() << " " << (1/MeV)*final4Mom1CMS.e() << " " 271 << (1/MeV)*final4Mom1CMS.vect().mag() < 263 << (1/MeV)*final4Mom1CMS.vect().mag() << " " 272 << (1/MeV)*final4Mom2CMS.x() << " " << 264 << (1/MeV)*final4Mom2CMS.x() << " " << (1/MeV)*final4Mom2CMS.y() << " " 273 << (1/MeV)*final4Mom2CMS.z() << " " << 265 << (1/MeV)*final4Mom2CMS.z() << " " << (1/MeV)*final4Mom2CMS.e() << " " 274 << (1/MeV)*final4Mom2CMS.vect().mag() < 266 << (1/MeV)*final4Mom2CMS.vect().mag() << endl; 275 */ 267 */ 276 // ------ end debug 268 // ------ end debug 277 269 278 theProducts->clear(); 270 theProducts->clear(); 279 theProducts->push_back(kt1); 271 theProducts->push_back(kt1); 280 theProducts->push_back(kt2); 272 theProducts->push_back(kt2); 281 return true; 273 return true; 282 } 274 } 283 275 284 276 285 277 286 G4ThreeVector G4Absorber::GetRandomDirection() 278 G4ThreeVector G4Absorber::GetRandomDirection() 287 { 279 { 288 G4double theta = 2.0*G4UniformRand()-1.0; 280 G4double theta = 2.0*G4UniformRand()-1.0; 289 theta = std::acos(theta); 281 theta = std::acos(theta); 290 G4double phi = G4UniformRand()*2*pi; 282 G4double phi = G4UniformRand()*2*pi; 291 G4ThreeVector direction(std::sin(theta)*std: 283 G4ThreeVector direction(std::sin(theta)*std::cos(phi), std::sin(theta)*std::sin(phi), std::cos(theta)); 292 return direction; 284 return direction; 293 } 285 } 294 286 295 287 296 288 297 289 298 290 299 291 300 292 301 293