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 #include <utility>
00027
00028 #include "G4QGSParticipants.hh"
00029 #include "globals.hh"
00030 #include "G4SystemOfUnits.hh"
00031 #include "G4LorentzVector.hh"
00032
00033
00034
00035
00036
00037 G4int G4QGSParticipants_NPart = 0;
00038
00039 G4QGSParticipants::G4QGSParticipants() : theDiffExcitaton(),
00040 ModelMode(SOFT),
00041
00042 nCutMax(7),ThresholdParameter(0.000*GeV),
00043 QGSMThreshold(3*GeV),theNucleonRadius(1.5*fermi)
00044
00045 {
00046 }
00047
00048 G4QGSParticipants::G4QGSParticipants(const G4QGSParticipants &right)
00049 : G4VParticipants(),ModelMode(right.ModelMode), nCutMax(right.nCutMax),
00050 ThresholdParameter(right.ThresholdParameter), QGSMThreshold(right.QGSMThreshold),
00051 theNucleonRadius(right.theNucleonRadius)
00052 {
00053 }
00054
00055
00056 G4QGSParticipants::~G4QGSParticipants()
00057 {
00058 }
00059
00060 void G4QGSParticipants::BuildInteractions(const G4ReactionProduct &thePrimary)
00061 {
00062
00063
00064 G4VSplitableHadron* aProjectile = SelectInteractions(thePrimary);
00065
00066
00067 SplitHadrons();
00068
00069
00070 PerformSoftCollisions();
00071
00072
00073 PerformDiffractiveCollisions();
00074
00075
00076 std::for_each(theInteractions.begin(), theInteractions.end(), DeleteInteractionContent());
00077 theInteractions.clear();
00078 std::for_each(theTargets.begin(), theTargets.end(), DeleteSplitableHadron());
00079 theTargets.clear();
00080 delete aProjectile;
00081 }
00082
00083 G4VSplitableHadron* G4QGSParticipants::SelectInteractions(const G4ReactionProduct &thePrimary)
00084 {
00085 G4VSplitableHadron* aProjectile = new G4QGSMSplitableHadron(thePrimary, TRUE);
00086 G4PomeronCrossSection theProbability(thePrimary.GetDefinition());
00087 G4double outerRadius = theNucleus->GetOuterRadius();
00088
00089
00090 theNucleus->StartLoop();
00091 G4Nucleon * pNucleon = theNucleus->GetNextNucleon();
00092 G4LorentzVector aPrimaryMomentum(thePrimary.GetMomentum(), thePrimary.GetTotalEnergy());
00093
00094
00095 G4double s_nucleus = (aPrimaryMomentum + pNucleon->Get4Momentum()).mag2();
00096 G4double ThresholdMass = thePrimary.GetMass() + pNucleon->GetDefinition()->GetPDGMass();
00097 ModelMode = SOFT;
00098 if (sqr(ThresholdMass + ThresholdParameter) > s_nucleus)
00099 {
00100 ModelMode = DIFFRACTIVE;
00101
00102 }
00103 if (sqr(ThresholdMass + QGSMThreshold) > s_nucleus)
00104 {
00105 ModelMode = DIFFRACTIVE;
00106 }
00107
00108
00109 std::for_each(theInteractions.begin(), theInteractions.end(), DeleteInteractionContent());
00110 theInteractions.clear();
00111 G4int totalCuts = 0;
00112
00113 #ifdef debug_QGS
00114 G4double eK = thePrimary.GetKineticEnergy()/GeV;
00115 #endif
00116 #ifdef debug_G4QGSParticipants
00117 G4double impactUsed = 0;
00118 G4LorentzVector intNuclMom;
00119 #endif
00120
00121 while(theInteractions.size() == 0)
00122 {
00123
00124 std::pair<G4double, G4double> theImpactParameter;
00125 theImpactParameter = theNucleus->ChooseImpactXandY(outerRadius+theNucleonRadius);
00126 G4double impactX = theImpactParameter.first;
00127 G4double impactY = theImpactParameter.second;
00128
00129
00130 theNucleus->StartLoop();
00131 G4int nucleonCount = 0;
00132 G4QGSParticipants_NPart = 0;
00133 #ifdef debug_G4QGSParticipants
00134 intNuclMom=aPrimaryMomentum;
00135 #endif
00136 while( (pNucleon = theNucleus->GetNextNucleon()) )
00137 {
00138 if(totalCuts>1.5*thePrimary.GetKineticEnergy()/GeV)
00139 {
00140 break;
00141 }
00142 nucleonCount++;
00143
00144 G4LorentzVector nucleonMomentum=pNucleon->Get4Momentum();
00145 nucleonMomentum.setE(nucleonMomentum.e()-pNucleon->GetBindingEnergy());
00146 G4double s_nucleon = (aPrimaryMomentum + nucleonMomentum).mag2();
00147 G4double Distance2 = sqr(impactX - pNucleon->GetPosition().x()) +
00148 sqr(impactY - pNucleon->GetPosition().y());
00149 G4double Probability = theProbability.GetInelasticProbability(s_nucleon, Distance2);
00150
00151 G4double rndNumber = G4UniformRand();
00152
00153 if (Probability > rndNumber)
00154 {
00155 #ifdef debug_G4QGSParticipants
00156 G4cout << "DEBUG p="<< Probability<<" r="<<rndNumber<<" d="<<std::sqrt(Distance2)<<G4endl;
00157 G4cout << " qgspart+ " << aPrimaryMomentum << " # " << aPrimaryMomentum.mag()
00158 << pNucleon->Get4Momentum() << " # " << (pNucleon->Get4Momentum()).mag()<< G4endl;
00159 intNuclMom += nucleonMomentum;
00160 #endif
00161 pNucleon->SetMomentum(nucleonMomentum);
00162 G4QGSMSplitableHadron* aTarget = new G4QGSMSplitableHadron(*pNucleon);
00163 G4QGSParticipants_NPart ++;
00164 theTargets.push_back(aTarget);
00165 pNucleon->Hit(aTarget);
00166 if ((theProbability.GetDiffractiveProbability(s_nucleon, Distance2)/Probability > G4UniformRand()
00167 &&(ModelMode==SOFT)) || (ModelMode==DIFFRACTIVE ))
00168 {
00169
00170 if(IsSingleDiffractive())
00171 {
00172 theSingleDiffExcitation.ExciteParticipants(aProjectile, aTarget);
00173 }
00174 else
00175 {
00176 theDiffExcitaton.ExciteParticipants(aProjectile, aTarget);
00177 }
00178 G4InteractionContent * aInteraction = new G4InteractionContent(aProjectile);
00179 aInteraction->SetTarget(aTarget);
00180 theInteractions.push_back(aInteraction);
00181 aInteraction->SetNumberOfDiffractiveCollisions(1);
00182 totalCuts += 1;
00183 }
00184 else
00185 {
00186
00187
00188 G4int nCut;
00189 G4double * running = new G4double[nCutMax];
00190 running[0] = 0;
00191 for(nCut = 0; nCut < nCutMax; nCut++)
00192 {
00193 running[nCut] = theProbability.GetCutPomeronProbability(s_nucleon, Distance2, nCut + 1);
00194 if(nCut!=0) running[nCut] += running[nCut-1];
00195 }
00196 G4double random = running[nCutMax-1]*G4UniformRand();
00197 for(nCut = 0; nCut < nCutMax; nCut++)
00198 {
00199 if(running[nCut] > random) break;
00200 }
00201 delete [] running;
00202 nCut = 0;
00203 aTarget->IncrementCollisionCount(nCut+1);
00204 aProjectile->IncrementCollisionCount(nCut+1);
00205 G4InteractionContent * aInteraction = new G4InteractionContent(aProjectile);
00206 aInteraction->SetTarget(aTarget);
00207 aInteraction->SetNumberOfSoftCollisions(nCut+1);
00208 theInteractions.push_back(aInteraction);
00209 totalCuts += nCut+1;
00210 #ifdef debug_G4QGSParticipants
00211 impactUsed=Distance2;
00212 #endif
00213 }
00214 }
00215 }
00216 #ifdef debug_G4QGSParticipants
00217 G4cout << G4endl<<"NUCLEONCOUNT "<<nucleonCount<<G4endl;
00218 G4cout << " Interact 4-Vect " << intNuclMom << G4endl;
00219 #endif
00220 }
00221 #ifdef debug_G4QGSParticipants
00222 G4cout << G4endl<<"CUTDEBUG "<< totalCuts <<G4endl;
00223 G4cout << "Impact Parameter used = "<<impactUsed<<G4endl;
00224 #endif
00225 return aProjectile;
00226 }
00227
00228 void G4QGSParticipants::PerformDiffractiveCollisions()
00229 {
00230
00231 unsigned int i;
00232 for(i = 0; i < theInteractions.size(); i++)
00233 {
00234 G4InteractionContent* anIniteraction = theInteractions[i];
00235 G4VSplitableHadron* aProjectile = anIniteraction->GetProjectile();
00236 G4Parton* aParton = aProjectile->GetNextParton();
00237 G4PartonPair * aPartonPair;
00238
00239 if (aParton)
00240 {
00241 aPartonPair = new G4PartonPair(aParton, aProjectile->GetNextAntiParton(),
00242 G4PartonPair::DIFFRACTIVE,
00243 G4PartonPair::PROJECTILE);
00244 #ifdef debug_G4QGSPart_PDiffColl
00245 G4cout << "DiffPair Pro " << aPartonPair->GetParton1()->GetPDGcode() << " "
00246 << aPartonPair->GetParton1()->Get4Momentum() << " "
00247 << aPartonPair->GetParton1()->GetX() << " " << G4endl;
00248 G4cout << " " << aPartonPair->GetParton2()->GetPDGcode() << " "
00249 << aPartonPair->GetParton2()->Get4Momentum() << " "
00250 << aPartonPair->GetParton2()->GetX() << " " << G4endl;
00251 #endif
00252 thePartonPairs.push_back(aPartonPair);
00253 }
00254
00255 G4VSplitableHadron* aTarget = anIniteraction->GetTarget();
00256 aParton = aTarget->GetNextParton();
00257 if (aParton)
00258 {
00259 aPartonPair = new G4PartonPair(aParton, aTarget->GetNextAntiParton(),
00260 G4PartonPair::DIFFRACTIVE,
00261 G4PartonPair::TARGET);
00262 #ifdef debug_G4QGSPart_PDiffColl
00263 G4cout << "DiffPair Tgt " << aPartonPair->GetParton1()->GetPDGcode() << " "
00264 << aPartonPair->GetParton1()->Get4Momentum() << " "
00265 << aPartonPair->GetParton1()->GetX() << " " << G4endl;
00266 G4cout << " " << aPartonPair->GetParton2()->GetPDGcode() << " "
00267 << aPartonPair->GetParton2()->Get4Momentum() << " "
00268 << aPartonPair->GetParton2()->GetX() << " " << G4endl;
00269 #endif
00270 thePartonPairs.push_back(aPartonPair);
00271 }
00272 }
00273 }
00274
00275 void G4QGSParticipants::PerformSoftCollisions()
00276 {
00277 std::vector<G4InteractionContent*>::iterator i;
00278 G4LorentzVector str4Mom;
00279 i = theInteractions.begin();
00280 while ( i != theInteractions.end() )
00281 {
00282 G4InteractionContent* anIniteraction = *i;
00283 G4PartonPair * aPair = NULL;
00284 if (anIniteraction->GetNumberOfSoftCollisions())
00285 {
00286 G4VSplitableHadron* pProjectile = anIniteraction->GetProjectile();
00287 G4VSplitableHadron* pTarget = anIniteraction->GetTarget();
00288 for (G4int j = 0; j < anIniteraction->GetNumberOfSoftCollisions(); j++)
00289 {
00290 aPair = new G4PartonPair(pTarget->GetNextParton(), pProjectile->GetNextAntiParton(),
00291 G4PartonPair::SOFT, G4PartonPair::TARGET);
00292 #ifdef debug_G4QGSPart_PSoftColl
00293 G4cout << "SoftPair " << aPair->GetParton1()->GetPDGcode() << " "
00294 << aPair->GetParton1()->Get4Momentum() << " "
00295 << aPair->GetParton1()->GetX() << " " << G4endl;
00296 G4cout << " " << aPair->GetParton2()->GetPDGcode() << " "
00297 << aPair->GetParton2()->Get4Momentum() << " "
00298 << aPair->GetParton2()->GetX() << " " << G4endl;
00299 #endif
00300 #ifdef debug_G4QGSParticipants
00301 str4Mom += aPair->GetParton1()->Get4Momentum();
00302 str4Mom += aPair->GetParton2()->Get4Momentum();
00303 #endif
00304 thePartonPairs.push_back(aPair);
00305 aPair = new G4PartonPair(pProjectile->GetNextParton(), pTarget->GetNextAntiParton(),
00306 G4PartonPair::SOFT, G4PartonPair::PROJECTILE);
00307 #ifdef debug_G4QGSPart_PSoftColl
00308 G4cout << "SoftPair " << aPair->GetParton1()->GetPDGcode() << " "
00309 << aPair->GetParton1()->Get4Momentum() << " "
00310 << aPair->GetParton1()->GetX() << " " << G4endl;
00311 G4cout << " " << aPair->GetParton2()->GetPDGcode() << " "
00312 << aPair->GetParton2()->Get4Momentum() << " "
00313 << aPair->GetParton2()->GetX() << " " << G4endl;
00314 #endif
00315 #ifdef debug_G4QGSParticipants
00316 str4Mom += aPair->GetParton1()->Get4Momentum();
00317 str4Mom += aPair->GetParton2()->Get4Momentum();
00318 #endif
00319 thePartonPairs.push_back(aPair);
00320 }
00321 delete *i;
00322 i=theInteractions.erase(i);
00323 } else {
00324 i++;
00325 }
00326 }
00327 #ifdef debug_G4QGSPart_PSoftColl
00328 G4cout << " string 4 mom " << str4Mom << G4endl;
00329 #endif
00330 }