[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 "G4MesonAbsorption.hh" |
---|
| 27 | #include "G4LorentzRotation.hh" |
---|
| 28 | #include "G4LorentzVector.hh" |
---|
| 29 | #include "Randomize.hh" |
---|
| 30 | #include "G4KineticTrackVector.hh" |
---|
| 31 | #include "G4CollisionInitialState.hh" |
---|
| 32 | #include <cmath> |
---|
| 33 | #include "G4PionPlus.hh" |
---|
| 34 | #include "G4PionMinus.hh" |
---|
| 35 | #include "G4ParticleDefinition.hh" |
---|
| 36 | #include "G4HadTmpUtil.hh" |
---|
| 37 | |
---|
| 38 | // first prototype |
---|
| 39 | |
---|
| 40 | const std::vector<G4CollisionInitialState *> & G4MesonAbsorption:: |
---|
| 41 | GetCollisions(G4KineticTrack * aProjectile, |
---|
| 42 | std::vector<G4KineticTrack *> & someCandidates, |
---|
| 43 | G4double aCurrentTime) |
---|
| 44 | { |
---|
| 45 | theCollisions.clear(); |
---|
| 46 | if(someCandidates.size() >1) |
---|
| 47 | { |
---|
| 48 | std::vector<G4KineticTrack *>::iterator j=someCandidates.begin(); |
---|
| 49 | for(; j != someCandidates.end(); ++j) |
---|
| 50 | { |
---|
| 51 | G4double collisionTime = GetTimeToAbsorption(*aProjectile, **j); |
---|
| 52 | if(collisionTime == DBL_MAX) |
---|
| 53 | { |
---|
| 54 | continue; |
---|
| 55 | } |
---|
| 56 | G4KineticTrackVector aTarget; |
---|
| 57 | aTarget.push_back(*j); |
---|
| 58 | FindAndFillCluster(aTarget, aProjectile, someCandidates); |
---|
| 59 | if(aTarget.size()>=2) |
---|
| 60 | { |
---|
| 61 | theCollisions.push_back( |
---|
| 62 | new G4CollisionInitialState(collisionTime+aCurrentTime, aProjectile, aTarget, this) ); |
---|
| 63 | } |
---|
| 64 | } |
---|
| 65 | } |
---|
| 66 | return theCollisions; |
---|
| 67 | } |
---|
| 68 | |
---|
| 69 | |
---|
| 70 | void G4MesonAbsorption:: |
---|
| 71 | FindAndFillCluster(G4KineticTrackVector & result, |
---|
| 72 | G4KineticTrack * aProjectile, std::vector<G4KineticTrack *> & someCandidates) |
---|
| 73 | { |
---|
| 74 | std::vector<G4KineticTrack *>::iterator j=someCandidates.begin(); |
---|
| 75 | G4KineticTrack * aTarget = result[0]; |
---|
| 76 | G4int chargeSum = G4lrint(aTarget->GetDefinition()->GetPDGCharge()); |
---|
| 77 | chargeSum+=G4lrint(aProjectile->GetDefinition()->GetPDGCharge()); |
---|
| 78 | G4ThreeVector firstBase = aTarget->GetPosition(); |
---|
| 79 | G4double min = DBL_MAX; |
---|
| 80 | G4KineticTrack * partner = 0; |
---|
| 81 | for(; j != someCandidates.end(); ++j) |
---|
| 82 | { |
---|
| 83 | if(*j == aTarget) continue; |
---|
| 84 | G4int cCharge = G4lrint((*j)->GetDefinition()->GetPDGCharge()); |
---|
| 85 | if (chargeSum+cCharge > 2) continue; |
---|
| 86 | if (chargeSum+cCharge < 0) continue; |
---|
| 87 | // get the one with the smallest distance. |
---|
| 88 | G4ThreeVector secodeBase = (*j)->GetPosition(); |
---|
| 89 | if((firstBase+secodeBase).mag()<min) |
---|
| 90 | { |
---|
| 91 | min=(firstBase+secodeBase).mag(); |
---|
| 92 | partner = *j; |
---|
| 93 | } |
---|
| 94 | } |
---|
| 95 | if(partner) result.push_back(partner); |
---|
| 96 | else result.clear(); |
---|
| 97 | } |
---|
| 98 | |
---|
| 99 | |
---|
| 100 | G4KineticTrackVector * G4MesonAbsorption:: |
---|
| 101 | GetFinalState(G4KineticTrack * projectile, |
---|
| 102 | std::vector<G4KineticTrack *> & targets) |
---|
| 103 | { |
---|
| 104 | // G4cout << "We have an absorption !!!!!!!!!!!!!!!!!!!!!!"<<G4endl; |
---|
| 105 | // Only 2-body absorption for the time being. |
---|
| 106 | // If insufficient, use 2-body absorption and clusterization to emulate 3-,4-,etc body absorption. |
---|
| 107 | G4LorentzVector thePro = projectile->Get4Momentum(); |
---|
| 108 | G4LorentzVector theT1 = targets[0]->Get4Momentum(); |
---|
| 109 | G4LorentzVector theT2 = targets[1]->Get4Momentum(); |
---|
| 110 | G4LorentzVector incoming = thePro + theT1 + theT2; |
---|
| 111 | G4double energyBalance = incoming.t(); |
---|
| 112 | G4int chargeBalance = G4lrint(projectile->GetDefinition()->GetPDGCharge() |
---|
| 113 | + targets[0]->GetDefinition()->GetPDGCharge() |
---|
| 114 | + targets[1]->GetDefinition()->GetPDGCharge()); |
---|
| 115 | |
---|
| 116 | G4int baryonBalance = projectile->GetDefinition()->GetBaryonNumber() |
---|
| 117 | + targets[0]->GetDefinition()->GetBaryonNumber() |
---|
| 118 | + targets[1]->GetDefinition()->GetBaryonNumber(); |
---|
| 119 | |
---|
| 120 | |
---|
| 121 | // boost all to MMS |
---|
| 122 | G4LorentzRotation toSPS((-1)*(thePro + theT1 + theT2).boostVector()); |
---|
| 123 | theT1 = toSPS * theT1; |
---|
| 124 | theT2 = toSPS * theT2; |
---|
| 125 | thePro = toSPS * thePro; |
---|
| 126 | G4LorentzRotation fromSPS(toSPS.inverse()); |
---|
| 127 | |
---|
| 128 | // rotate to z |
---|
| 129 | G4LorentzRotation toZ; |
---|
| 130 | G4LorentzVector Ptmp=projectile->Get4Momentum(); |
---|
| 131 | toZ.rotateZ(-1*Ptmp.phi()); |
---|
| 132 | toZ.rotateY(-1*Ptmp.theta()); |
---|
| 133 | theT1 = toZ * theT1; |
---|
| 134 | theT2 = toZ * theT2; |
---|
| 135 | thePro = toZ * thePro; |
---|
| 136 | G4LorentzRotation toLab(toZ.inverse()); |
---|
| 137 | |
---|
| 138 | // Get definitions |
---|
| 139 | G4ParticleDefinition * d1 = targets[0]->GetDefinition(); |
---|
| 140 | G4ParticleDefinition * d2 = targets[1]->GetDefinition(); |
---|
| 141 | if(0.5>G4UniformRand()) |
---|
| 142 | { |
---|
| 143 | G4ParticleDefinition * temp; |
---|
| 144 | temp=d1;d1=d2;d2=temp; |
---|
| 145 | } |
---|
| 146 | G4ParticleDefinition * dp = projectile->GetDefinition(); |
---|
| 147 | if(dp->GetPDGCharge()<-.5) |
---|
| 148 | { |
---|
| 149 | if(d1->GetPDGCharge()>.5) |
---|
| 150 | { |
---|
| 151 | if(d2->GetPDGCharge()>.5 && 0.5>G4UniformRand()) |
---|
| 152 | { |
---|
| 153 | d2 = G4Neutron::NeutronDefinition(); |
---|
| 154 | } |
---|
| 155 | else |
---|
| 156 | { |
---|
| 157 | d1 = G4Neutron::NeutronDefinition(); |
---|
| 158 | } |
---|
| 159 | } |
---|
| 160 | else if(d2->GetPDGCharge()>.5) |
---|
| 161 | { |
---|
| 162 | d2 = G4Neutron::NeutronDefinition(); |
---|
| 163 | } |
---|
| 164 | } |
---|
| 165 | else if(dp->GetPDGCharge()>.5) |
---|
| 166 | { |
---|
| 167 | if(d1->GetPDGCharge()<.5) |
---|
| 168 | { |
---|
| 169 | if(d2->GetPDGCharge()<.5 && 0.5>G4UniformRand()) |
---|
| 170 | { |
---|
| 171 | d2 = G4Proton::ProtonDefinition(); |
---|
| 172 | } |
---|
| 173 | else |
---|
| 174 | { |
---|
| 175 | d1 = G4Proton::ProtonDefinition(); |
---|
| 176 | } |
---|
| 177 | } |
---|
| 178 | else if(d2->GetPDGCharge()<.5) |
---|
| 179 | { |
---|
| 180 | d2 = G4Proton::ProtonDefinition(); |
---|
| 181 | } |
---|
| 182 | } |
---|
| 183 | |
---|
| 184 | // calculate the momenta. |
---|
| 185 | G4double M = (thePro+theT1+theT2).mag(); |
---|
| 186 | G4double m1 = d1->GetPDGMass(); |
---|
| 187 | G4double m2 = d2->GetPDGMass(); |
---|
| 188 | G4double m = std::sqrt(M*M-m1*m1-m2*m2); |
---|
| 189 | G4double p = std::sqrt((m*m*m*m - 4.*m1*m1 * m2*m2)/(4.*(M*M))); |
---|
| 190 | G4double costh = 2.*G4UniformRand()-1.; |
---|
| 191 | G4double phi = 2.*pi*G4UniformRand(); |
---|
| 192 | G4ThreeVector pFinal(p*std::sin(std::acos(costh))*std::cos(phi), p*std::sin(std::acos(costh))*std::sin(phi), p*costh); |
---|
| 193 | |
---|
| 194 | // G4cout << "testing p "<<p-pFinal.mag()<<G4endl; |
---|
| 195 | // construct the final state particles lorentz momentum. |
---|
| 196 | G4double eFinal1 = std::sqrt(m1*m1+pFinal.mag2()); |
---|
| 197 | G4LorentzVector final1(pFinal, eFinal1); |
---|
| 198 | G4double eFinal2 = std::sqrt(m2*m2+pFinal.mag2()); |
---|
| 199 | G4LorentzVector final2(-1.*pFinal, eFinal2); |
---|
| 200 | |
---|
| 201 | // rotate back. |
---|
| 202 | final1 = toLab * final1; |
---|
| 203 | final2 = toLab * final2; |
---|
| 204 | |
---|
| 205 | // boost back to LAB |
---|
| 206 | final1 = fromSPS * final1; |
---|
| 207 | final2 = fromSPS * final2; |
---|
| 208 | |
---|
| 209 | // make the final state |
---|
| 210 | G4KineticTrack * f1 = |
---|
| 211 | new G4KineticTrack(d1, 0., targets[0]->GetPosition(), final1); |
---|
| 212 | G4KineticTrack * f2 = |
---|
| 213 | new G4KineticTrack(d2, 0., targets[1]->GetPosition(), final2); |
---|
| 214 | G4KineticTrackVector * result = new G4KineticTrackVector; |
---|
| 215 | result->push_back(f1); |
---|
| 216 | result->push_back(f2); |
---|
| 217 | |
---|
| 218 | for(size_t hpw=0; hpw<result->size(); hpw++) |
---|
| 219 | { |
---|
| 220 | energyBalance-=result->operator[](hpw)->Get4Momentum().t(); |
---|
| 221 | chargeBalance-=G4lrint(result->operator[](hpw)->GetDefinition()->GetPDGCharge()); |
---|
| 222 | baryonBalance-=result->operator[](hpw)->GetDefinition()->GetBaryonNumber(); |
---|
| 223 | } |
---|
| 224 | if(getenv("AbsorptionEnergyBalanceCheck")) |
---|
| 225 | std::cout << "DEBUGGING energy balance B: " |
---|
| 226 | <<energyBalance<<" " |
---|
| 227 | <<chargeBalance<<" " |
---|
| 228 | <<baryonBalance<<" " |
---|
| 229 | <<G4endl; |
---|
| 230 | return result; |
---|
| 231 | } |
---|
| 232 | |
---|
| 233 | |
---|
| 234 | G4double G4MesonAbsorption:: |
---|
| 235 | GetTimeToAbsorption(const G4KineticTrack& trk1, const G4KineticTrack& trk2) |
---|
| 236 | { |
---|
| 237 | if(trk1.GetDefinition()!=G4PionPlus::PionPlusDefinition() && |
---|
| 238 | trk1.GetDefinition()!=G4PionMinus::PionMinusDefinition() && |
---|
| 239 | trk2.GetDefinition()!=G4PionPlus::PionPlusDefinition() && |
---|
| 240 | trk2.GetDefinition()!=G4PionMinus::PionMinusDefinition()) |
---|
| 241 | { |
---|
| 242 | return DBL_MAX; |
---|
| 243 | } |
---|
| 244 | G4double time = DBL_MAX; |
---|
| 245 | G4double sqrtS = (trk1.Get4Momentum() + trk2.Get4Momentum()).mag(); |
---|
| 246 | |
---|
| 247 | // Check whether there is enough energy for elastic scattering |
---|
| 248 | // (to put the particles on to mass shell |
---|
| 249 | |
---|
| 250 | if (trk1.GetActualMass() + trk2.GetActualMass() < sqrtS) |
---|
| 251 | { |
---|
| 252 | G4LorentzVector mom1 = trk1.GetTrackingMomentum(); |
---|
| 253 | G4ThreeVector position = trk1.GetPosition() - trk2.GetPosition(); |
---|
| 254 | if ( mom1.mag2() < -1.*eV ) |
---|
| 255 | { |
---|
| 256 | G4cout << "G4MesonAbsorption::GetTimeToInteraction(): negative m2:" << mom1.mag2() << G4endl; |
---|
| 257 | } |
---|
| 258 | G4ThreeVector velocity = mom1.vect()/mom1.e() * c_light; |
---|
| 259 | G4double collisionTime = - (position * velocity) / (velocity * velocity); |
---|
| 260 | |
---|
| 261 | if (collisionTime > 0) |
---|
| 262 | { |
---|
| 263 | G4LorentzVector mom2(0,0,0,trk2.Get4Momentum().mag()); |
---|
| 264 | G4LorentzRotation toCMSFrame((-1)*(mom1 + mom2).boostVector()); |
---|
| 265 | mom1 = toCMSFrame * mom1; |
---|
| 266 | mom2 = toCMSFrame * mom2; |
---|
| 267 | |
---|
| 268 | G4LorentzVector coordinate1(trk1.GetPosition(), 100.); |
---|
| 269 | G4LorentzVector coordinate2(trk2.GetPosition(), 100.); |
---|
| 270 | G4ThreeVector pos = ((toCMSFrame * coordinate1).vect() - |
---|
| 271 | (toCMSFrame * coordinate2).vect()); |
---|
| 272 | G4ThreeVector mom = mom1.vect() - mom2.vect(); |
---|
| 273 | |
---|
| 274 | G4double distance = pos * pos - (pos*mom) * (pos*mom) / (mom*mom); |
---|
| 275 | |
---|
| 276 | // global optimization |
---|
| 277 | static const G4double maxCrossSection = 500*millibarn; |
---|
| 278 | if(pi*distance>maxCrossSection) return time; |
---|
| 279 | // charged particles special optimization |
---|
| 280 | static const G4double maxChargedCrossSection = 200*millibarn; |
---|
| 281 | if(std::abs(trk1.GetDefinition()->GetPDGCharge())>0.1 && |
---|
| 282 | std::abs(trk2.GetDefinition()->GetPDGCharge())>0.1 && |
---|
| 283 | pi*distance>maxChargedCrossSection) return time; |
---|
| 284 | // neutrons special optimization |
---|
| 285 | if(( trk1.GetDefinition() == G4Neutron::Neutron() || |
---|
| 286 | trk1.GetDefinition() == G4Neutron::Neutron() ) && |
---|
| 287 | sqrtS>1.91*GeV && pi*distance>maxChargedCrossSection) return time; |
---|
| 288 | |
---|
| 289 | G4double totalCrossSection = AbsorptionCrossSection(trk1,trk2); |
---|
| 290 | if ( totalCrossSection > 0 ) |
---|
| 291 | { |
---|
| 292 | if (distance <= totalCrossSection / pi) |
---|
| 293 | { |
---|
| 294 | time = collisionTime; |
---|
| 295 | } |
---|
| 296 | } |
---|
| 297 | } |
---|
| 298 | } |
---|
| 299 | return time; |
---|
| 300 | } |
---|
| 301 | |
---|
| 302 | G4double G4MesonAbsorption:: |
---|
| 303 | AbsorptionCrossSection(const G4KineticTrack & aT, const G4KineticTrack & bT) |
---|
| 304 | { |
---|
| 305 | G4double t = 0; |
---|
| 306 | if(aT.GetDefinition()==G4PionPlus::PionPlusDefinition() || |
---|
| 307 | aT.GetDefinition()==G4PionMinus::PionMinusDefinition() ) |
---|
| 308 | { |
---|
| 309 | t = aT.Get4Momentum().t()-aT.Get4Momentum().mag()/MeV; |
---|
| 310 | } |
---|
| 311 | else if(bT.GetDefinition()==G4PionPlus::PionPlusDefinition() || |
---|
| 312 | bT.GetDefinition()!=G4PionMinus::PionMinusDefinition()) |
---|
| 313 | { |
---|
| 314 | t = bT.Get4Momentum().t()-bT.Get4Momentum().mag()/MeV; |
---|
| 315 | } |
---|
| 316 | |
---|
| 317 | static G4double it [26] = |
---|
| 318 | {0,4,50,5.5,75,8,95,10,120,11.5,140,12,160,11.5,180,10,190,8,210,6,235,4,260,3,300,2}; |
---|
| 319 | |
---|
| 320 | if(t>it[24]) |
---|
| 321 | { |
---|
| 322 | theCross = 0; |
---|
| 323 | } |
---|
| 324 | else |
---|
| 325 | { |
---|
| 326 | G4int count = 0; |
---|
| 327 | while(t>it[count])count+=2; |
---|
| 328 | G4double x1 = it[count-2]; |
---|
| 329 | G4double x2 = it[count]; |
---|
| 330 | G4double y1 = it[count-1]; |
---|
| 331 | G4double y2 = it[count+1]; |
---|
| 332 | theCross = y1+(y2-y1)/(x2-x1)*(t-x1); |
---|
| 333 | // G4cout << "Printing the absorption crosssection " |
---|
| 334 | // <<x1<< " "<<x2<<" "<<t<<" "<<y1<<" "<<y2<<" "<<0.5*theCross<<G4endl; |
---|
| 335 | } |
---|
| 336 | return .5*theCross*millibarn; |
---|
| 337 | } |
---|