00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027 #include "G4MesonAbsorption.hh"
00028 #include "G4PhysicalConstants.hh"
00029 #include "G4SystemOfUnits.hh"
00030 #include "G4LorentzRotation.hh"
00031 #include "G4LorentzVector.hh"
00032 #include "Randomize.hh"
00033 #include "G4KineticTrackVector.hh"
00034 #include "G4CollisionInitialState.hh"
00035 #include <cmath>
00036 #include "G4PionPlus.hh"
00037 #include "G4PionMinus.hh"
00038 #include "G4ParticleDefinition.hh"
00039 #include "G4HadTmpUtil.hh"
00040
00041
00042
00043 const std::vector<G4CollisionInitialState *> & G4MesonAbsorption::
00044 GetCollisions(G4KineticTrack * aProjectile,
00045 std::vector<G4KineticTrack *> & someCandidates,
00046 G4double aCurrentTime)
00047 {
00048 theCollisions.clear();
00049 if(someCandidates.size() >1)
00050 {
00051 std::vector<G4KineticTrack *>::iterator j=someCandidates.begin();
00052 for(; j != someCandidates.end(); ++j)
00053 {
00054 G4double collisionTime = GetTimeToAbsorption(*aProjectile, **j);
00055 if(collisionTime == DBL_MAX)
00056 {
00057 continue;
00058 }
00059 G4KineticTrackVector aTarget;
00060 aTarget.push_back(*j);
00061 FindAndFillCluster(aTarget, aProjectile, someCandidates);
00062 if(aTarget.size()>=2)
00063 {
00064 theCollisions.push_back(
00065 new G4CollisionInitialState(collisionTime+aCurrentTime, aProjectile, aTarget, this) );
00066 }
00067 }
00068 }
00069 return theCollisions;
00070 }
00071
00072
00073 void G4MesonAbsorption::
00074 FindAndFillCluster(G4KineticTrackVector & result,
00075 G4KineticTrack * aProjectile, std::vector<G4KineticTrack *> & someCandidates)
00076 {
00077 std::vector<G4KineticTrack *>::iterator j=someCandidates.begin();
00078 G4KineticTrack * aTarget = result[0];
00079 G4int chargeSum = G4lrint(aTarget->GetDefinition()->GetPDGCharge());
00080 chargeSum+=G4lrint(aProjectile->GetDefinition()->GetPDGCharge());
00081 G4ThreeVector firstBase = aTarget->GetPosition();
00082 G4double min = DBL_MAX;
00083 G4KineticTrack * partner = 0;
00084 for(; j != someCandidates.end(); ++j)
00085 {
00086 if(*j == aTarget) continue;
00087 G4int cCharge = G4lrint((*j)->GetDefinition()->GetPDGCharge());
00088 if (chargeSum+cCharge > 2) continue;
00089 if (chargeSum+cCharge < 0) continue;
00090
00091 G4ThreeVector secodeBase = (*j)->GetPosition();
00092 if((firstBase+secodeBase).mag()<min)
00093 {
00094 min=(firstBase+secodeBase).mag();
00095 partner = *j;
00096 }
00097 }
00098 if(partner) result.push_back(partner);
00099 else result.clear();
00100 }
00101
00102
00103 G4KineticTrackVector * G4MesonAbsorption::
00104 GetFinalState(G4KineticTrack * projectile,
00105 std::vector<G4KineticTrack *> & targets)
00106 {
00107
00108
00109
00110 G4LorentzVector thePro = projectile->Get4Momentum();
00111 G4LorentzVector theT1 = targets[0]->Get4Momentum();
00112 G4LorentzVector theT2 = targets[1]->Get4Momentum();
00113 G4LorentzVector incoming = thePro + theT1 + theT2;
00114 G4double energyBalance = incoming.t();
00115 G4int chargeBalance = G4lrint(projectile->GetDefinition()->GetPDGCharge()
00116 + targets[0]->GetDefinition()->GetPDGCharge()
00117 + targets[1]->GetDefinition()->GetPDGCharge());
00118
00119 G4int baryonBalance = projectile->GetDefinition()->GetBaryonNumber()
00120 + targets[0]->GetDefinition()->GetBaryonNumber()
00121 + targets[1]->GetDefinition()->GetBaryonNumber();
00122
00123
00124
00125 G4LorentzRotation toSPS((-1)*(thePro + theT1 + theT2).boostVector());
00126 theT1 = toSPS * theT1;
00127 theT2 = toSPS * theT2;
00128 thePro = toSPS * thePro;
00129 G4LorentzRotation fromSPS(toSPS.inverse());
00130
00131
00132 G4LorentzRotation toZ;
00133 G4LorentzVector Ptmp=projectile->Get4Momentum();
00134 toZ.rotateZ(-1*Ptmp.phi());
00135 toZ.rotateY(-1*Ptmp.theta());
00136 theT1 = toZ * theT1;
00137 theT2 = toZ * theT2;
00138 thePro = toZ * thePro;
00139 G4LorentzRotation toLab(toZ.inverse());
00140
00141
00142 G4ParticleDefinition * d1 = targets[0]->GetDefinition();
00143 G4ParticleDefinition * d2 = targets[1]->GetDefinition();
00144 if(0.5>G4UniformRand())
00145 {
00146 G4ParticleDefinition * temp;
00147 temp=d1;d1=d2;d2=temp;
00148 }
00149 G4ParticleDefinition * dp = projectile->GetDefinition();
00150 if(dp->GetPDGCharge()<-.5)
00151 {
00152 if(d1->GetPDGCharge()>.5)
00153 {
00154 if(d2->GetPDGCharge()>.5 && 0.5>G4UniformRand())
00155 {
00156 d2 = G4Neutron::NeutronDefinition();
00157 }
00158 else
00159 {
00160 d1 = G4Neutron::NeutronDefinition();
00161 }
00162 }
00163 else if(d2->GetPDGCharge()>.5)
00164 {
00165 d2 = G4Neutron::NeutronDefinition();
00166 }
00167 }
00168 else if(dp->GetPDGCharge()>.5)
00169 {
00170 if(d1->GetPDGCharge()<.5)
00171 {
00172 if(d2->GetPDGCharge()<.5 && 0.5>G4UniformRand())
00173 {
00174 d2 = G4Proton::ProtonDefinition();
00175 }
00176 else
00177 {
00178 d1 = G4Proton::ProtonDefinition();
00179 }
00180 }
00181 else if(d2->GetPDGCharge()<.5)
00182 {
00183 d2 = G4Proton::ProtonDefinition();
00184 }
00185 }
00186
00187
00188 G4double M_sq = (thePro+theT1+theT2).mag2();
00189 G4double m1_sq = sqr(d1->GetPDGMass());
00190 G4double m2_sq = sqr(d2->GetPDGMass());
00191 G4double m_sq = M_sq-m1_sq-m2_sq;
00192 G4double p = std::sqrt((m_sq*m_sq - 4.*m1_sq * m2_sq)/(4.*M_sq));
00193 G4double costh = 2.*G4UniformRand()-1.;
00194 G4double phi = 2.*pi*G4UniformRand();
00195 G4ThreeVector pFinal(p*std::sin(std::acos(costh))*std::cos(phi), p*std::sin(std::acos(costh))*std::sin(phi), p*costh);
00196
00197
00198
00199 G4double eFinal1 = std::sqrt(m1_sq+pFinal.mag2());
00200 G4LorentzVector final1(pFinal, eFinal1);
00201 G4double eFinal2 = std::sqrt(m2_sq+pFinal.mag2());
00202 G4LorentzVector final2(-1.*pFinal, eFinal2);
00203
00204
00205 final1 = toLab * final1;
00206 final2 = toLab * final2;
00207
00208
00209 final1 = fromSPS * final1;
00210 final2 = fromSPS * final2;
00211
00212
00213 G4KineticTrack * f1 =
00214 new G4KineticTrack(d1, 0., targets[0]->GetPosition(), final1);
00215 G4KineticTrack * f2 =
00216 new G4KineticTrack(d2, 0., targets[1]->GetPosition(), final2);
00217 G4KineticTrackVector * result = new G4KineticTrackVector;
00218 result->push_back(f1);
00219 result->push_back(f2);
00220
00221 for(size_t hpw=0; hpw<result->size(); hpw++)
00222 {
00223 energyBalance-=result->operator[](hpw)->Get4Momentum().t();
00224 chargeBalance-=G4lrint(result->operator[](hpw)->GetDefinition()->GetPDGCharge());
00225 baryonBalance-=result->operator[](hpw)->GetDefinition()->GetBaryonNumber();
00226 }
00227 if(getenv("AbsorptionEnergyBalanceCheck"))
00228 std::cout << "DEBUGGING energy balance B: "
00229 <<energyBalance<<" "
00230 <<chargeBalance<<" "
00231 <<baryonBalance<<" "
00232 <<G4endl;
00233 return result;
00234 }
00235
00236
00237 G4double G4MesonAbsorption::
00238 GetTimeToAbsorption(const G4KineticTrack& trk1, const G4KineticTrack& trk2)
00239 {
00240 if(trk1.GetDefinition()!=G4PionPlus::PionPlusDefinition() &&
00241 trk1.GetDefinition()!=G4PionMinus::PionMinusDefinition() &&
00242 trk2.GetDefinition()!=G4PionPlus::PionPlusDefinition() &&
00243 trk2.GetDefinition()!=G4PionMinus::PionMinusDefinition())
00244 {
00245 return DBL_MAX;
00246 }
00247 G4double time = DBL_MAX;
00248 G4double sqrtS = (trk1.Get4Momentum() + trk2.Get4Momentum()).mag();
00249
00250
00251
00252
00253 if (trk1.GetActualMass() + trk2.GetActualMass() < sqrtS)
00254 {
00255 G4LorentzVector mom1 = trk1.GetTrackingMomentum();
00256 G4ThreeVector position = trk1.GetPosition() - trk2.GetPosition();
00257 if ( mom1.mag2() < -1.*eV )
00258 {
00259 G4cout << "G4MesonAbsorption::GetTimeToInteraction(): negative m2:" << mom1.mag2() << G4endl;
00260 }
00261 G4ThreeVector velocity = mom1.vect()/mom1.e() * c_light;
00262 G4double collisionTime = - (position * velocity) / (velocity * velocity);
00263
00264 if (collisionTime > 0)
00265 {
00266 G4LorentzVector mom2(0,0,0,trk2.Get4Momentum().mag());
00267 G4LorentzRotation toCMSFrame((-1)*(mom1 + mom2).boostVector());
00268 mom1 = toCMSFrame * mom1;
00269 mom2 = toCMSFrame * mom2;
00270
00271 G4LorentzVector coordinate1(trk1.GetPosition(), 100.);
00272 G4LorentzVector coordinate2(trk2.GetPosition(), 100.);
00273 G4ThreeVector pos = ((toCMSFrame * coordinate1).vect() -
00274 (toCMSFrame * coordinate2).vect());
00275 G4ThreeVector mom = mom1.vect() - mom2.vect();
00276
00277 G4double distance = pos * pos - (pos*mom) * (pos*mom) / (mom*mom);
00278
00279
00280 static const G4double maxCrossSection = 500*millibarn;
00281 if(pi*distance>maxCrossSection) return time;
00282
00283 static const G4double maxChargedCrossSection = 200*millibarn;
00284 if(std::abs(trk1.GetDefinition()->GetPDGCharge())>0.1 &&
00285 std::abs(trk2.GetDefinition()->GetPDGCharge())>0.1 &&
00286 pi*distance>maxChargedCrossSection) return time;
00287
00288 if(( trk1.GetDefinition() == G4Neutron::Neutron() ||
00289 trk2.GetDefinition() == G4Neutron::Neutron() ) &&
00290 sqrtS>1.91*GeV && pi*distance>maxChargedCrossSection) return time;
00291
00292 G4double totalCrossSection = AbsorptionCrossSection(trk1,trk2);
00293 if ( totalCrossSection > 0 )
00294 {
00295 if (distance <= totalCrossSection / pi)
00296 {
00297 time = collisionTime;
00298 }
00299 }
00300 }
00301 }
00302 return time;
00303 }
00304
00305 G4double G4MesonAbsorption::
00306 AbsorptionCrossSection(const G4KineticTrack & aT, const G4KineticTrack & bT)
00307 {
00308 G4double t = 0;
00309 if(aT.GetDefinition()==G4PionPlus::PionPlusDefinition() ||
00310 aT.GetDefinition()==G4PionMinus::PionMinusDefinition() )
00311 {
00312 t = aT.Get4Momentum().t()-aT.Get4Momentum().mag()/MeV;
00313 }
00314 else if(bT.GetDefinition()==G4PionPlus::PionPlusDefinition() ||
00315 bT.GetDefinition()!=G4PionMinus::PionMinusDefinition())
00316 {
00317 t = bT.Get4Momentum().t()-bT.Get4Momentum().mag()/MeV;
00318 }
00319
00320 static G4double it [26] =
00321 {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};
00322
00323 G4double aCross(0);
00324 if(t<=it[24])
00325 {
00326 G4int count = 0;
00327 while(t>it[count])count+=2;
00328 G4double x1 = it[count-2];
00329 G4double x2 = it[count];
00330 G4double y1 = it[count-1];
00331 G4double y2 = it[count+1];
00332 aCross = y1+(y2-y1)/(x2-x1)*(t-x1);
00333
00334
00335 }
00336 return .5*aCross*millibarn;
00337 }