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
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049 #include "G4QIonIonCollision.hh"
00050 #include "G4PhysicalConstants.hh"
00051 #include "G4SystemOfUnits.hh"
00052
00053
00054
00055
00056 G4int G4QIonIonCollision::nCutMax=7;
00057 G4double G4QIonIonCollision::stringTension=1.*GeV/fermi;
00058 G4double G4QIonIonCollision::tubeDensity =1./fermi;
00059
00060 G4double G4QIonIonCollision::widthOfPtSquare=-0.75*GeV*GeV;
00061
00062 G4QIonIonCollision::G4QIonIonCollision(G4QNucleus &pNucleus, const G4QNucleus &tNucleus)
00063 {
00064 static const G4double mProt= G4QPDGCode(2212).GetMass();
00065 static const G4double mPi0= G4QPDGCode(111).GetMass();
00066 theWorld= G4QCHIPSWorld::Get();
00067 theResult = new G4QHadronVector;
00068 G4bool stringsInitted=false;
00069 G4LorentzRotation toZ;
00070 G4LorentzVector proj4M=pNucleus.Get4Momentum();
00071
00072 G4double targM=tNucleus.GetGSMass();
00073 #ifdef edebug
00074 G4cout<<"G4QIn::Constr:*Called*,pA="<<pNucleus<<proj4M<<",tA="<<tNucleus<<targM<<G4endl;
00075 #endif
00076 G4int pZ=pNucleus.GetZ();
00077 G4int pN=pNucleus.GetN();
00078 G4int pA=pZ+pN;
00079 G4int pPDG=90000000+pZ*1000+pN;
00080 G4int tZ=tNucleus.GetZ();
00081 G4int tN=tNucleus.GetN();
00082 G4int tA=tZ+tN;
00083 G4int tPDG=90000000+tZ*1000+tN;
00084 toZ.rotateZ(-proj4M.phi());
00085 toZ.rotateY(-proj4M.theta());
00086 G4LorentzVector zProj4M=toZ*proj4M;
00087 pNucleus.Set4Momentum(zProj4M);
00088 #ifdef edebug
00089 G4int totChg=pZ+tZ;
00090 G4int totBaN=pA+tA;
00091 G4LorentzVector tgLS4M(0.,0.,0.,targM);
00092 G4LorentzVector totLS4M=proj4M+tgLS4M;
00093 G4LorentzVector totZLS4M=zProj4M+tgLS4M;
00094 G4cout<<"-EMC-G4QIonIonCollision::Constr: tLS4M="<<totLS4M<<",tZLS4M="<<totZLS4M<<G4endl;
00095
00096 #endif
00097 G4LorentzRotation toLab(toZ.inverse());
00098 G4QInteractionVector theInteractions;
00099 G4QPartonPairVector thePartonPairs;
00100 G4int ModelMode=SOFT;
00101 theTargNucleus=G4QNucleus(tZ,tN);
00102 theTargNucleus.InitByPDG(tPDG);
00103 theTargNucleus.Init3D();
00104 #ifdef debug
00105 G4cout<<"G4QIonIonCollision::Constr:TargNuclMom="<<theTargNucleus.Get4Momentum()<<G4endl;
00106 #endif
00107 theProjNucleus=G4QNucleus(pZ,pN);
00108 theProjNucleus.InitByPDG(pPDG);
00109 theProjNucleus.Init3D();
00110 #ifdef debug
00111 G4cout<<"G4QIonIonCollision::Constr:ProjNuclMom="<<theProjNucleus.Get4Momentum()<<G4endl;
00112 #endif
00113 #ifdef edebug
00114 G4LorentzVector sumP1=theProjNucleus.GetNucleons4Momentum();
00115 G4LorentzVector sumT1=theTargNucleus.GetNucleons4Momentum();
00116 G4cout<<"-EMC-G4QIonIonCollision::Construct: pNuc4M="<<sumP1<<", tNuc4M="<<sumT1<<G4endl;
00117 #endif
00118 G4ThreeVector theCMVelocity(0.,0.,0.);
00119 G4ThreeVector theALVelocity(0.,0.,0.);
00120
00121 proj4M = pNucleus.Get4Momentum();
00122 G4double pz_per_projectile = proj4M.pz()/pA;
00123 G4double e_per_projectile = proj4M.e()/pA+targM/tA;
00124 G4double vz = pz_per_projectile/e_per_projectile;
00125 #ifdef debug
00126 G4cout<<"G4QIonIonCollision::Construct: (КщеЯ)Projectile4M="<<proj4M<<", Vz="<<vz
00127 <<", pA="<<pA<<", pE="<<e_per_projectile<<G4endl;
00128 #endif
00129 theCMVelocity.setZ(vz);
00130 theProjNucleus.DoLorentzBoost(-theCMVelocity);
00131 theTargNucleus.DoLorentzBoost(-theCMVelocity);
00132 G4LorentzVector prN4M=theProjNucleus.Get4Momentum();
00133 G4double avz=prN4M.pz()/prN4M.e();
00134 theALVelocity.setZ(avz);
00135 #ifdef edebug
00136 G4LorentzVector sumP2=theProjNucleus.GetNucleons4Momentum();
00137 G4LorentzVector sumT2=theTargNucleus.GetNucleons4Momentum();
00138 G4cout<<"-EMC-G4QIonIonCollision::Construct: Boosted, vCM="<<vz<<", vAL="<<avz<<", tN4M="
00139 <<sumT2<<", pN4M="<<sumP2<<G4endl;
00140 #endif
00141 G4int maxCuts = 7;
00142
00143
00144
00145 G4double outerRadius = theProjNucleus.GetOuterRadius()+theTargNucleus.GetOuterRadius();
00146 G4QProbability theProbability(2212);
00147
00148 #ifdef debug
00149 G4cout<<"G4QIonIonCollision::Construct: theIntSize="<<theInteractions.size()<<G4endl;
00150 #endif
00151 G4int attCnt=-1;
00152 G4int maxAtt=27;
00153
00154
00155 G4QHadron* pNucleon;
00156 G4QHadron* tNucleon;
00157 G4QNucleus curProjNucleus;
00158 G4QNucleus curTargNucleus;
00159 while(!theInteractions.size() && ++attCnt < maxAtt)
00160 {
00161
00162
00163 G4int nint=theInteractions.size();
00164 for(G4int ni=0; ni < nint; ++ni)
00165 {
00166 G4QInteraction* curInt = theInteractions[ni];
00167 delete curInt->GetProjectile();
00168 delete curInt->GetTarget();
00169 delete curInt;
00170 }
00171 theInteractions.clear();
00172
00173 if(attCnt)
00174 {
00175 curProjNucleus.DeleteNucleons();
00176 curTargNucleus.DeleteNucleons();
00177 }
00178 curProjNucleus = theProjNucleus;
00179 curTargNucleus = theTargNucleus;
00180
00181 std::pair<G4double, G4double> theImpactParameter;
00182 theImpactParameter = curProjNucleus.ChooseImpactXandY(outerRadius);
00183 G4double impactX = theImpactParameter.first;
00184 G4double impactY = theImpactParameter.second;
00185 #ifdef debug
00186 G4cout<<"G4QIonIonCollision::Con:At#="<<attCnt<<",X="<<impactX<<",Y="<<impactY<<G4endl;
00187 #endif
00188 curProjNucleus.StartLoop();
00189 G4int pnCnt=0;
00190 G4int pnA=curProjNucleus.GetA();
00191 #ifdef debu
00192 G4cout<<"G4QIonIonCollision::Constr: Before the WHILE over pNucleons,pA="<<pnA<<G4endl;
00193 #endif
00194 while((pNucleon=curProjNucleus.GetNextNucleon()) && pnCnt < pnA)
00195 {
00196 ++pnCnt;
00197 G4QInteractionVector curInteractions;
00198 G4LorentzVector pNuc4M=pNucleon->Get4Momentum();
00199 G4ThreeVector pNucR=pNucleon->GetPosition();
00200 G4double pImpX = impactX+pNucR.x();
00201 G4double pImpY = impactY+pNucR.y();
00202 G4double prEn=proj4M.e();
00203 G4int proB=pNucleus.GetBaryonNumber();
00204 if (proB>0) prEn-=pNucleus.GetMass();
00205 else if(proB<0) prEn+=mProt;
00206 #ifdef debug
00207 G4double impactUsed = 0.;
00208 G4cout<<"G4QIonIonCollision::Construct: estimated energy, prEn="<<prEn<<G4endl;
00209 #endif
00210 G4int totalCuts = 0;
00211
00212 G4int tnA=curTargNucleus.GetA();
00213 G4double pImp=std::sqrt(pImpX*pImpX+pImpY*pImpY);
00214 G4double eflen=curTargNucleus.GetThickness(pImp);
00215 maxEn=eflen*stringTension;
00216 maxNuc=static_cast<int>(eflen*tubeDensity+0.5);
00217 #ifdef debug
00218 G4cout<<"G4QIonIonCollision::Con:pE="<<prEn<<" < mE="<<maxEn<<",mN="<<maxNuc<<G4endl;
00219 #endif
00220 if(prEn < maxEn)
00221 {
00222 G4QHadron* aProjectile = new G4QHadron(*pNucleon);
00223 curTargNucleus.StartLoop();
00224 tNucleon=curTargNucleus.GetNextNucleon();
00225 G4QHadron* aTarget = new G4QHadron(*tNucleon);
00226 G4QInteraction* anInteraction = new G4QInteraction(aProjectile);
00227 anInteraction->SetTarget(aTarget);
00228 anInteraction->SetNumberOfDINRCollisions(1);
00229 curInteractions.push_back(anInteraction);
00230 curProjNucleus.DoLorentzBoost(-theALVelocity);
00231 curProjNucleus.SubtractNucleon(pNucleon);
00232 curProjNucleus.DoLorentzBoost(theALVelocity);
00233 curTargNucleus.DoLorentzBoost(theCMVelocity);
00234 curTargNucleus.SubtractNucleon(tNucleon);
00235 curTargNucleus.DoLorentzBoost(-theCMVelocity);
00236 #ifdef debug
00237 G4cout<<"G4QIonIonCollision::Construct: DINR interaction is created"<<G4endl;
00238 #endif
00239 break;
00240 }
00241
00242 curTargNucleus.StartLoop();
00243 G4int tnCnt = 0;
00244 #ifdef debu
00245 G4cout<<"G4QIonIonCollision::Constr: Before WHILE over tNucleons, tA="<<tnA<<G4endl;
00246 #endif
00247 while((tNucleon=curTargNucleus.GetNextNucleon()) && tnCnt<tnA && totalCuts<maxCuts)
00248 {
00249 ++tnCnt;
00250 G4LorentzVector tNuc4M=tNucleon->Get4Momentum();
00251 #ifdef debug
00252 G4cout<<"G4QIonIonCollision::Constr: OuterR="<<outerRadius<<", mC="<<maxCuts
00253 <<", pA="<<curProjNucleus<<", tA="<<curTargNucleus<<G4endl;
00254 #endif
00255
00256 G4double s_value = (tNuc4M + pNuc4M).mag2();
00257 G4double ThresholdMass = pNucleon->GetMass() + tNucleon->GetMass();
00258 #ifdef debug
00259 G4cout<<"G4QInel::Constr: s="<<s_value<<", ThreshM="<<sqr(ThresholdMass)<<G4endl;
00260 #endif
00261 ModelMode = SOFT;
00262 if (s_value < 0.)
00263 {
00264 G4cerr<<"*G4QInelast::Constr:s="<<s_value<<",pN4M="<<pNuc4M<<",tN4M="<<tNuc4M<<G4endl;
00265 G4Exception("G4QIonIonCollision::Construct:","72",FatalException,"LowEn(NegS)");
00266 }
00267 if (s_value < sqr(ThresholdMass))
00268 {
00269 #ifdef debug
00270 G4cout<<"G4QIonIonCollision::Constr: ***OnlyDiffraction***, ThM="<<ThresholdMass
00271 <<">sqrt(s)="<<std::sqrt(s_value)<<" -> only Diffraction is possible"<<G4endl;
00272 #endif
00273 ModelMode = DIFFRACTIVE;
00274 }
00275 G4ThreeVector tNucR=tNucleon->GetPosition();
00276 G4double dImpX = pImpX-tNucR.x();
00277 G4double dImpY = pImpY-tNucR.y();
00278 G4double Distance2=dImpX*dImpX+dImpY*dImpY;
00279 #ifdef sdebug
00280 G4cout<<"G4QIonIonCollision::Construct: s="<<s_value<<", D2="<<Distance2<<G4endl;
00281 #endif
00282
00283 if(s_value<=10000.)
00284 {
00285 G4cout<<"-Warning-G4QIonIonCollision::Construct: s < .01 GeV^2, p4M="
00286 <<pNucleon->Get4Momentum()<<",t4M="<<tNucleon->Get4Momentum()<<G4endl;
00287 continue;
00288 }
00289 G4double Probability = theProbability.GetPomInelProbability(s_value, Distance2);
00290
00291 #ifdef sdebug
00292 G4cout<<"G4QIonIonCollision::Construct: Probubility="<<Probability<<G4endl;
00293 #endif
00294 G4double rndNumber = G4UniformRand();
00295
00296 #ifdef sdebug
00297 G4cout<<"G4QIonIonCollision::Construct: NLOOP prob="<<Probability<<", rndm="
00298 <<rndNumber<<", d="<<std::sqrt(Distance2)<<G4endl;
00299 #endif
00300 if (Probability > rndNumber)
00301 {
00302 G4QHadron* aProjectile = new G4QHadron(*pNucleon);
00303 G4QHadron* aTarget = new G4QHadron(*tNucleon);
00304 #ifdef edebug
00305 G4cout<<"--->EMC-->G4QIonIonCollision::Construct: TargNucleon is filled, 4M/PDG="
00306 <<aTarget->Get4Momentum()<<aTarget->GetPDGCode()<<G4endl;
00307 #endif
00308
00309 curProjNucleus.DoLorentzBoost(-theALVelocity);
00310 curProjNucleus.SubtractNucleon(pNucleon);
00311 curProjNucleus.DoLorentzBoost(theALVelocity);
00312 curTargNucleus.DoLorentzBoost(theCMVelocity);
00313 curTargNucleus.SubtractNucleon(tNucleon);
00314 curTargNucleus.DoLorentzBoost(-theCMVelocity);
00315 if((theProbability.GetPomDiffProbability(s_value,Distance2)/Probability >
00316 G4UniformRand() && ModelMode==SOFT ) || ModelMode==DIFFRACTIVE)
00317 {
00318
00319 if(IsSingleDiffractive()) ExciteSingDiffParticipants(aProjectile, aTarget);
00320 else ExciteDiffParticipants(aProjectile, aTarget);
00321 G4QInteraction* anInteraction = new G4QInteraction(aProjectile);
00322 anInteraction->SetTarget(aTarget);
00323 anInteraction->SetNumberOfDiffractiveCollisions(1);
00324 curInteractions.push_back(anInteraction);
00325
00326 totalCuts++;
00327 #ifdef debug
00328 G4cout<<"G4QIonIonCollision::Constr:NLOOP DiffInteract,tC="<<totalCuts<<G4endl;
00329 #endif
00330 }
00331 else
00332 {
00333
00334
00335 G4int nCut;
00336 G4double* running = new G4double[nCutMax];
00337 for(nCut = 0; nCut < nCutMax; nCut++)
00338 {
00339 running[nCut]= theProbability.GetCutPomProbability(s_value, Distance2, nCut+1);
00340 if(nCut) running[nCut] += running[nCut-1];
00341 }
00342 G4double random = running[nCutMax-1]*G4UniformRand();
00343 for(nCut = 0; nCut < nCutMax; nCut++) if(running[nCut] > random) break;
00344 delete [] running;
00345 #ifdef debug
00346 G4cout<<"G4QIonIonCollision::Construct: NLOOP-Soft Chosen nCut="<<nCut<<G4endl;
00347 #endif
00348
00349
00350
00351 aTarget->IncrementCollisionCount(nCut+1);
00352 aProjectile->IncrementCollisionCount(nCut+1);
00353 G4QInteraction* anInteraction = new G4QInteraction(aProjectile);
00354 anInteraction->SetTarget(aTarget);
00355 anInteraction->SetNumberOfSoftCollisions(nCut+1);
00356 curInteractions.push_back(anInteraction);
00357 totalCuts += nCut+1;
00358 #ifdef debug
00359 G4cout<<"G4QIonIonCollision::Constr:NLOOP SoftInteract,tC="<<totalCuts<<G4endl;
00360 impactUsed=Distance2;
00361 #endif
00362 }
00363 }
00364 }
00365 G4int nCurInt=curInteractions.size();
00366 for(G4int ni=0; ni < nCurInt; ++ni) theInteractions.push_back(curInteractions[ni]);
00367 curInteractions.clear();
00368 }
00369 }
00370 theProjNucleus.DeleteNucleons();
00371 theTargNucleus.DeleteNucleons();
00372 theProjNucleus = curProjNucleus;
00373 theTargNucleus = curTargNucleus;
00374 curProjNucleus.DeleteNucleons();
00375 curTargNucleus.DeleteNucleons();
00376 G4int nInt=theInteractions.size();
00377 #ifdef debug
00378 G4cout<<"G4QIonIonCollision::Constr: #ofInteractions = "<<nInt<<", #OfDINR = "
00379 <<theInteractions[0]->GetNumberOfDINRCollisions()<<G4endl;
00380 #endif
00381 if(!nInt || (nInt==1 && theInteractions[0]->GetNumberOfDINRCollisions()==1))
00382 {
00383 G4QHadron* aTarget;
00384 G4QHadron* aProjectile;
00385 if(nInt)
00386 {
00387 aTarget=theInteractions[0]->GetTarget();
00388 aProjectile=theInteractions[0]->GetProjectile();
00389 delete theInteractions[0];
00390 theInteractions.clear();
00391 }
00392 else
00393 {
00394 theProjNucleus.StartLoop();
00395 pNucleon=theProjNucleus.GetNextNucleon();
00396 aProjectile = new G4QHadron(*pNucleon);
00397 theProjNucleus.DoLorentzBoost(-theALVelocity);
00398 theProjNucleus.SubtractNucleon(pNucleon);
00399 theProjNucleus.DoLorentzBoost(theALVelocity);
00400 theTargNucleus.StartLoop();
00401 tNucleon=theTargNucleus.GetNextNucleon();
00402 aTarget = new G4QHadron(*tNucleon);
00403 theTargNucleus.DoLorentzBoost(theCMVelocity);
00404 theTargNucleus.SubtractNucleon(tNucleon);
00405 theTargNucleus.DoLorentzBoost(-theCMVelocity);
00406 }
00407 G4QContent QQC=aTarget->GetQC()+aProjectile->GetQC();
00408 G4LorentzVector Q4M=aTarget->Get4Momentum()+aProjectile->Get4Momentum();
00409 delete aTarget;
00410 delete aProjectile;
00411
00412 Q4M.boost(theCMVelocity);
00413 Q4M=toLab*Q4M;
00414 G4Quasmon* stringQuasmon = new G4Quasmon(QQC, Q4M);
00415 theQuasmons.push_back(stringQuasmon);
00416 theTargNucleus.DoLorentzBoost(theCMVelocity);
00417 theTargNucleus.DoLorentzRotation(toLab);
00418 return;
00419 }
00420
00421
00422
00423 for(G4int i=0; i<nInt; i++)
00424 {
00425 theInteractions[i]->SplitHadrons();
00426 #ifdef edebug
00427 G4QHadron* projH=theInteractions[i]->GetProjectile();
00428 G4QHadron* targH=theInteractions[i]->GetTarget();
00429 G4LorentzVector pSP(0.,0.,0.,0.);
00430 G4LorentzVector tSP(0.,0.,0.,0.);
00431 std::list<G4QParton*> projCP=projH->GetColor();
00432 std::list<G4QParton*> projAC=projH->GetAntiColor();
00433 std::list<G4QParton*> targCP=targH->GetColor();
00434 std::list<G4QParton*> targAC=targH->GetAntiColor();
00435 std::list<G4QParton*>::iterator picp = projCP.begin();
00436 std::list<G4QParton*>::iterator pecp = projCP.end();
00437 std::list<G4QParton*>::iterator piac = projAC.begin();
00438 std::list<G4QParton*>::iterator peac = projAC.end();
00439 std::list<G4QParton*>::iterator ticp = targCP.begin();
00440 std::list<G4QParton*>::iterator tecp = targCP.end();
00441 std::list<G4QParton*>::iterator tiac = targAC.begin();
00442 std::list<G4QParton*>::iterator teac = targAC.end();
00443 for(; picp!=pecp&& piac!=peac&& ticp!=tecp&& tiac!=teac; ++picp,++piac,++ticp,++tiac)
00444 {
00445 pSP+=(*picp)->Get4Momentum();
00446 pSP+=(*piac)->Get4Momentum();
00447 tSP+=(*ticp)->Get4Momentum();
00448 tSP+=(*tiac)->Get4Momentum();
00449 }
00450 G4cout<<"-EMC-G4QIonIonCollision::Construct: Interaction#"<<i<<",dP4M="
00451 <<projH->Get4Momentum()-pSP<<",dT4M="<<targH->Get4Momentum()-tSP<<G4endl;
00452 #endif
00453 }
00454
00455
00456
00457 G4QInteractionVector::iterator it;
00458 #ifdef debug
00459 G4cout<<"G4QIonIonCollision::Constr: Creation ofSoftCollisionPartonPair STARTS"<<G4endl;
00460 #endif
00461 G4bool rep=true;
00462 while(rep && theInteractions.size())
00463 {
00464 for(it = theInteractions.begin(); it != theInteractions.end(); ++it)
00465 {
00466 G4QInteraction* anIniteraction = *it;
00467 G4QPartonPair* aPair=0;
00468 G4int nSoftCollisions = anIniteraction->GetNumberOfSoftCollisions();
00469 #ifdef debug
00470 G4cout<<"G4QIonIonCollision::Construct: #0f SoftCollisions ="<<nSoftCollisions<<G4endl;
00471 #endif
00472 if (nSoftCollisions)
00473 {
00474 G4QHadron* pProjectile = anIniteraction->GetProjectile();
00475 G4QHadron* pTarget = anIniteraction->GetTarget();
00476 for (G4int j = 0; j < nSoftCollisions; j++)
00477 {
00478 aPair = new G4QPartonPair(pTarget->GetNextParton(),
00479 pProjectile->GetNextAntiParton(),
00480 G4QPartonPair::SOFT, G4QPartonPair::TARGET);
00481 thePartonPairs.push_back(aPair);
00482 aPair = new G4QPartonPair(pProjectile->GetNextParton(),
00483 pTarget->GetNextAntiParton(),
00484 G4QPartonPair::SOFT, G4QPartonPair::PROJECTILE);
00485 thePartonPairs.push_back(aPair);
00486 #ifdef debug
00487 G4cout<<"--->G4QIonIonCollision::Constr: SOFT, 2 parton pairs are filled"<<G4endl;
00488 #endif
00489 }
00490 delete *it;
00491 it=theInteractions.erase(it);
00492 if( it != theInteractions.begin() )
00493 {
00494 it--;
00495 rep=false;
00496 }
00497 else
00498 {
00499 rep=true;
00500 break;
00501 }
00502 }
00503 else rep=false;
00504 }
00505 }
00506 #ifdef debug
00507 G4cout<<"G4QIonIonCollision::Constr: -> Parton pairs for SOFT strings are made"<<G4endl;
00508 #endif
00509
00510
00511
00512 for(unsigned i = 0; i < theInteractions.size(); i++)
00513 {
00514
00515 G4QInteraction* anIniteraction = theInteractions[i];
00516 G4QPartonPair* aPartonPair;
00517 #ifdef debug
00518 G4cout<<"G4QIonIonCollision::Constr: CreationOfDiffractivePartonPairs, i="<<i<<G4endl;
00519 #endif
00520
00521 G4QHadron* aProjectile = anIniteraction->GetProjectile();
00522 G4QParton* aParton = aProjectile->GetNextParton();
00523 if (aParton)
00524 {
00525 aPartonPair = new G4QPartonPair(aParton, aProjectile->GetNextAntiParton(),
00526 G4QPartonPair::DIFFRACTIVE,
00527 G4QPartonPair::PROJECTILE);
00528 thePartonPairs.push_back(aPartonPair);
00529 #ifdef debug
00530 G4cout<<"G4QIonIonCollision::Constr: proj Diffractive PartonPair is filled"<<G4endl;
00531 #endif
00532 }
00533
00534 G4QHadron* aTarget = anIniteraction->GetTarget();
00535 aParton = aTarget->GetNextParton();
00536 if (aParton)
00537 {
00538 aPartonPair = new G4QPartonPair(aParton, aTarget->GetNextAntiParton(),
00539 G4QPartonPair::DIFFRACTIVE, G4QPartonPair::TARGET);
00540 thePartonPairs.push_back(aPartonPair);
00541 #ifdef debug
00542 G4cout<<"G4QIonIonCollision::Constr: target Diffractive PartonPair's filled"<<G4endl;
00543 #endif
00544 }
00545 }
00546 #ifdef debug
00547 G4cout<<"G4QIonIonCollision::Construct: DiffractivePartonPairs are created"<<G4endl;
00548 #endif
00549
00550
00551
00552 std::for_each(theInteractions.begin(),theInteractions.end(), DeleteQInteraction());
00553 theInteractions.clear();
00554 #ifdef debug
00555 G4cout<<"G4QIonIonCollision::Construct: Temporary objects are cleaned up"<<G4endl;
00556 #endif
00557
00558 theProjNucleus.DoLorentzBoost(theCMVelocity);
00559 theTargNucleus.DoLorentzBoost(theCMVelocity);
00560
00561 #ifdef debug
00562 G4cout<<">>........>>G4QIonIonCollision::Construct: >>..>> Strings are created "<<G4endl;
00563 #endif
00564 G4QPartonPair* aPair;
00565 G4QString* aString=0;
00566 while(thePartonPairs.size())
00567 {
00568 aPair = thePartonPairs.back();
00569 thePartonPairs.pop_back();
00570 #ifdef debug
00571 G4cout<<"G4QIonIonCollision::Construct:StringType="<<aPair->GetCollisionType()<<G4endl;
00572 #endif
00573 aString= new G4QString(aPair);
00574 #ifdef debug
00575 G4cout<<"G4QIonIonCollision::Construct:NewString4M="<<aString->Get4Momentum()<<G4endl;
00576 #endif
00577 aString->Boost(theCMVelocity);
00578 strings.push_back(aString);
00579 stringsInitted=true;
00580 delete aPair;
00581 }
00582 #ifdef edebug
00583 G4LorentzVector sum=theProjNucleus.Get4Momentum()+theTargNucleus.Get4Momentum();
00584 G4int rChg=totChg-theProjNucleus.GetZ()-theTargNucleus.GetZ();
00585 G4int rBaN=totBaN-theProjNucleus.GetA()-theTargNucleus.GetA();
00586 G4int nStrings=strings.size();
00587 G4cout<<"-EMC-G4QIonIonCollision::Constr:#ofString="<<nStrings<<",tNuc4M="<<sum<<G4endl;
00588 for(G4int i=0; i<nStrings; i++)
00589 {
00590 G4QString* prString=strings[i];
00591 G4LorentzVector strI4M=prString->Get4Momentum();
00592 sum+=strI4M;
00593 G4int sChg=prString->GetCharge();
00594 G4int sBaN=prString->GetBaryonNumber();
00595 G4int LPDG=prString->GetLeftParton()->GetPDGCode();
00596 G4int RPDG=prString->GetRightParton()->GetPDGCode();
00597 G4QContent LQC =prString->GetLeftParton()->GetQC();
00598 G4QContent RQC =prString->GetRightParton()->GetQC();
00599 rChg-=sChg;
00600 rBaN-=sBaN;
00601 G4cout<<"-EMC-G4QIonIonCollision::Construct: String#"<<i<<", 4M="<<strI4M<<", LPDG="
00602 <<LPDG<<LQC<<",RPDG="<<RPDG<<RQC<<", Ch="<<sChg<<", BN="<<sBaN<<G4endl;
00603 }
00604 G4cout<<"-EMC-G4QInel::Constr: r4M="<<sum-totZLS4M<<", rC="<<rChg<<", rB="<<rBaN<<G4endl;
00605 #endif
00606 if(!stringsInitted)
00607 {
00608 G4cerr<<"*****G4QIonIonCollision::Construct:**** No strings are created *****"<<G4endl;
00609 G4Exception("G4QIonIonCollision::Construct:","72",FatalException,"No Strings created");
00610 }
00611 #ifdef debug
00612 G4cout<<"G4QIonIonCollision::Constr:BeforeRotation, #0fStrings="<<strings.size()<<G4endl;
00613 #endif
00614
00615
00616
00617 for(unsigned astring=0; astring < strings.size(); astring++)
00618 strings[astring]->LorentzRotate(toLab);
00619 theProjNucleus.DoLorentzRotation(toLab);
00620 theTargNucleus.DoLorentzRotation(toLab);
00621
00622 #ifdef edebug
00623 G4LorentzVector sm=theProjNucleus.Get4Momentum()+theTargNucleus.Get4Momentum();
00624 G4int rCg=totChg-theProjNucleus.GetZ()-theTargNucleus.GetZ();
00625 G4int rBC=totBaN-theProjNucleus.GetA()-theTargNucleus.GetA();
00626 G4int nStrs=strings.size();
00627 G4cout<<"-EMCLS-G4QIonIonCollision::Constr:#ofS="<<nStrings<<",Nuc4M(E=M)="<<sum<<G4endl;
00628 for(G4int i=0; i<nStrs; i++)
00629 {
00630 G4LorentzVector strI4M=strings[i]->Get4Momentum();
00631 sm+=strI4M;
00632 G4int sChg=strings[i]->GetCharge();
00633 rCg-=sChg;
00634 G4int sBaN=strings[i]->GetBaryonNumber();
00635 rBC-=sBaN;
00636 G4cout<<"-EMCLS-G4QInel::Construct:String#"<<i<<",4M="<<strI4M<<strI4M.m()<<",Charge="
00637 <<sChg<<",BaryN="<<sBaN<<G4endl;
00638 }
00639 G4cout<<"-EMCLS-...G4QInel::Constr: r4M="<<sm-totLS4M<<", rC="<<rCg<<",rB="<<rBC<<G4endl;
00640 #endif
00641
00642
00643
00644 SwapPartons();
00645
00646
00647
00648 G4int problem=0;
00649 G4QStringVector::iterator ist;
00650 G4bool con=true;
00651 while(con && strings.size())
00652 {
00653 for(ist = strings.begin(); ist < strings.end(); ++ist)
00654 {
00655 G4bool bad=true;
00656 G4LorentzVector cS4M=(*ist)->Get4Momentum();
00657 G4double cSM2=cS4M.m2();
00658 G4QParton* cLeft=(*ist)->GetLeftParton();
00659 G4QParton* cRight=(*ist)->GetRightParton();
00660 G4int cLT=cLeft->GetType();
00661 G4int cRT=cRight->GetType();
00662 G4int cLPDG=cLeft->GetPDGCode();
00663 G4int cRPDG=cRight->GetPDGCode();
00664 G4int aLPDG=0;
00665 G4int aRPDG=0;
00666 if (cLPDG > 7) aLPDG= cLPDG/100;
00667 else if(cLPDG <-7) aLPDG=(-cLPDG)/100;
00668 if (cRPDG > 7) aRPDG= cRPDG/100;
00669 else if(cRPDG <-7) aRPDG=(-cRPDG)/100;
00670 G4int L1=0;
00671 G4int L2=0;
00672 if(aLPDG)
00673 {
00674 L1=aLPDG/10;
00675 L2=aLPDG%10;
00676 }
00677 G4int R1=0;
00678 G4int R2=0;
00679 if(aRPDG)
00680 {
00681 R1=aRPDG/10;
00682 R2=aRPDG%10;
00683 }
00684 G4double cSM=cSM2;
00685 if(cSM2>0.) cSM=std::sqrt(cSM2);
00686 #ifdef debug
00687 G4cout<<"G4QIonIonCollision::Construct: Needs Fusion? cLPDG="<<cLPDG<<",cRPDG="<<cRPDG
00688 <<",cM(cM2 if neg)="<<cSM<<G4endl;
00689 #endif
00690 if(cSM>0.)
00691 {
00692 G4bool single=true;
00693 G4double miM=0.;
00694 if(cLT==2 && cRT==2)
00695 {
00696 if(L1!=R1 && L1!=R2 && L2!=R1 && L2!=R2)
00697 {
00698 single=false;
00699 G4QPDGCode tmp;
00700 std::pair<G4int,G4int> paB=tmp.MakeTwoBaryons(L1, L2, R1, R2);
00701 miM=G4QPDGCode(paB.first).GetMass()+G4QPDGCode(paB.second).GetMass();
00702 }
00703 }
00704 if(single) miM=G4QPDGCode((*ist)->GetQC().GetSPDGCode()).GetMass() + mPi0;
00705
00706 #ifdef debug
00707 G4cout<<"G4QInel::Const:*IsItGood? realM="<<std::sqrt(cSM2)<<" > GSM="<<miM<<G4endl;
00708 #endif
00709 if(std::sqrt(cSM2) > miM) bad=false;
00710 }
00711 if(bad)
00712 {
00713 #ifdef debug
00714 G4cout<<"G4QInel::Const:TryFuse,L1="<<L1<<",L2="<<L2<<",R1="<<R1<<",R2="<<R2<<G4endl;
00715 #endif
00716 G4int cST=cLT+cRT;
00717 G4double excess=-DBL_MAX;
00718 G4double maxiM2=-DBL_MAX;
00719 G4QStringVector::iterator sst;
00720 G4QStringVector::iterator pst;
00721 G4int sLPDG=0;
00722 G4int sRPDG=0;
00723 G4int sOrd=0;
00724 G4bool minC=true;
00725 if(cSM2>0.) minC=false;
00726 for(pst = strings.begin(); pst < strings.end(); pst++) if(pst != ist)
00727 {
00728 G4LorentzVector pS4M=(*pst)->Get4Momentum()+cS4M;
00729 G4int nLPDG=0;
00730 G4int nRPDG=0;
00731 G4double pExcess=-DBL_MAX;
00732 G4double pSM2=pS4M.m2();
00733 #ifdef debug
00734 G4cout<<"->G4QIonIonCollision::Construct: sum4M="<<pS4M<<", M2="<<pSM2<<G4endl;
00735 #endif
00736
00737
00738 G4QParton* pLeft=(*pst)->GetLeftParton();
00739 G4QParton* pRight=(*pst)->GetRightParton();
00740 G4int pLT=pLeft->GetType();
00741 G4int pRT=pRight->GetType();
00742 G4int pLPDG=pLeft->GetPDGCode();
00743 G4int pRPDG=pRight->GetPDGCode();
00744 G4int LPDG=0;
00745 G4int RPDG=0;
00746 if (pLPDG > 7) LPDG= pLPDG/100;
00747 else if(pLPDG <-7) LPDG=(-pLPDG)/100;
00748 if (pRPDG > 7) RPDG= pRPDG/100;
00749 else if(pRPDG <-7) RPDG=(-pRPDG)/100;
00750 G4int pL1=0;
00751 G4int pL2=0;
00752 if(LPDG)
00753 {
00754 pL1=LPDG/10;
00755 pL2=LPDG%10;
00756 }
00757 G4int pR1=0;
00758 G4int pR2=0;
00759 if(RPDG)
00760 {
00761 pR1=RPDG/10;
00762 pR2=RPDG%10;
00763 }
00764 G4int pST=pLT+pRT;
00765 #ifdef debug
00766 G4cout<<"G4QInel::Construct: Partner/w pLPDG="<<pLPDG<<", pRPDG="<<pRPDG<<", pM2="
00767 <<pSM2<<G4endl;
00768 #endif
00769
00770 G4int tf=0;
00771 G4int af=0;
00772 if (cST==2 && pST==2)
00773 {
00774 tf=1;
00775 af=1;
00776 }
00777 else if(cST==2 && pST==3)
00778 {
00779 tf=2;
00780 if (pLPDG > 7 &&
00781 ( (cLPDG<0 && (-cLPDG==pL1 || -cLPDG==pL2 || -cLPDG==pRPDG) ) ||
00782 (cRPDG<0 && (-cRPDG==pL1 || -cRPDG==pL2 || -cRPDG==pRPDG) )
00783 )
00784 ) af=1;
00785 else if(pRPDG > 7 &&
00786 ( (cLPDG<0 && (-cLPDG==pR1 || -cLPDG==pR2 || -cLPDG==pLPDG) ) ||
00787 (cRPDG<0 && (-cRPDG==pR1 || -cRPDG==pR2 || -cRPDG==pLPDG) )
00788 )
00789 ) af=2;
00790 else if(pLPDG <-7 &&
00791 ( (cLPDG>0 && ( cLPDG==pL1 || cLPDG==pL2 || cLPDG==-pRPDG) ) ||
00792 (cRPDG>0 && ( cRPDG==pL1 || cRPDG==pL2 || cRPDG==-pRPDG) )
00793 )
00794 ) af=3;
00795 else if(pRPDG <-7 &&
00796 ( (cLPDG>0 && ( cLPDG==pR1 || cLPDG==pR2 || cLPDG==-pLPDG) ) ||
00797 (cRPDG>0 && ( cRPDG==pR1 || cRPDG==pR2 || cRPDG==-pLPDG) )
00798 )
00799 ) af=4;
00800 #ifdef debug
00801 else G4cout<<"G4QIonIonCollision::Constr:2(QaQ+QDiQ/aQaDiQ) Can't fuse"<<G4endl;
00802 #endif
00803 }
00804 else if(cST==3 && pST==2)
00805 {
00806 tf=3;
00807 if (cLPDG > 7 &&
00808 ( (pLPDG<0 && (-pLPDG==L1 || -pLPDG==L2) ) ||
00809 (pRPDG<0 && (-pRPDG==L1 || -pRPDG==L2) )
00810 )
00811 ) af=1;
00812 else if(cRPDG > 7 &&
00813 ( (pLPDG<0 && (-pLPDG==R1 || -pLPDG==R2) ) ||
00814 (pRPDG<0 && (-pRPDG==R1 || -pRPDG==R2) )
00815 )
00816 ) af=2;
00817 else if(cLPDG <-7 &&
00818 ( (pLPDG>0 && ( pLPDG==L1 || pLPDG==L2) ) ||
00819 (pRPDG>0 && ( pRPDG==L1 || pRPDG==L2) )
00820 )
00821 ) af=3;
00822 else if(cRPDG <-7 &&
00823 ( (pLPDG>0 && ( pLPDG==R1 || pLPDG==R2) ) ||
00824 (pRPDG>0 && ( pRPDG==R1 || pRPDG==R2) )
00825 )
00826 ) af=4;
00827 #ifdef debug
00828 else G4cout<<"G4QIonIonCollision::Constr:3(QDiQ/aQaDiQ+QaQ) Can't fuse"<<G4endl;
00829 #endif
00830 }
00831 else if(cST==2 && pST==4)
00832 {
00833 tf=4;
00834 if (pLPDG > 7)
00835 {
00836 if ( (-cLPDG==pL1 || -cLPDG==pL2) && (cRPDG==pR1 || cRPDG==pR2) ) af=1;
00837 else if( (-cRPDG==pL1 || -cRPDG==pL2) && (cLPDG==pR1 || cLPDG==pR2) ) af=2;
00838 }
00839 else if(pRPDG > 7)
00840 {
00841 if ( (-cRPDG==pR1 || -cRPDG==pR2) && (cLPDG==pL1 || cLPDG==pL2) ) af=3;
00842 else if( (-cLPDG==pR1 || -cLPDG==pR2) && (cRPDG==pL1 || cRPDG==pL2) ) af=4;
00843 }
00844 #ifdef debug
00845 else G4cout<<"-G4QIonIonCollision::Constr: 4 (QaQ+aQDiQDiQ) Can't fuse"<<G4endl;
00846 #endif
00847 }
00848 else if(cST==4 && pST==2)
00849 {
00850 tf=5;
00851 if (cLPDG > 7)
00852 {
00853 if ( (-pLPDG==L1 || -pLPDG==L2) && (pRPDG==R1 || pRPDG==R2) ) af=1;
00854 else if( (-pRPDG==L1 || -pRPDG==L2) && (pLPDG==R1 || pLPDG==R2) ) af=2;
00855 }
00856 else if(cRPDG > 7)
00857 {
00858 if ( (-pRPDG==R1 || -pRPDG==R2) && (pLPDG==L1 || pLPDG==L2) ) af=3;
00859 else if( (-pLPDG==R1 || -pLPDG==R2) && (pRPDG==L1 || pRPDG==L2) ) af=4;
00860 }
00861 #ifdef debug
00862 else G4cout<<"-G4QIonIonCollision::Constr: 5 (aQDiQDiQ+QaQ) Can't fuse"<<G4endl;
00863 #endif
00864 }
00865 else if(cST==3 && pST==3)
00866 {
00867 tf=6;
00868 if(pLPDG > 7)
00869 {
00870 if (cLPDG<-7 && (-cRPDG==pL1 || -cRPDG==pL2) && (pRPDG==L1 || pRPDG==L2))
00871 af=1;
00872 else if(cRPDG<-7 && (-cLPDG==pL1 || -cLPDG==pL2) && (pRPDG==R1 || pRPDG==R2))
00873 af=2;
00874 }
00875 else if(pRPDG > 7)
00876 {
00877 if (cLPDG<-7 && (-cRPDG==pR1 || -cRPDG==pR2) && (pLPDG==L1 || pLPDG==L2))
00878 af=3;
00879 else if(cRPDG<-7 && (-cLPDG==pR1 || -cLPDG==pR2) && (pLPDG==R1 || pLPDG==R2))
00880 af=4;
00881 }
00882 else if(cLPDG > 7)
00883 {
00884 if (pLPDG<-7 && (-pRPDG==L1 || -pRPDG==L2) && (cRPDG==pL1 || cRPDG==pL2))
00885 af=5;
00886 else if(pRPDG<-7 && (-pLPDG==L1 || -pLPDG==L2) && (cRPDG==pR1 || cRPDG==pR2))
00887 af=6;
00888 }
00889 else if(cRPDG > 7)
00890 {
00891 if (pLPDG<-7 && (-pRPDG==R1 || -pRPDG==R2) && (cLPDG==pL1 || cLPDG==pL2))
00892 af=7;
00893 else if(pRPDG<-7 && (-pLPDG==R1 || -pLPDG==R2) && (cLPDG==pR1 || cLPDG==pR2))
00894 af=8;
00895 }
00896 #ifdef debug
00897 else G4cout<<"-G4QIonIonCollision::Constr: 6 (QDiQ+aQaDiQ) Can't fuse"<<G4endl;
00898 #endif
00899 }
00900 #ifdef debug
00901 G4cout<<"G4QIonIonCollision::Constr: **Possibility**, tf="<<tf<<",af="<<af<<G4endl;
00902 #endif
00903 if(tf && af)
00904 {
00905
00906 G4int order=0;
00907 switch (tf)
00908 {
00909 case 1:
00910 if (cLPDG > 0 && pLPDG > 0)
00911 {
00912 order= 1;
00913 if (cLPDG > pLPDG) nLPDG=cLPDG*1000+pLPDG*100+1;
00914 else if(cLPDG < pLPDG) nLPDG=pLPDG*1000+cLPDG*100+1;
00915 else nLPDG=pLPDG*1000+cLPDG*100+3;
00916 if (cRPDG < pRPDG) nRPDG=cRPDG*1000+pRPDG*100-1;
00917 else if(cRPDG > pRPDG) nRPDG=pRPDG*1000+cRPDG*100-1;
00918 else nRPDG=pRPDG*1000+cRPDG*100-3;
00919 }
00920 else if(cLPDG < 0 && pLPDG < 0)
00921 {
00922 order= 1;
00923 if (cRPDG > pRPDG) nRPDG=cRPDG*1000+pRPDG*100+1;
00924 else if(cRPDG < pRPDG) nRPDG=pRPDG*1000+cRPDG*100+1;
00925 else nRPDG=pRPDG*1000+cRPDG*100+3;
00926 if (cLPDG < pLPDG) nLPDG=cLPDG*1000+pLPDG*100-1;
00927 else if(cLPDG > pLPDG) nLPDG=pLPDG*1000+cLPDG*100-1;
00928 else nLPDG=pLPDG*1000+cLPDG*100-3;
00929 }
00930 else if(cRPDG > 0 && pLPDG > 0)
00931 {
00932 order=-1;
00933 if (cRPDG > pLPDG) nLPDG=cRPDG*1000+pLPDG*100+1;
00934 else if(cRPDG < pLPDG) nLPDG=pLPDG*1000+cRPDG*100+1;
00935 else nLPDG=pLPDG*1000+cRPDG*100+3;
00936 if (cLPDG < pRPDG) nRPDG=cLPDG*1000+pRPDG*100-1;
00937 else if(cLPDG > pRPDG) nRPDG=pRPDG*1000+cLPDG*100-1;
00938 else nRPDG=pRPDG*1000+cLPDG*100-3;
00939 }
00940 else if(cRPDG < 0 && pLPDG < 0)
00941 {
00942 order=-1;
00943 if (cLPDG > pRPDG) nRPDG=cLPDG*1000+pRPDG*100+1;
00944 else if(cLPDG < pRPDG) nRPDG=pRPDG*1000+cLPDG*100+1;
00945 else nRPDG=pRPDG*1000+cLPDG*100+3;
00946 if (cRPDG < pLPDG) nLPDG=cRPDG*1000+pLPDG*100-1;
00947 else if(cRPDG > pLPDG) nLPDG=pLPDG*1000+cRPDG*100-1;
00948 else nLPDG=pLPDG*1000+cRPDG*100-3;
00949 }
00950 break;
00951 case 2:
00952 switch (af)
00953 {
00954 case 1:
00955 if(cLPDG < 0)
00956 {
00957 order= 1;
00958 if(-cLPDG==pRPDG)
00959 {
00960 nLPDG=pLPDG;
00961 nRPDG=cRPDG;
00962 }
00963 else
00964 {
00965 if (cRPDG > pRPDG) nRPDG=cRPDG*1000+pRPDG*100+1;
00966 else if(cRPDG < pRPDG) nRPDG=pRPDG*1000+cRPDG*100+1;
00967 else nRPDG=pRPDG*1000+cRPDG*100+3;
00968 if (-cLPDG == pL1) nLPDG=pL2;
00969 else nLPDG=pL1;
00970 }
00971 }
00972 else
00973 {
00974 order=-1;
00975 if(-cRPDG==pRPDG)
00976 {
00977 nLPDG=pLPDG;
00978 nRPDG=cLPDG;
00979 }
00980 else
00981 {
00982 if (cLPDG > pRPDG) nRPDG=cLPDG*1000+pRPDG*100+1;
00983 else if(cLPDG < pRPDG) nRPDG=pRPDG*1000+cLPDG*100+1;
00984 else nRPDG=pRPDG*1000+cLPDG*100+3;
00985 if (-cRPDG == pL1) nLPDG=pL2;
00986 else nLPDG=pL1;
00987 }
00988 }
00989 break;
00990 case 2:
00991 if(cLPDG < 0)
00992 {
00993 order=-1;
00994 if(-cLPDG==pLPDG)
00995 {
00996 nLPDG=cRPDG;
00997 nRPDG=pRPDG;
00998 }
00999 else
01000 {
01001 if (cRPDG > pLPDG) nLPDG=cRPDG*1000+pLPDG*100+1;
01002 else if(cRPDG < pLPDG) nLPDG=pLPDG*1000+cRPDG*100+1;
01003 else nLPDG=pLPDG*1000+cRPDG*100+3;
01004 if (-cLPDG == pR1) nRPDG=pR2;
01005 else nRPDG=pR1;
01006 }
01007 }
01008 else
01009 {
01010 order= 1;
01011 if(-cRPDG==pLPDG)
01012 {
01013 nLPDG=cLPDG;
01014 nRPDG=pRPDG;
01015 }
01016 else
01017 {
01018 if (cLPDG > pLPDG) nLPDG=cLPDG*1000+pLPDG*100+1;
01019 else if(cLPDG < pLPDG) nLPDG=pLPDG*1000+cLPDG*100+1;
01020 else nLPDG=pLPDG*1000+cLPDG*100+3;
01021 if (-cRPDG == pR1) nRPDG=pR2;
01022 else nRPDG=pR1;
01023 }
01024 }
01025 break;
01026 case 3:
01027 if(cLPDG > 0)
01028 {
01029 order= 1;
01030 if(cLPDG==-pRPDG)
01031 {
01032 nLPDG=pLPDG;
01033 nRPDG=cRPDG;
01034 }
01035 else
01036 {
01037 if (cRPDG < pRPDG) nRPDG=cRPDG*1000+pRPDG*100-1;
01038 else if(cRPDG > pRPDG) nRPDG=pRPDG*1000+cRPDG*100-1;
01039 else nRPDG=pRPDG*1000+cRPDG*100-3;
01040 if ( cLPDG == pL1) nLPDG=-pL2;
01041 else nLPDG=-pL1;
01042 }
01043 }
01044 else
01045 {
01046 order=-1;
01047 if(cRPDG==-pRPDG)
01048 {
01049 nLPDG=pLPDG;
01050 nRPDG=cLPDG;
01051 }
01052 else
01053 {
01054 if (cLPDG < pRPDG) nRPDG=cLPDG*1000+pRPDG*100-1;
01055 else if(cLPDG > pRPDG) nRPDG=pRPDG*1000+cLPDG*100-1;
01056 else nRPDG=pRPDG*1000+cLPDG*100-3;
01057 if ( cRPDG == pL1) nLPDG=-pL2;
01058 else nLPDG=-pL1;
01059 }
01060 }
01061 break;
01062 case 4:
01063 if(cLPDG > 0)
01064 {
01065 order=-1;
01066 if(cLPDG==-pLPDG)
01067 {
01068 nLPDG=cRPDG;
01069 nRPDG=pRPDG;
01070 }
01071 else
01072 {
01073 if (cRPDG < pLPDG) nLPDG=cRPDG*1000+pLPDG*100-1;
01074 else if(cRPDG > pLPDG) nLPDG=pLPDG*1000+cRPDG*100-1;
01075 else nLPDG=pLPDG*1000+cRPDG*100-3;
01076 if ( cLPDG == pR1) nRPDG=-pR2;
01077 else nRPDG=-pR1;
01078 }
01079 }
01080 else
01081 {
01082 order= 1;
01083 if(cRPDG==-pLPDG)
01084 {
01085 nLPDG=cLPDG;
01086 nRPDG=pRPDG;
01087 }
01088 else
01089 {
01090 if (cLPDG < pLPDG) nLPDG=cLPDG*1000+pLPDG*100-1;
01091 else if(cLPDG > pLPDG) nLPDG=pLPDG*1000+cLPDG*100-1;
01092 else nLPDG=pLPDG*1000+cLPDG*100-3;
01093 if ( cRPDG == pR1) nRPDG=-pR2;
01094 else nRPDG=-pR1;
01095 }
01096 }
01097 break;
01098 }
01099 break;
01100 case 3:
01101 switch (af)
01102 {
01103 case 1:
01104 if(pLPDG < 0)
01105 {
01106 order= 1;
01107 if (pRPDG > cRPDG) nRPDG=pRPDG*1000+cRPDG*100+1;
01108 else if(pRPDG < cRPDG) nRPDG=cRPDG*1000+pRPDG*100+1;
01109 else nRPDG=cRPDG*1000+pRPDG*100+3;
01110 if (-pLPDG == L1) nLPDG=L2;
01111 else nLPDG=L1;
01112 }
01113 else
01114 {
01115 order=-1;
01116 if (pLPDG > cRPDG) nLPDG=pLPDG*1000+cRPDG*100+1;
01117 else if(pLPDG < cRPDG) nLPDG=cRPDG*1000+pLPDG*100+1;
01118 else nLPDG=cRPDG*1000+pLPDG*100+3;
01119 if (-pRPDG == L1) nRPDG=L2;
01120 else nRPDG=L1;
01121 }
01122 break;
01123 case 2:
01124 if(pLPDG < 0)
01125 {
01126 order=-1;
01127 if (pRPDG > cLPDG) nRPDG=pRPDG*1000+cLPDG*100+1;
01128 else if(pRPDG < cLPDG) nRPDG=cLPDG*1000+pRPDG*100+1;
01129 else nRPDG=cLPDG*1000+pRPDG*100+3;
01130 if (-pLPDG == R1) nLPDG=R2;
01131 else nLPDG=R1;
01132 }
01133 else
01134 {
01135 order= 1;
01136 if (pLPDG > cLPDG) nLPDG=pLPDG*1000+cLPDG*100+1;
01137 else if(pLPDG < cLPDG) nLPDG=cLPDG*1000+pLPDG*100+1;
01138 else nLPDG=cLPDG*1000+pLPDG*100+3;
01139 if (-pRPDG == R1) nRPDG=R2;
01140 else nRPDG=R1;
01141 }
01142 break;
01143 case 3:
01144 if(pLPDG > 0)
01145 {
01146 order= 1;
01147 if (pRPDG < cRPDG) nRPDG=pRPDG*1000+cRPDG*100-1;
01148 else if(pRPDG > cRPDG) nRPDG=cRPDG*1000+pRPDG*100-1;
01149 else nRPDG=cRPDG*1000+pRPDG*100-3;
01150 if ( pLPDG == L1) nLPDG=-L2;
01151 else nLPDG=-L1;
01152 }
01153 else
01154 {
01155 order=-1;
01156 if (pLPDG < cRPDG) nLPDG=pLPDG*1000+cRPDG*100-1;
01157 else if(pLPDG > cRPDG) nLPDG=cRPDG*1000+pLPDG*100-1;
01158 else nLPDG=cRPDG*1000+pLPDG*100-3;
01159 if ( pRPDG == L1) nRPDG=-L2;
01160 else nRPDG=-L1;
01161 }
01162 break;
01163 case 4:
01164 if(pLPDG > 0)
01165 {
01166 order=-1;
01167 if (pRPDG < cLPDG) nRPDG=pRPDG*1000+cLPDG*100-1;
01168 else if(pRPDG > cLPDG) nRPDG=cLPDG*1000+pRPDG*100-1;
01169 else nRPDG=cLPDG*1000+pRPDG*100-3;
01170 if ( pLPDG == R1) nLPDG=-R2;
01171 else nLPDG=-R1;
01172 }
01173 else
01174 {
01175 order= 1;
01176 if (pLPDG < cLPDG) nLPDG=pLPDG*1000+cLPDG*100-1;
01177 else if(pLPDG > cLPDG) nLPDG=cLPDG*1000+pLPDG*100-1;
01178 else nLPDG=cLPDG*1000+pLPDG*100-3;
01179 if ( pRPDG == R1) nRPDG=-R2;
01180 else nRPDG=-R1;
01181 }
01182 break;
01183 }
01184 break;
01185 case 4:
01186 switch (af)
01187 {
01188 case 1:
01189 order= 1;
01190 if(-cLPDG == pL1) nLPDG= pL2;
01191 else nLPDG= pL1;
01192 if( cRPDG == pR1) nRPDG=-pR2;
01193 else nRPDG=-pR1;
01194 break;
01195 case 2:
01196 order=-1;
01197 if(-cRPDG == pL1) nLPDG= pL2;
01198 else nLPDG= pL1;
01199 if( cLPDG == pR1) nRPDG=-pR2;
01200 else nRPDG=-pR1;
01201 break;
01202 case 3:
01203 order= 1;
01204 if( cLPDG == pL1) nLPDG=-pL2;
01205 else nLPDG=-pL1;
01206 if(-cRPDG == pR1) nRPDG= pR2;
01207 else nRPDG= pR1;
01208 break;
01209 case 4:
01210 order=-1;
01211 if( cRPDG == pL1) nLPDG=-pL2;
01212 else nLPDG=-pL1;
01213 if(-cLPDG == pR1) nRPDG= pR2;
01214 else nRPDG= pR1;
01215 break;
01216 }
01217 break;
01218 case 5:
01219 switch (af)
01220 {
01221 case 1:
01222 order= 1;
01223 if(-pLPDG == L1) nLPDG= L2;
01224 else nLPDG= L1;
01225 if( pRPDG == R1) nRPDG=-R2;
01226 else nRPDG=-R1;
01227 break;
01228 case 2:
01229 order=-1;
01230 if(-pRPDG == L1) nRPDG= L2;
01231 else nRPDG= L1;
01232 if( pLPDG == R1) nLPDG=-R2;
01233 else nLPDG=-R1;
01234 break;
01235 case 3:
01236 order= 1;
01237 if( pLPDG == L1) nLPDG=-L2;
01238 else nLPDG=-L1;
01239 if(-pRPDG == R1) nRPDG= R2;
01240 else nRPDG= R1;
01241 break;
01242 case 4:
01243 order=-1;
01244 if( pRPDG == L1) nRPDG=-L2;
01245 else nRPDG=-L1;
01246 if(-pLPDG == R1) nLPDG= R2;
01247 else nLPDG= R1;
01248 break;
01249 }
01250 break;
01251 case 6:
01252 switch (af)
01253 {
01254 case 1:
01255 order=-1;
01256 if(-cRPDG == pL1) nLPDG= pL2;
01257 else nLPDG= pL1;
01258 if( pRPDG == L1) nRPDG= -L2;
01259 else nRPDG= -L1;
01260 break;
01261 case 2:
01262 order= 1;
01263 if(-cLPDG == pL1) nLPDG= pL2;
01264 else nLPDG= pL1;
01265 if( pRPDG == R1) nRPDG= -R2;
01266 else nRPDG= -R1;
01267 break;
01268 case 3:
01269 order= 1;
01270 if(-cRPDG == pR1) nRPDG= pR2;
01271 else nRPDG= pR1;
01272 if( pLPDG == L1) nLPDG= -L2;
01273 else nLPDG= -L1;
01274 break;
01275 case 4:
01276 order=-1;
01277 if(-cLPDG == pR1) nRPDG= pR2;
01278 else nRPDG= pR1;
01279 if( pLPDG == R1) nLPDG= -R2;
01280 else nLPDG= -R1;
01281 break;
01282 case 5:
01283 order=-1;
01284 if(-pRPDG == L1) nRPDG= L2;
01285 else nRPDG= L1;
01286 if( cRPDG == pL1) nLPDG=-pL2;
01287 else nLPDG=-pL1;
01288 break;
01289 case 6:
01290 order= 1;
01291 if(-pLPDG == L1) nLPDG= L2;
01292 else nLPDG= L1;
01293 if( cRPDG == pR1) nRPDG=-pR2;
01294 else nRPDG=-pR1;
01295 break;
01296 case 7:
01297 order= 1;
01298 if(-pRPDG == R1) nRPDG= R2;
01299 else nRPDG= R1;
01300 if( cLPDG == pL1) nLPDG=-pL2;
01301 else nLPDG=-pL1;
01302 break;
01303 case 8:
01304 order=-1;
01305 if(-pLPDG == R1) nLPDG= R2;
01306 else nLPDG= R1;
01307 if( cLPDG == pR1) nRPDG=-pR2;
01308 else nRPDG=-pR1;
01309 break;
01310 }
01311 break;
01312 }
01313 if(!order) G4cerr<<"-Warning-G4QInel::Constr: t="<<tf<<", a="<<af<<", cL="<<cLPDG
01314 <<", cR="<<cRPDG<<", pL="<<pLPDG<<", pR="<<pRPDG<<G4endl;
01315 else
01316 {
01317
01318 G4int LT=1;
01319 if(std::abs(nLPDG) > 7) ++LT;
01320 G4int RT=1;
01321 if(std::abs(nRPDG) > 7) ++RT;
01322 G4double minM=0.;
01323 G4bool sing=true;
01324 if(cLT==2 && cRT==2)
01325 {
01326 aLPDG=0;
01327 aRPDG=0;
01328 if(cLPDG>0)
01329 {
01330 aLPDG=nLPDG/100;
01331 aRPDG=(-nRPDG)/100;
01332 }
01333 else
01334 {
01335 aRPDG=nRPDG/100;
01336 aLPDG=(-nLPDG)/100;
01337 }
01338 G4int nL1=aLPDG/10;
01339 G4int nL2=aLPDG%10;
01340 G4int nR1=aRPDG/10;
01341 G4int nR2=aRPDG%10;
01342 if(nL1!=nR1 && nL1!=nR2 && nL2!=nR1 && nL2!=nR2)
01343 {
01344 #ifdef debug
01345 G4cout<<"G4QIonIonCollis::Constr:aLPDG="<<aLPDG<<", aRPDG="<<aRPDG<<G4endl;
01346 #endif
01347 sing=false;
01348 G4QPDGCode tmp;
01349 std::pair<G4int,G4int> pB=tmp.MakeTwoBaryons(nL1, nL2, nR1, nR2);
01350 minM=G4QPDGCode(pB.first).GetMass()+G4QPDGCode(pB.second).GetMass();
01351 }
01352 }
01353 if(sing)
01354 {
01355 std::pair<G4int,G4int> newPair = std::make_pair(nLPDG,nRPDG);
01356 G4QContent newStQC(newPair);
01357 #ifdef debug
01358 G4cout<<"G4QIn::Con: LPDG="<<nLPDG<<",RPDG="<<nRPDG<<",QC="<<newStQC<<G4endl;
01359 #endif
01360 G4int minPDG=newStQC.GetSPDGCode();
01361 minM=G4QPDGCode(minPDG).GetMass() + mPi0;
01362 }
01363
01364 G4bool win=false;
01365 G4double pSM=0.;
01366 if(pSM2>0.) pSM=std::sqrt(pSM2);
01367 if(minC && pSM2 > maxiM2)
01368 {
01369 maxiM2=pSM2;
01370 win=true;
01371 }
01372 else if(!minC || pSM > minM)
01373 {
01374 pExcess=pSM-minM;
01375 if(minC || pExcess > excess)
01376 {
01377 minC=false;
01378 excess=pExcess;
01379 win=true;
01380 }
01381 }
01382 if(win)
01383 {
01384 sst=pst;
01385 sLPDG=nLPDG;
01386 sRPDG=nRPDG;
01387 sOrd=order;
01388 }
01389 }
01390 }
01391
01392 }
01393 if(sOrd)
01394 {
01395 G4LorentzVector cL4M=cLeft->Get4Momentum();
01396 G4LorentzVector cR4M=cRight->Get4Momentum();
01397 G4QParton* pLeft=(*sst)->GetLeftParton();
01398 G4QParton* pRight=(*sst)->GetRightParton();
01399 G4LorentzVector pL4M=pLeft->Get4Momentum();
01400 G4LorentzVector pR4M=pRight->Get4Momentum();
01401 #ifdef debug
01402 G4cout<<"G4QIonIonCollis::Constr:cS4M="<<cS4M<<" fused/w pS4M="<<pL4M+pR4M<<G4endl;
01403 #endif
01404 if(sOrd>0)
01405 {
01406 pL4M+=cL4M;
01407 pR4M+=cR4M;
01408 }
01409 else
01410 {
01411 pL4M+=cR4M;
01412 pR4M+=cL4M;
01413 }
01414 pLeft->SetPDGCode(sLPDG);
01415 pLeft->Set4Momentum(pL4M);
01416 pRight->SetPDGCode(sRPDG);
01417 pRight->Set4Momentum(pR4M);
01418 delete (*ist);
01419 strings.erase(ist);
01420 #ifdef debug
01421 G4LorentzVector ss4M=pL4M+pR4M;
01422 G4cout<<"G4QIonIonCollision::Constr:Created,4M="<<ss4M<<",m2="<<ss4M.m2()<<G4endl;
01423 #endif
01424 if( ist != strings.begin() )
01425 {
01426 ist--;
01427 con=false;
01428 #ifdef debug
01429 G4cout<<"G4QIonIonCollision::Construct: *** IST Decremented ***"<<G4endl;
01430 #endif
01431 }
01432 else
01433 {
01434 con=true;
01435 #ifdef debug
01436 G4cout<<"G4QIonIonCollision::Construct: *** IST Begin ***"<<G4endl;
01437 #endif
01438 break;
01439 }
01440 }
01441 else
01442 {
01443 #ifdef debug
01444 G4cout<<"-Warning-G4QInel::Const: S4M="<<cS4M<<",M2="<<cSM2<<" Leave ASIS"<<G4endl;
01445 #endif
01446 ++problem;
01447 con=false;
01448 }
01449 }
01450 else con=false;
01451 }
01452 #ifdef debug
01453 G4cout<<"G4QIonIonCollision::Construct: *** IST While *** , con="<<con<<G4endl;
01454 #endif
01455 }
01456 #ifdef edebug
01457
01458 G4LorentzVector t4M=theProjNucleus.Get4Momentum()+theTargNucleus.Get4Momentum();
01459 G4int rC=totChg-theProjNucleus.GetZ()-theTargNucleus.GetZ();
01460 G4int rB=totBaN-theProjNucleus.GetA()-theTargNucleus.GetA();
01461 G4int nStr=strings.size();
01462 G4cout<<"-EMCLS-G4QIn::Const: AfterSUPPRESION #ofS="<<nStr<<",tNuc4M(E=M)="<<sum<<G4endl;
01463 for(G4int i=0; i<nStr; i++)
01464 {
01465 G4LorentzVector strI4M=strings[i]->Get4Momentum();
01466 t4M+=strI4M;
01467 G4int sChg=strings[i]->GetCharge();
01468 rC-=sChg;
01469 G4int sBaN=strings[i]->GetBaryonNumber();
01470 rB-=sBaN;
01471 G4cout<<"-EMCLS-G4QIonIonCollision::Construct: St#"<<i<<", 4M="<<strI4M<<", M="
01472 <<strI4M.m()<<", C="<<sChg<<", B="<<sBaN<<G4endl;
01473 }
01474 G4cout<<"-EMCLS-G4QInel::Construct: r4M="<<t4M-totLS4M<<", rC="<<rC<<", rB="<<rB<<G4endl;
01475 #endif
01476
01477
01478
01479 #ifdef debug
01480 G4cout<<"G4QIonIonCollision::Construct: problem="<<problem<<G4endl;
01481 #endif
01482 if(problem)
01483 {
01484 G4int nOfStr=strings.size();
01485 #ifdef debug
01486 G4cout<<"G4QIonIonCollision::Constr:SecurityDiQaDiQReduction, #OfStr="<<nOfStr<<G4endl;
01487 #endif
01488 for (G4int astring=0; astring < nOfStr; astring++)
01489 {
01490 G4QString* curString=strings[astring];
01491 G4QParton* cLeft=curString->GetLeftParton();
01492 G4QParton* cRight=curString->GetRightParton();
01493 G4int LT=cLeft->GetType();
01494 G4int RT=cRight->GetType();
01495 G4int sPDG=cLeft->GetPDGCode();
01496 G4int nPDG=cRight->GetPDGCode();
01497 if(LT==2 && RT==2)
01498 {
01499 #ifdef debug
01500 G4cout<<"G4QIonIonCollis::Constr:TrySelfReduString,L="<<sPDG<<",R="<<nPDG<<G4endl;
01501 #endif
01502 if( cLeft->ReduceDiQADiQ(cLeft, cRight) )
01503 {
01504 sPDG=cLeft->GetPDGCode();
01505 nPDG=cRight->GetPDGCode();
01506 #ifdef debug
01507 G4cout<<"+G4QInel::Const:#"<<astring<<" Reduced, L="<<sPDG<<", R="<<nPDG<<G4endl;
01508 #endif
01509 }
01510 #ifdef debug
01511 else G4cout<<"--*--G4QInel::Const: #"<<astring<<" DQ-aDQ reduction Failed"<<G4endl;
01512 #endif
01513 }
01514 else if(sPDG==3 && nPDG==-3)
01515 {
01516 sPDG= 1;
01517 nPDG=-1;
01518 cLeft->SetPDGCode(sPDG);
01519 cRight->SetPDGCode(nPDG);
01520 }
01521 else if(sPDG==-3 && nPDG==3)
01522 {
01523 sPDG=-1;
01524 nPDG= 1;
01525 cLeft->SetPDGCode(sPDG);
01526 cRight->SetPDGCode(nPDG);
01527 }
01528 }
01529 SwapPartons();
01530 }
01531 #ifdef edebug
01532 G4LorentzVector u4M=theProjNucleus.Get4Momentum()+theTargNucleus.Get4Momentum();
01533 G4int rCh=totChg-theProjNucleus.GetZ()-theTargNucleus.GetZ();
01534 G4int rBa=totBaN-theProjNucleus.GetA()-theTargNucleus.GetA();
01535 G4int nStri=strings.size();
01536 G4cout<<"-EMCLS-G4QIn::Const: FinalConstruct, #ofSt="<<nStri<<",tN4M(E=M)="<<t4M<<G4endl;
01537 for(G4int i=0; i<nStri; i++)
01538 {
01539 G4LorentzVector strI4M=strings[i]->Get4Momentum();
01540 u4M+=strI4M;
01541 G4int sChg=strings[i]->GetCharge();
01542 rCh-=sChg;
01543 G4int sBaN=strings[i]->GetBaryonNumber();
01544 rBa-=sBaN;
01545 G4cout<<"-EMCLS-G4QIonIonCollision::Construct: St#"<<i<<", 4M="<<strI4M<<", M="
01546 <<strI4M.m()<<", C="<<sChg<<", B="<<sBaN<<G4endl;
01547 }
01548 G4cout<<"-EMCLS-G4QInel::Construct: r4M="<<u4M-totLS4M<<",rC="<<rCh<<",rB="<<rBa<<G4endl;
01549 #endif
01550 }
01551
01552 G4QIonIonCollision::~G4QIonIonCollision()
01553 {
01554 std::for_each(strings.begin(), strings.end(), DeleteQString() );
01555 }
01556
01557 G4QHadronVector* G4QIonIonCollision::Fragment()
01558 {
01559 #ifdef debug
01560 G4cout<<"*******>G4QIonIonCollision::Fragment: ***Called***, Res="<<theResult<<G4endl;
01561 #endif
01562 G4int striNum=strings.size();
01563 G4int hadrNum=theResult->size();
01564 #ifdef edebug
01565 G4int nQm=theQuasmons.size();
01566 G4LorentzVector totLS4M=theProjNucleus.Get4Momentum()+theTargNucleus.Get4Momentum();
01567 G4int totChg=theProjNucleus.GetZ()+theTargNucleus.GetZ();
01568 G4int totBaN=theTargNucleus.GetA()+theTargNucleus.GetA();
01569 G4cout<<"-EMCLS-G4QIn::Fragment: CHECKRecovery, #ofS="<<striNum<<", Nuc4M(E=M)="<<totLS4M
01570 <<",#Q="<<nQm<<",#H="<<hadrNum<<G4endl;
01571 for(G4int i=0; i < striNum; i++)
01572 {
01573 G4LorentzVector strI4M=strings[i]->Get4Momentum();
01574 totLS4M+=strI4M;
01575 G4int sChg=strings[i]->GetCharge();
01576 totChg+=sChg;
01577 G4int sBaN=strings[i]->GetBaryonNumber();
01578 totBaN+=sBaN;
01579 G4cout<<"-EMCLS-G4QIonIonCollision::Fragment:String#"<<i<<", 4M="<<strI4M<<", M="
01580 <<strI4M.m()<<", C="<<sChg<<", B="<<sBaN<<G4endl;
01581 }
01582 for(G4int i=0; i < nQm; i++)
01583 {
01584 G4LorentzVector hI4M=theQuasmons[i]->Get4Momentum();
01585 totLS4M+=hI4M;
01586 G4int hChg=theQuasmons[i]->GetCharge();
01587 totChg+=hChg;
01588 G4int hBaN=theQuasmons[i]->GetBaryonNumber();
01589 totBaN+=hBaN;
01590 G4cout<<"-EMCLS-G4QIonIonCollision::Fragment: Quasmon#"<<i<<", 4M="<<hI4M<<", C="<<hChg
01591 <<", B="<<hBaN<<G4endl;
01592 }
01593 for(G4int i=0; i < hadrNum; i++)
01594 {
01595 G4LorentzVector hI4M=(*theResult)[i]->Get4Momentum();
01596 totLS4M+=hI4M;
01597 G4int hChg=(*theResult)[i]->GetCharge();
01598 totChg+=hChg;
01599 G4int hBaN=(*theResult)[i]->GetBaryonNumber();
01600 totBaN+=hBaN;
01601 G4cout<<"-EMCLS-G4QIn::Fragment:H#"<<i<<",4M="<<hI4M<<",C="<<hChg<<",B="<<hBaN<<G4endl;
01602 }
01603 #endif
01604 #ifdef debug
01605 G4cout<<"***>G4QIonIonCollision::Fragm: #OfStr="<<striNum<<", #OfRes="<<hadrNum<<G4endl;
01606 #endif
01607 if(!striNum && hadrNum)
01608 {
01609 #ifdef debug
01610 G4cout<<"***>G4QIonIonCollision::Fragm:**Quasi-Elastic**, #OfResult="<<hadrNum<<G4endl;
01611 #endif
01612 return theResult;
01613 }
01614 else if(striNum) Breeder();
01615 else
01616 {
01617 if(hadrNum)
01618 {
01619 for(G4int ih=0; ih<hadrNum; ih++) delete (*theResult)[ih];
01620 theResult->clear();
01621 }
01622 G4LorentzVector rp4M=theProjNucleus.Get4Momentum();
01623 G4int rpPDG=theProjNucleus.GetPDG();
01624 G4QHadron* resPNuc = new G4QHadron(rpPDG,rp4M);
01625 theResult->push_back(resPNuc);
01626 G4LorentzVector rt4M=theTargNucleus.Get4Momentum();
01627 G4int rtPDG=theTargNucleus.GetPDG();
01628 G4QHadron* resTNuc = new G4QHadron(rtPDG,rt4M);
01629 theResult->push_back(resTNuc);
01630 }
01631 G4int nQuas=theQuasmons.size();
01632 G4int theRS=theResult->size();
01633 #ifdef debug
01634 G4cout<<"***>G4QIonIonCollision::Fragm:beforeEnv, #OfQ="<<nQuas<<",#OfR="<<theRS<<G4endl;
01635 #endif
01636 if(nQuas && theRS)
01637 {
01638
01639 G4QHadron* resNuc = (*theResult)[theRS-1];
01640 G4LorentzVector resNuc4M = resNuc->Get4Momentum();
01641 G4int resNucPDG= resNuc->GetPDGCode();
01642 G4QNucleus theEnv(resNucPDG);
01643 delete resNuc;
01644 theResult->pop_back();
01645 --theRS;
01646 #ifdef pdebug
01647 G4cout<<"G4QIonIonCollision::Fragm:#OfRemainingHadron="<<theRS<<", A="<<theEnv<<G4endl;
01648 #endif
01649
01650 for(G4int j=theRS-1; j>-2; --j)
01651 {
01652 G4LorentzVector qsum4M=resNuc4M;
01653 G4QContent qsumQC=theEnv.GetQCZNS();
01654 #ifdef pdebug
01655 G4cout<<"G4QIonIonCollis::Fragm: rN4M"<<qsum4M<<qsum4M.m()<<",rNQC="<<qsumQC<<G4endl;
01656 #endif
01657 G4Quasmon* firstQ=0;
01658 G4LorentzVector first4M;
01659 G4QContent firstQC;
01660 for(G4int i=0; i<nQuas; ++i)
01661 {
01662 G4Quasmon* curQuasm=theQuasmons[i];
01663 G4LorentzVector cur4M=curQuasm->Get4Momentum();
01664 G4QContent curQC=curQuasm->GetQC();
01665 qsum4M+=cur4M;
01666 qsumQC+=curQC;
01667 #ifdef pdebug
01668 G4cout<<"G4QIonIonCollis::Fragm: Q#"<<i<<", QQC="<<curQC<<", sQC="<<qsumQC<<G4endl;
01669 #endif
01670 if(!i)
01671 {
01672 firstQ =curQuasm;
01673 first4M=cur4M;
01674 firstQC=curQC;
01675 }
01676 }
01677 G4int miPDG=qsumQC.GetSPDGCode();
01678 G4double gsM=0.;
01679 if(miPDG == 10)
01680 {
01681 G4QChipolino QCh(qsumQC);
01682 gsM=QCh.GetQPDG1().GetMass()+QCh.GetQPDG2().GetMass();
01683
01684
01685 }
01686
01687
01688
01689
01690
01691
01692 else if(miPDG < 80000000 && std::abs(miPDG)%10 > 2)
01693 gsM=theWorld->GetQParticle(G4QPDGCode(miPDG))->MinMassOfFragm();
01694 else gsM=G4QPDGCode(miPDG).GetMass();
01695 G4double reM=qsum4M.m();
01696 #ifdef pdebug
01697 G4cout<<"G4QIonIonCollision::Fragm: PDG="<<miPDG<<", rM="<<reM<<",GSM="<<gsM<<G4endl;
01698 #endif
01699 if(reM > gsM) break;
01700 if(j > -1)
01701 {
01702 G4QHadron* cH = (*theResult)[j];
01703 G4LorentzVector h4M = cH->Get4Momentum();
01704 G4QContent hQC = cH->GetQC();
01705 firstQ->Set4Momentum(first4M+h4M);
01706 firstQ->SetQC(firstQC+hQC);
01707 delete cH;
01708 theResult->pop_back();
01709 #ifdef pdebug
01710 G4cout<<"G4QInel::Fragm: H#"<<j<<", hQC="<<hQC<<",hPDG="<<cH->GetPDGCode()<<G4endl;
01711 #endif
01712 }
01713 else
01714 {
01715 G4cerr<<"***G4QIonIonCollision::Fra:PDG="<<miPDG<<",M="<<reM<<",GSM="<<gsM<<G4endl;
01716 G4Exception("G4QIonIonCollision::Fragment:","27",FatalException,"Can'tRecoverGSM");
01717 }
01718 }
01719 G4double nucE=resNuc4M.e();
01720 if(nucE<1.E-12) nucE=0.;
01721 G4ThreeVector nucVel(0.,0.,0.);
01722 G4QHadronVector* output=0;
01723 G4QEnvironment* pan= new G4QEnvironment(theEnv);
01724 #ifdef pdebug
01725 G4cout<<"G4QIonIonCollision::Fragm: nucE="<<nucE<<", nQ="<<nQuas<<G4endl;
01726 #endif
01727 if(nucE) nucVel=resNuc4M.vect()/nucE;
01728 for(G4int i=0; i<nQuas; ++i)
01729 {
01730 G4Quasmon* curQuasm=theQuasmons[i];
01731 if(nucE) curQuasm->Boost(-nucVel);
01732 pan->AddQuasmon(curQuasm);
01733 #ifdef pdebug
01734 G4LorentzVector cQ4M=curQuasm->Get4Momentum();
01735 G4cout<<"G4QIonIonCollision::Fragment: Quasmon# "<<i<<" added, 4M="<<cQ4M<<G4endl;
01736 #endif
01737 }
01738 try
01739 {
01740 delete output;
01741 output = pan->Fragment();
01742 }
01743 catch (G4QException& error)
01744 {
01745 G4cerr<<"***G4QIonIonCollision::Fragment: G4QE Exception is catched"<<G4endl;
01746
01747 G4Exception("G4QIonIonCollision::Fragment()", "HAD_CHPS_0027",
01748 FatalException, "CHIPSCrash");
01749 }
01750 delete pan;
01751 if(output)
01752 {
01753 G4int nOut=output->size();
01754 for(G4int j=0; j<nOut; j++)
01755 {
01756 G4QHadron* curHadron=(*output)[j];
01757 if(nucE) curHadron->Boost(nucVel);
01758 theResult->push_back(curHadron);
01759 }
01760 delete output;
01761 }
01762 }
01763 else if(!striNum) G4cout<<"-Warning-G4QIonIonCollision::Fragment:NothingWasDone"<<G4endl;
01764 #ifdef debug
01765 G4cout<<"=--=>G4QIonIonCollision::Fragment: Final #OfResult="<<theResult->size()<<G4endl;
01766 #endif
01767 G4int nQ =theQuasmons.size();
01768 if(nQ) theQuasmons.clear();
01769 #ifdef edebug
01770 G4LorentzVector f4M(0.,0.,0.,0.);
01771 G4int fCh=totChg;
01772 G4int fBN=totBaN;
01773 G4int nHd=theResult->size();
01774 G4cout<<"-EMCLS-G4QIonIonCollision::Fragment: #ofHadr="<<nHd<<",#OfQuasm="<<nQ<<G4endl;
01775 for(G4int i=0; i<nHd; i++)
01776 {
01777 G4LorentzVector hI4M=(*theResult)[i]->Get4Momentum();
01778 f4M+=hI4M;
01779 G4int hChg=(*theResult)[i]->GetCharge();
01780 fCh-=hChg;
01781 G4int hBaN=(*theResult)[i]->GetBaryonNumber();
01782 fBN-=hBaN;
01783 G4cout<<"-EMCLS-G4QIonIonCollision::Fragment: Hadron#"<<i<<", 4M="<<hI4M<<", PDG="
01784 <<(*theResult)[i]->GetPDGCode()<<", C="<<hChg<<", B="<<hBaN<<G4endl;
01785 }
01786 G4cout<<"-EMCLS-G4QInel::Fragm:, r4M="<<f4M-totLS4M<<", rC="<<fCh<<",rB="<<fBN<<G4endl;
01787 #endif
01788 return theResult;
01789 }
01790
01791 void G4QIonIonCollision::Breeder()
01792 {
01793 static const G4double eps = 0.001;
01794
01795
01796
01797 #ifdef edebug
01798 G4LorentzVector totLS4M=theProjNucleus.Get4Momentum()+theTargNucleus.Get4Momentum();
01799 G4int totChg=theProjNucleus.GetZ()+theTargNucleus.GetZ();
01800 G4int totBaN=theProjNucleus.GetA()+theTargNucleus.GetA();
01801 G4int nStri=strings.size();
01802 G4cout<<"-EMCLS-G4QIn::Breed: CHECKRecovery #ofS="<<nStri<<",N4M(E=M)="<<totLS4M<<G4endl;
01803 for(G4int i=0; i<nStri; i++)
01804 {
01805 G4LorentzVector strI4M=strings[i]->Get4Momentum();
01806 totLS4M+=strI4M;
01807 G4int sChg=strings[i]->GetCharge();
01808 totChg+=sChg;
01809 G4int sBaN=strings[i]->GetBaryonNumber();
01810 totBaN+=sBaN;
01811 G4cout<<"-EMCLS-G4QIonIonCollision::Breeder: St#"<<i<<", 4M="<<strI4M<<", M="
01812 <<strI4M.m()<<", C="<<sChg<<", B="<<sBaN<<G4endl;
01813 }
01814 #endif
01815 G4int nOfStr=strings.size();
01816 #ifdef debug
01817 G4cout<<"G4QIonIonCollision::Breeder: BeforeFragmentation, #OfStr="<<nOfStr<<G4endl;
01818 #endif
01819 G4LorentzVector ft4M(0.,0.,0.,0.);
01820 G4QContent ftQC(0,0,0,0,0,0);
01821 G4bool ftBad=false;
01822 for(G4int i=0; i < nOfStr; ++i)
01823 {
01824 G4LorentzVector pS4M=strings[i]->Get4Momentum();
01825 ft4M+=pS4M;
01826 G4QContent pSQC=strings[i]->GetQC();
01827 ftQC+=pSQC;
01828 if(pS4M.m2() < 0.) ftBad=true;
01829 #ifdef debug
01830 G4cout<<">G4QIonIonCollision::Breed:1stTest,S#"<<i<<",4M="<<pS4M<<",QC="<<pSQC<<G4endl;
01831 #endif
01832 }
01833 if(ftBad)
01834 {
01835 G4Quasmon* stringQuasmon = new G4Quasmon(ftQC, ft4M);
01836 #ifdef debug
01837 G4cout<<"->G4QIonIonCollision::Breed:*TotQ*,QC="<<ftQC<<",4M="<<ft4M<<ft4M.m()<<G4endl;
01838 #endif
01839 theQuasmons.push_back(stringQuasmon);
01840 G4LorentzVector rp4M=theProjNucleus.Get4Momentum();
01841 G4int rpPDG=theProjNucleus.GetPDG();
01842 G4QHadron* resPNuc = new G4QHadron(rpPDG,rp4M);
01843 theResult->push_back(resPNuc);
01844 G4LorentzVector rt4M=theTargNucleus.Get4Momentum();
01845 G4int rtPDG=theTargNucleus.GetPDG();
01846 G4QHadron* resTNuc = new G4QHadron(rtPDG,rt4M);
01847 theResult->push_back(resTNuc);
01848 return;
01849 }
01850 for (G4int astring=0; astring < nOfStr; astring++)
01851 {
01852 #ifdef edebug
01853 G4LorentzVector sum=theProjNucleus.Get4Momentum()+theTargNucleus.Get4Momentum();
01854 G4int rChg=totChg-theProjNucleus.GetZ()-theTargNucleus.GetZ();
01855 G4int rBaN=totBaN-theProjNucleus.GetA()-theTargNucleus.GetA();
01856 G4int nOfHadr=theResult->size();
01857 G4cout<<"-EMCLS-G4QIonIonCollision::Breed:#ofSt="<<nOfStr<<",#ofHad="<<nOfHadr<<G4endl;
01858 for(G4int i=astring; i<nOfStr; i++)
01859 {
01860 G4LorentzVector strI4M=strings[i]->Get4Momentum();
01861 sum+=strI4M;
01862 G4int sChg=strings[i]->GetCharge();
01863 rChg-=sChg;
01864 G4int sBaN=strings[i]->GetBaryonNumber();
01865 rBaN-=sBaN;
01866 G4cout<<"-EMCLS-G4QI::Breed:S#"<<i<<",4M="<<strI4M<<",C="<<sChg<<",B="<<sBaN<<G4endl;
01867 }
01868 for(G4int i=0; i<nOfHadr; i++)
01869 {
01870 G4LorentzVector hI4M=(*theResult)[i]->Get4Momentum();
01871 sum+=hI4M;
01872 G4int hChg=(*theResult)[i]->GetCharge();
01873 rChg-=hChg;
01874 G4int hBaN=(*theResult)[i]->GetBaryonNumber();
01875 rBaN-=hBaN;
01876 G4cout<<"-EMCLS-G4QIn::Breed: H#"<<i<<",4M="<<hI4M<<",C="<<hChg<<",B="<<hBaN<<G4endl;
01877 }
01878 G4cout<<"....-EMCLS-G4QInel::Br:r4M="<<sum-totLS4M<<",rC="<<rChg<<",rB="<<rBaN<<G4endl;
01879 #endif
01880 G4QString* curString=strings[astring];
01881 if(!curString->GetDirection()) continue;
01882 #ifdef edebug
01883 G4int curStrChg = curString->GetCharge();
01884 G4int curStrBaN = curString->GetBaryonNumber();
01885 #endif
01886 G4LorentzVector curString4M = curString->Get4Momentum();
01887 #ifdef debug
01888 G4cout<<"=--=>G4QIonIonCollision::Breeder: String#"<<astring<<",s4M/m="<<curString4M
01889 <<curString4M.m()<<", LPart="<<curString->GetLeftParton()->GetPDGCode()
01890 <<", RPart="<<curString->GetRightParton()->GetPDGCode()<<G4endl;
01891 #endif
01892 G4QHadronVector* theHadrons = 0;
01893 theHadrons=curString->FragmentString(true);
01894 if (!theHadrons)
01895 {
01896
01897 G4QParton* cLeft=curString->GetLeftParton();
01898 G4QParton* cRight=curString->GetRightParton();
01899 G4int sPDG=cLeft->GetPDGCode();
01900 G4int nPDG=cRight->GetPDGCode();
01901 G4int LT=cLeft->GetType();
01902 G4int RT=cRight->GetType();
01903 G4int LS=LT+RT;
01904 if(LT==2 && RT==2)
01905 {
01906 #ifdef debug
01907 G4cout<<"G4QIonIonCollision::Breed:TryReduceString, L="<<sPDG<<",R="<<nPDG<<G4endl;
01908 #endif
01909 if( cLeft->ReduceDiQADiQ(cLeft, cRight) )
01910 {
01911 LT=1;
01912 RT=1;
01913 LS=2;
01914 sPDG=cLeft->GetPDGCode();
01915 nPDG=cRight->GetPDGCode();
01916 #ifdef debug
01917 G4cout<<"G4QIonIonCollision::Breed:AfterReduction,L="<<sPDG<<",R="<<nPDG<<G4endl;
01918 #endif
01919 theHadrons=curString->FragmentString(true);
01920 cLeft=curString->GetLeftParton();
01921 cRight=curString->GetRightParton();
01922 #ifdef debug
01923 G4cout<<"G4QInel::Breed:L="<<cLeft->Get4Momentum()<<",R="<<cRight->Get4Momentum()
01924 <<G4endl;
01925 #endif
01926 }
01927 #ifdef debug
01928 else G4cout<<"^G4QIonIonCollision::Breed: DQ-aDQ reduction to Q-aQ Failed"<<G4endl;
01929 #endif
01930 }
01931 #ifdef debug
01932 G4cout<<"G4QIonIonCollision::Breed:AfterRedAttempt, theH="<<theHadrons<<", L4M="
01933 <<cLeft->Get4Momentum()<<", R4M="<<cRight->Get4Momentum()<<G4endl;
01934 #endif
01935 unsigned next=astring+1;
01936 if (!theHadrons)
01937 {
01938 G4int fusionDONE=0;
01939 if(next < strings.size())
01940 {
01941 G4int fustr=0;
01942 G4int swap=0;
01943 G4double Vmin=DBL_MAX;
01944 G4int dPDG=nPDG;
01945 G4int qPDG=sPDG;
01946 if(dPDG<-99 || (dPDG>0&&dPDG<7) || qPDG>99 || (qPDG<0 && qPDG>-7))
01947 {
01948 swap=qPDG;
01949 qPDG=dPDG;
01950 dPDG=swap;
01951 }
01952 if(dPDG>99) dPDG/=100;
01953 if(qPDG<-99) qPDG=-(-qPDG)/100;
01954 #ifdef debug
01955 G4cout<<"G4QIonIonCollision::Breed:TryFuseStringS, q="<<qPDG<<", a="<<dPDG
01956 <<", n="<<next<<G4endl;
01957 #endif
01958 G4ThreeVector curV=curString4M.vect()/curString4M.e();
01959 G4int reduce=0;
01960 G4int restr=0;
01961 G4int MPS=0;
01962 for (restr=next; restr < nOfStr; restr++)
01963 {
01964 reduce=0;
01965 G4QString* reString=strings[restr];
01966 G4QParton* Left=reString->GetLeftParton();
01967 G4QParton* Right=reString->GetRightParton();
01968 G4int uPDG=Left->GetPDGCode();
01969 G4int mPDG=Right->GetPDGCode();
01970 G4int PLT =Left->GetType();
01971 G4int PRT =Right->GetType();
01972 G4int aPDG=mPDG;
01973 G4int rPDG=uPDG;
01974 if(aPDG<-99 || (aPDG>0 && aPDG<7) || rPDG>99 || (rPDG<0 && rPDG>-7))
01975 {
01976 swap=rPDG;
01977 rPDG=aPDG;
01978 aPDG=swap;
01979 }
01980 if(aPDG > 99) aPDG/=100;
01981 if(rPDG <-99) rPDG=-(-rPDG)/100;
01982
01983 #ifdef debug
01984 G4cout<<"G4Qnel::Breed: TryReduce #"<<restr<<", q="<<rPDG<<",a="<<aPDG<<G4endl;
01985 #endif
01986 if(LT==2 && RT==2 && PLT==2 && PRT==2)
01987 {
01988 G4int cQ1=(-qPDG)/10;
01989 G4int cQ2=(-qPDG)%10;
01990 G4int cA1=dPDG/10;
01991 G4int cA2=dPDG%10;
01992 G4int pQ1=(-rPDG)/10;
01993 G4int pQ2=(-rPDG)%10;
01994 G4int pA1=aPDG/10;
01995 G4int pA2=aPDG%10;
01996 #ifdef debug
01997 G4cout<<"G4QIonIonCollision::Breeder: cQ="<<cQ1<<","<<cQ2<<", cA="<<cA1<<","
01998 <<cA2<<", pQ="<<pQ1<<","<<pQ2<<", pA="<<pA1<<","<<pA2<<G4endl;
01999 #endif
02000 G4bool iQA = (cA1==pQ1 || cA1==pQ2 || cA2==pQ1 || cA2==pQ2);
02001 G4bool iAQ = (cQ1==pA1 || cQ1==pA2 || cQ2==pA1 || cQ2==pA2);
02002 if(iQA) reduce++;
02003 if(iAQ) reduce++;
02004 if (reduce==2)
02005 {
02006 if(sPDG>0 && uPDG<0)
02007 {
02008 std::pair<G4int,G4int> resLL=ReducePair(sPDG/100, (-uPDG)/100);
02009 G4int newCL=resLL.first;
02010 G4int newPL=resLL.second;
02011 if(!newCL || !newPL)
02012 {
02013 G4cerr<<"*G4QIonIonCollision::Breed:CL="<<newCL<<",PL="<<newPL<<G4endl;
02014 G4Exception("G4QIonIonCollision::Breeder:","72",FatalException,"2LL-");
02015 }
02016 std::pair<G4int,G4int> resRR=ReducePair((-nPDG)/100, mPDG/100);
02017 G4int newCR=resRR.first;
02018 G4int newPR=resRR.second;
02019 if(!newCR || !newPR)
02020 {
02021 G4cerr<<"*G4QIonIonCollision::Breed:CR="<<newCR<<",PR="<<newPR<<G4endl;
02022 G4Exception("G4QIonIonCollision::Breeder:","72",FatalException,"2RR-");
02023 }
02024 cLeft->SetPDGCode(newCL);
02025 cRight->SetPDGCode(-newCR);
02026 Left->SetPDGCode(-newPL);
02027 Right->SetPDGCode(newPR);
02028 break;
02029 }
02030 else if(sPDG<0 && uPDG>0)
02031 {
02032 std::pair<G4int,G4int> resLL=ReducePair((-sPDG)/100, uPDG/100);
02033 G4int newCL=resLL.first;
02034 G4int newPL=resLL.second;
02035 if(!newCL || !newPL)
02036 {
02037 G4cerr<<"*G4QIonIonCollision::Breed:CL="<<newCL<<",PL="<<newPL<<G4endl;
02038 G4Exception("G4QIonIonCollision::Breeder:","72",FatalException,"2LL+");
02039 }
02040 std::pair<G4int,G4int> resRR=ReducePair(nPDG/100, (-mPDG)/100);
02041 G4int newCR=resRR.first;
02042 G4int newPR=resRR.second;
02043 if(!newCR || !newPR)
02044 {
02045 G4cerr<<"*G4QIonIonCollision::Breed:CR="<<newCR<<",PR="<<newPR<<G4endl;
02046 G4Exception("G4QIonIonCollision::Breeder:","72",FatalException,"2RR+");
02047 }
02048 cLeft->SetPDGCode(-newCL);
02049 cRight->SetPDGCode(newCR);
02050 Left->SetPDGCode(newPL);
02051 Right->SetPDGCode(-newPR);
02052 break;
02053 }
02054 else if(sPDG>0 && mPDG<0)
02055 {
02056 std::pair<G4int,G4int> resLL=ReducePair(sPDG/100, (-mPDG)/100);
02057 G4int newCL=resLL.first;
02058 G4int newPR=resLL.second;
02059 if(!newCL || !newPR)
02060 {
02061 G4cerr<<"*G4QIonIonCollision::Breed:CL="<<newCL<<",PR="<<newPR<<G4endl;
02062 G4Exception("G4QIonIonCollision::Breeder:","72",FatalException,"2-LR");
02063 }
02064 std::pair<G4int,G4int> resRR=ReducePair((-nPDG)/100, uPDG/100);
02065 G4int newCR=resRR.first;
02066 G4int newPL=resRR.second;
02067 if(!newCR || !newPR)
02068 {
02069 G4cerr<<"*G4QIonIonCollision::Breed:CR="<<newCR<<",PL="<<newPL<<G4endl;
02070 G4Exception("G4QIonIonCollision::Breeder:","72",FatalException,"2-LR");
02071 }
02072 cLeft->SetPDGCode(newCL);
02073 cRight->SetPDGCode(-newCR);
02074 Left->SetPDGCode(newPL);
02075 Right->SetPDGCode(-newPR);
02076 break;
02077 }
02078 else
02079 {
02080 std::pair<G4int,G4int> resLL=ReducePair((-sPDG)/100, mPDG/100);
02081 G4int newCL=resLL.first;
02082 G4int newPR=resLL.second;
02083 if(!newCL || !newPR)
02084 {
02085 G4cerr<<"*G4QIonIonCollision::Breed:CL="<<newCL<<",PR="<<newPR<<G4endl;
02086 G4Exception("G4QIonIonCollision::Breeder:","72",FatalException,"2-RL");
02087 }
02088 std::pair<G4int,G4int> resRR=ReducePair(nPDG/100, (-uPDG)/100);
02089 G4int newCR=resRR.first;
02090 G4int newPL=resRR.second;
02091 if(!newCR || !newPR)
02092 {
02093 G4cerr<<"*G4QIonIonCollision::Breed:CR="<<newCR<<",PL="<<newPL<<G4endl;
02094 G4Exception("G4QIonIonCollision::Breeder:","72",FatalException,"2-RL");
02095 }
02096 cLeft->SetPDGCode(-newCL);
02097 cRight->SetPDGCode(newCR);
02098 Left->SetPDGCode(-newPL);
02099 Right->SetPDGCode(newPR);
02100 break;
02101 }
02102 }
02103 }
02104
02105
02106 #ifdef debug
02107 G4cout<<"G4QInel::Breed: TryFuse/w #"<<restr<<",q="<<rPDG<<",a="<<aPDG<<G4endl;
02108 #endif
02109 G4int PLS=PLT+PRT;
02110 if( (LS==2 && PLS==2) ||
02111 ( ( (LS==2 && PLS==3) || (LS==3 && PLS==2) ) &&
02112 ( (aPDG> 7 && (-dPDG==aPDG/10 || -dPDG==aPDG%10) ) ||
02113 (dPDG> 7 && (-aPDG==dPDG/10 || -aPDG==dPDG%10) ) ||
02114 (rPDG<-7 && (qPDG==(-rPDG)/10 || qPDG==(-rPDG)%10) ) ||
02115 (qPDG<-7 && (rPDG==(-qPDG)/10 || rPDG==(-qPDG)%10) )
02116
02117 )
02118 ) ||
02119 ( ( LS==2 && PLS==4 &&
02120 (aPDG> 7 && (-dPDG == aPDG/10 || -dPDG == aPDG%10) ) &&
02121 (rPDG<-7 && (qPDG==(-rPDG)/10 || qPDG==(-rPDG)%10) )
02122 ) ||
02123 ( LS==4 && PLS==2 &&
02124 (dPDG> 7 && (-aPDG == dPDG/10 || -aPDG == dPDG%10) ) &&
02125 (qPDG<-7 && (rPDG==(-qPDG)/10 || rPDG==(-qPDG)%10) )
02126 )
02127 ) ||
02128 ( LS==3 && PLS==3 &&
02129 ( (aPDG> 7 && (-dPDG == aPDG/10 || -dPDG == aPDG%10) &&
02130 qPDG<-7 && (rPDG==(-qPDG)/10 || rPDG==(-qPDG)%10)
02131 ) ||
02132 (dPDG> 7 && (-aPDG == dPDG/10 || -aPDG == dPDG%10) &&
02133 rPDG<-7 && (dPDG==(-rPDG)/10 || dPDG==(-rPDG)%10)
02134 )
02135 )
02136 )
02137 )
02138 {
02139 G4LorentzVector reString4M = reString->Get4Momentum();
02140 G4ThreeVector reV = reString4M.vect()/reString4M.e();
02141 G4double dV=(curV-reV).mag2();
02142 #ifdef debug
02143 G4cout<<"G4QIonIonCollision::Breeder: StringCand#"<<restr<<", q="<<rPDG
02144 <<", a="<<aPDG<<", L="<<uPDG<<", R="<<mPDG<<",dV="<<dV<<G4endl;
02145 #endif
02146 if(dV < Vmin)
02147 {
02148 Vmin=dV;
02149 fustr=restr;
02150 MPS=PLS;
02151 }
02152 }
02153 }
02154 if(reduce==2)
02155 {
02156 #ifdef debug
02157 G4cout<<"-G4QIonIonCollision::Breed:Reduced #"<<astring<<" & #"<<restr<<G4endl;
02158 #endif
02159 astring--;
02160 continue;
02161 }
02162 if(fustr)
02163 {
02164 #ifdef debug
02165 G4cout<<"G4QInel::Breeder: StPartner#"<<fustr<<", LT="<<LT<<",RT="<<RT<<G4endl;
02166 #endif
02167 G4QString* fuString=strings[fustr];
02168 G4QParton* Left=fuString->GetLeftParton();
02169 G4QParton* Right=fuString->GetRightParton();
02170 G4int uPDG=Left->GetPDGCode();
02171 G4int mPDG=Right->GetPDGCode();
02172 G4int rPDG=uPDG;
02173 G4int aPDG=mPDG;
02174 if(aPDG<-99 || (aPDG>0 && aPDG<7) || rPDG>99 || (rPDG<0 && rPDG>-7))
02175 {
02176 swap=rPDG;
02177 rPDG=aPDG;
02178 aPDG=swap;
02179 }
02180 if(aPDG > 99) aPDG/=100;
02181 if(rPDG <-99) rPDG=-(-rPDG)/100;
02182
02183 G4LorentzVector resL4M(0.,0.,0.,0.);
02184 G4LorentzVector resR4M(0.,0.,0.,0.);
02185 G4LorentzVector L4M=cLeft->Get4Momentum();
02186 G4LorentzVector R4M=cRight->Get4Momentum();
02187 #ifdef debug
02188 G4cout<<"G4QIonIonCollision::Breeder:BeforeFuDir,sL="<<sPDG<<",nR="<<nPDG
02189 <<",uL="<<uPDG<<",mR="<<mPDG<<",L4M="<<L4M<<",R4M="<<R4M<<G4endl;
02190 #endif
02191 G4LorentzVector PL4M=Left->Get4Momentum();
02192 G4LorentzVector PR4M=Right->Get4Momentum();
02193 fusionDONE=AnnihilationOrder(LS, MPS, uPDG, mPDG, sPDG, nPDG);
02194 if (fusionDONE>0)
02195 {
02196 if(fusionDONE>1)
02197 {
02198 if ( (uPDG<0 || nPDG<0) && -uPDG==nPDG ) Left->SetPDGCode(sPDG);
02199 else if( (mPDG<0 || sPDG<0) && -mPDG==sPDG ) Right->SetPDGCode(nPDG);
02200 else if( (uPDG<0 || sPDG<0) && -uPDG==sPDG ) Left->SetPDGCode(nPDG);
02201 else if( (mPDG<0 || nPDG<0) && -mPDG==nPDG ) Right->SetPDGCode(sPDG);
02202 }
02203 {
02204 Left->SetPDGCode(SumPartonPDG(uPDG, sPDG));
02205 Right->SetPDGCode(SumPartonPDG(mPDG, nPDG));
02206 }
02207 Left->Set4Momentum(L4M+PL4M);
02208 Right->Set4Momentum(R4M+PR4M);
02209 #ifdef debug
02210 G4cout<<"G4QIonIonCollision::Breeder:LL/RR s4M="<<fuString->Get4Momentum()
02211 <<",S="<<L4M+PL4M+R4M+PR4M<<", L="<<Left->Get4Momentum()<<", R="
02212 <<Right->Get4Momentum()<<G4endl;
02213 #endif
02214 }
02215 else if(fusionDONE<0)
02216 {
02217 Left->SetPDGCode(SumPartonPDG(uPDG, nPDG));
02218 Left->Set4Momentum(L4M+PR4M);
02219 Right->SetPDGCode(SumPartonPDG(mPDG, sPDG));
02220 Right->Set4Momentum(R4M+PL4M);
02221 #ifdef debug
02222 G4cout<<"G4QIonIonCollision::Breeder:LR/RL s4M="<<fuString->Get4Momentum()
02223 <<",S="<<L4M+PL4M+R4M+PR4M<<", L="<<Left->Get4Momentum()<<", R="
02224 <<Right->Get4Momentum()<<G4endl;
02225 #endif
02226 }
02227 #ifdef debug
02228 else G4cout<<"-Warning-G4QIonIonCollision::Breeder: WrongStringFusion"<<G4endl;
02229 #endif
02230 #ifdef edebug
02231 G4cout<<"#EMC#G4QIonIonCollision::Breed:StringFused,F="<<fusionDONE<<",L="<<L4M
02232 <<",R="<<R4M<<",pL="<<PL4M<<",pR="<<PR4M<<",nL="<<Left->Get4Momentum()
02233 <<",nR="<<Right->Get4Momentum()<<",S="<<fuString->Get4Momentum()<<G4endl;
02234 #endif
02235 if(fusionDONE)
02236 {
02237 #ifdef debug
02238 G4cout<<"###G4QIonIonCollision::Breed: Str#"<<astring<<" fused/w Str#"<<fustr
02239 <<"->"<<fuString->Get4Momentum()<<fuString->Get4Momentum().m()<<G4endl;
02240 #endif
02241 continue;
02242 }
02243 }
02244 #ifdef debug
02245 else
02246 {
02247
02248 G4cerr<<"**G4QIonIonCollision::Breed:*NoPart*M="<<curString->Get4Momentum().m()
02249 <<", F="<<fusionDONE<<", LPDG="<<curString->GetLeftParton()->GetPDGCode()
02250 <<", RPDG="<<curString->GetRightParton()->GetPDGCode()<<G4endl;
02251 }
02252 #endif
02253 }
02254 if(!fusionDONE)
02255 {
02256 G4int nHadr=theResult->size();
02257 #ifdef debug
02258 G4cout<<"+++4QFragmentation::Breeder:*TryHadr* M="<<curString->Get4Momentum().m()
02259 <<", nH="<<nHadr<<", LPDG="<<curString->GetLeftParton()->GetPDGCode()
02260 <<", RPDG="<<curString->GetRightParton()->GetPDGCode()<<G4endl;
02261 #endif
02262
02263 while( (nHadr=theResult->size()) && !theHadrons)
02264 {
02265 #ifdef edebug
02266 for(G4int i=0; i<nHadr; i++)
02267 {
02268 G4LorentzVector h4M=(*theResult)[i]->Get4Momentum();
02269 G4int hPDG=(*theResult)[i]->GetPDGCode();
02270 G4QContent hQC=(*theResult)[i]->GetQC();
02271 G4cout<<"-EMC-G4QInel::Breed:H#"<<i<<",4M="<<h4M<<hQC<<",PDG="<<hPDG<<G4endl;
02272 }
02273 #endif
02274 G4int fusDONE=0;
02275 G4int fuhad=-1;
02276 G4int newPDG=0;
02277 G4int secPDG=0;
02278 G4double maM2=-DBL_MAX;
02279 G4LorentzVector selH4M(0.,0.,0.,0.);
02280 G4QHadron* selHP=0;
02281 G4QString* cString=strings[astring];
02282 G4LorentzVector cString4M = cString->Get4Momentum();
02283 cLeft=cString->GetLeftParton();
02284 cRight=cString->GetRightParton();
02285 G4int sumT=cLeft->GetType()+cRight->GetType();
02286 sPDG=cLeft->GetPDGCode();
02287 nPDG=cRight->GetPDGCode();
02288 G4int cMax=0;
02289 for (G4int reh=0; reh < nHadr; reh++)
02290 {
02291 G4QHadron* curHadr=(*theResult)[reh];
02292 G4QContent curQC=curHadr->GetQC();
02293 G4int partPDG1=curQC.AddParton(sPDG);
02294 G4int partPDG2=curQC.AddParton(nPDG);
02295 #ifdef debug
02296 G4cout<<"G4QIonIonCollision::Breeder: Hadron # "<<reh<<", QC="<<curQC
02297 <<", P1="<<partPDG1<<", P2="<<partPDG2<<G4endl;
02298 #endif
02299 if(partPDG1 || partPDG2)
02300 {
02301 G4int cCur=1;
02302 if(sumT>3 && partPDG1 && partPDG2) cCur=2;
02303 G4LorentzVector curHadr4M = curHadr->Get4Momentum();
02304 G4double M2=(cString4M+curHadr4M).m2();
02305 #ifdef debug
02306 G4cout<<"G4QIonIonCollision::Breeder:*IN*Hadron#"<<reh<<",M2="<<M2<<G4endl;
02307 #endif
02308 if( (sumT<4 || cCur>=cMax) && M2 > maM2)
02309 {
02310 maM2=M2;
02311 fuhad=reh;
02312 if(partPDG1)
02313 {
02314 fusDONE= 1;
02315 newPDG=partPDG1;
02316 if(partPDG2)
02317 {
02318 secPDG=partPDG2;
02319 cMax=cCur;
02320 }
02321 }
02322 else
02323 {
02324 fusDONE=-1;
02325 newPDG=partPDG2;
02326 }
02327 #ifdef debug
02328 G4cout<<"G4QInel::Br:*Selected*,P1="<<partPDG1<<",P2="<<partPDG2<<G4endl;
02329 #endif
02330 selH4M=curHadr4M;
02331 selHP=curHadr;
02332 }
02333 }
02334 }
02335 #ifdef debug
02336 G4cout<<"G4QIonIonCollision::Breeder: fuh="<<fuhad<<",fus="<<fusDONE<<G4endl;
02337 #endif
02338 if(fuhad>-1)
02339 {
02340 if (fusDONE>0)
02341 {
02342 cLeft->SetPDGCode(newPDG);
02343 G4LorentzVector newLeft=cLeft->Get4Momentum()+selH4M;
02344 cLeft->Set4Momentum(newLeft);
02345 if(secPDG && cMax>1)
02346 {
02347 #ifdef debug
02348 G4cout<<"G4QInel::Br:TryReduce, nPDG="<<newPDG<<",sPDG="<<secPDG<<G4endl;
02349 #endif
02350 cLeft->ReduceDiQADiQ(cLeft, cRight);
02351 }
02352 #ifdef debug
02353 G4cout<<"G4QIonIonCollision::Breed: Left, s4M="<<curString->Get4Momentum()
02354 <<", L4M="<<newLeft<<", R4M="<<cRight->Get4Momentum()<<", h4M="
02355 <<selH4M<<", newL="<<newPDG<<", oldR="<<cRight->GetPDGCode()<<G4endl;
02356 #endif
02357 }
02358 else if(fusDONE<0)
02359 {
02360 cRight->SetPDGCode(newPDG);
02361 G4LorentzVector newRight=cRight->Get4Momentum()+selH4M;
02362 cRight->Set4Momentum(newRight);
02363 #ifdef debug
02364 G4cout<<"G4QIonIonCollision::Breed: Right, s4M="<<curString->Get4Momentum()
02365 <<", L4M="<<cLeft->Get4Momentum()<<", R4M="<<newRight<<", h4M="
02366 <<selH4M<<", newR="<<newPDG<<", oldL="<<cLeft->GetPDGCode()<<G4endl;
02367 #endif
02368 }
02369 #ifdef debug
02370 else G4cout<<"-G4QIonIonCollision::Breed: Wrong String+HadronFusion"<<G4endl;
02371 #endif
02372 #ifdef debug
02373 if(fusDONE) G4cout<<"####G4QIonIonCollision::Breeder: String #"<<astring
02374 <<" is fused with Hadron #"<<fuhad
02375 <<", new4Mom="<<curString->Get4Momentum()
02376 <<", M2="<<curString->Get4Momentum().m2()
02377 <<", QC="<<curString->GetQC()<<G4endl;
02378 #endif
02379 }
02380 else
02381 {
02382 #ifdef debug
02383 G4cerr<<"**G4QIonIonCollision::Breed:*NoH*,M="<<curString->Get4Momentum().m()
02384 <<", LPDG="<<curString->GetLeftParton()->GetPDGCode()
02385 <<", RPDG="<<curString->GetRightParton()->GetPDGCode()<<G4endl;
02386
02387
02388 #endif
02389 break;
02390 }
02391
02392 #ifdef debug
02393 G4cout<<"G4QIonIonCollision::Breed: before HR, nH="<<theResult->size()<<G4endl;
02394 G4int icon=0;
02395 #endif
02396 G4QHadronVector::iterator ih;
02397 #ifdef debug
02398 G4bool found=false;
02399 #endif
02400 for(ih = theResult->begin(); ih != theResult->end(); ih++)
02401 {
02402 #ifdef debug
02403 G4cout<<"G4QInelast::Breeder:#"<<icon<<", i="<<(*ih)<<", sH="<<selHP<<G4endl;
02404 icon++;
02405 #endif
02406 if((*ih)==selHP)
02407 {
02408 #ifdef debug
02409 G4cout<<"G4QInel::Breed: *HadronFound*, PDG="<<selHP->GetPDGCode()<<G4endl;
02410 #endif
02411 G4LorentzVector p4M=selHP->Get4Momentum();
02412 curString4M+=p4M;
02413 #ifdef edebug
02414 G4int Chg=selHP->GetCharge();
02415 G4int BaN=selHP->GetBaryonNumber();
02416 curStrChg+=Chg;
02417 curStrBaN+=BaN;
02418 G4cout<<"-EMC->>>G4QIonIonCollision::Breed: S+=H, 4M="<<curString4M<<", M="
02419 <<curString4M.m()<<", Charge="<<curStrChg<<", B="<<curStrBaN<<G4endl;
02420 #endif
02421 delete selHP;
02422 theResult->erase(ih);
02423 #ifdef debug
02424 found=true;
02425 #endif
02426 break;
02427 }
02428 }
02429 #ifdef debug
02430 if(!found) G4cout<<"*G4QIonIonCollision::Breed:nH="<<theResult->size()<<G4endl;
02431 #endif
02432
02433 theHadrons=curString->FragmentString(true);
02434 #ifdef debug
02435 G4cout<<"G4QInel::Breeder: tH="<<theHadrons<<",nH="<<theResult->size()<<G4endl;
02436 #endif
02437 }
02438 #ifdef debug
02439 G4cout<<"*G4QIonIonCollision::Breed: *CanTryToDecay?* nH="<<theHadrons<<", next="
02440 <<next<<" =? nS="<<strings.size()<<", nR="<<theResult->size()<<G4endl;
02441 #endif
02442 if(!theHadrons && next == strings.size() && !(theResult->size()))
02443 {
02444 G4QContent miQC=curString->GetQC();
02445 G4int miPDG=miQC.GetSPDGCode();
02446 if(miPDG == 10)
02447 {
02448 G4QChipolino QCh(miQC);
02449 G4QPDGCode h1QPDG=QCh.GetQPDG1();
02450 G4double h1M =h1QPDG.GetMass();
02451 G4QPDGCode h2QPDG=QCh.GetQPDG2();
02452 G4double h2M =h2QPDG.GetMass();
02453 G4double ttM =curString4M.m();
02454 if(h1M+h2M<ttM+eps)
02455 {
02456 G4LorentzVector h14M(0.,0.,0.,h1M);
02457 G4LorentzVector h24M(0.,0.,0.,h2M);
02458 if(std::fabs(ttM-h1M-h2M)<=eps)
02459 {
02460 G4double part1=h1M/(h1M+h2M);
02461 h14M=part1*curString4M;
02462 h24M=curString4M-h14M;
02463 }
02464 else
02465 {
02466 if(!G4QHadron(curString4M).DecayIn2(h14M,h24M))
02467 {
02468 G4cerr<<"***G4QIonIonCollision::Breeder: tM="<<ttM<<"->h1="<<h1QPDG
02469 <<"("<<h1M<<")+h2="<<h1QPDG<<"("<<h2M<<")="<<h1M+h2M<<G4endl;
02470 G4Exception("G4QIonIonCollision::Breeder:","72",FatalException,"QDec");
02471 }
02472 }
02473 G4QHadron* h1H = new G4QHadron(h1QPDG.GetPDGCode(),h14M);
02474 theResult->push_back(h1H);
02475 #ifdef debug
02476 G4LorentzVector f4M=h1H->Get4Momentum();
02477 G4int fPD=h1H->GetPDGCode();
02478 G4int fCg=h1H->GetCharge();
02479 G4int fBN=h1H->GetBaryonNumber();
02480 G4cout<<"-EMC->>G4QIonIonCollision::Breed:String=HadrChiPro1's filled,f4M="
02481 <<f4M<<", fPDG="<<fPD<<", fCg="<<fCg<<", fBN="<<fBN<<G4endl;
02482 #endif
02483 G4QHadron* h2H = new G4QHadron(h2QPDG.GetPDGCode(),h24M);
02484 theResult->push_back(h2H);
02485 #ifdef debug
02486 G4LorentzVector s4M=h2H->Get4Momentum();
02487 G4int sPD=h2H->GetPDGCode();
02488 G4int sCg=h2H->GetCharge();
02489 G4int sBN=h2H->GetBaryonNumber();
02490 G4cout<<"-EMC->>G4QIonIonCollision::Breed:String=HadrChiPro2's filled,s4M="
02491 <<s4M<<", sPDG="<<sPD<<", sCg="<<sCg<<", sBN="<<sBN<<G4endl;
02492 #endif
02493 #ifdef edebug
02494 G4cout<<"-EMC-..Chi..G4QIonIonCollision::Breeder: DecayCHECK, Ch4M="
02495 <<curString4M<<", d4M="<<curString4M-h14M-h24M<<G4endl;
02496 #endif
02497 break;
02498 }
02499 else
02500 {
02501 G4Quasmon* stringQuasmon = new G4Quasmon(miQC, curString4M);
02502 theQuasmons.push_back(stringQuasmon);
02503 break;
02504 }
02505 }
02506 else if(miPDG)
02507 {
02508 if (miPDG>0 && miPDG%10 < 3) miPDG+=2;
02509 else if(miPDG<0 && (-miPDG)%10< 3) miPDG-=2;
02510 G4Quasmon Quasm;
02511 G4QHadron* sHad = new G4QHadron(miPDG,curString4M);
02512 G4QHadronVector* tmpQHadVec=Quasm.DecayQHadron(sHad);
02513 G4int tmpN=tmpQHadVec->size();
02514 #ifdef debug
02515 G4cout<<"G4QIonIonCollision::Breeder: Decay the Last, Res#H="<<tmpN<<G4endl;
02516 #endif
02517 if(tmpN>1)
02518 {
02519 for(G4int aH=0; aH < tmpN; aH++)
02520 {
02521 theResult->push_back((*tmpQHadVec)[aH]);
02522 #ifdef debug
02523 G4QHadron* prodH =(*tmpQHadVec)[aH];
02524 G4LorentzVector p4M=prodH->Get4Momentum();
02525 G4int PDG=prodH->GetPDGCode();
02526 G4int Chg=prodH->GetCharge();
02527 G4int BaN=prodH->GetBaryonNumber();
02528 G4cout<<"-EMC->>G4QIonIonCollis::Breed:String=Hadr,H#"<<aH<<" filled,4M="
02529 <<p4M<<", PDG="<<PDG<<", Chg="<<Chg<<", BaN="<<BaN<<G4endl;
02530 #endif
02531 }
02532 }
02533 else
02534 {
02535 G4Quasmon* stringQuasmon = new G4Quasmon(miQC, curString4M);
02536 #ifdef debug
02537 G4cout<<"G4QIonIonCollision::Breeder:==> to Quasm="<<miQC<<curString4M
02538 <<", pNuc="<<theProjNucleus<<theProjNucleus.Get4Momentum()<<", tNuc="
02539 <<theTargNucleus<<theTargNucleus.Get4Momentum()<<", NString="
02540 <<strings.size()<<", nR="<<theResult->size()<<", nQ="
02541 <<theQuasmons.size()<<G4endl;
02542 #endif
02543 theQuasmons.push_back(stringQuasmon);
02544 delete sHad;
02545 tmpQHadVec->clear();
02546 delete tmpQHadVec;
02547 break;
02548 }
02549 tmpQHadVec->clear();
02550 delete tmpQHadVec;
02551 break;
02552 }
02553 }
02554 }
02555 }
02556
02557 #ifdef debug
02558 G4cout<<"G4QIonIonCollision::Breeder: theH="<<theHadrons<<"?=0, next="<<next<<G4endl;
02559 #endif
02560 if(!theHadrons && next < strings.size())
02561 {
02562
02563 G4QContent miQC=curString->GetQC();
02564 G4int miPDG=miQC.GetSPDGCode();
02565 #ifdef debug
02566 G4cout<<">>>G4QIonIonCollision::Breeder: SQC="<<miQC<<", miSPDG="<<miPDG<<G4endl;
02567 #endif
02568 G4double miM=0.;
02569 if(miPDG!=10) miM=G4QPDGCode(miPDG).GetMass();
02570 else
02571 {
02572 G4QChipolino QCh(miQC);
02573 miM=QCh.GetQPDG1().GetMass()+QCh.GetQPDG2().GetMass();
02574 }
02575 G4double cM2=curString4M.m2();
02576 #ifdef debug
02577 G4cout<<">>>G4QIonIonCollision::Breeder: minMass="<<miM<<", realM2="<<cM2<<G4endl;
02578 #endif
02579 G4double cM=0.;
02580 if(cM2>0.)
02581 {
02582 cM=std::sqrt(cM2);
02583 if(std::fabs(cM-miM) < eps)
02584 {
02585 if(miPDG!=10)
02586 {
02587 G4QHadron* sHad = new G4QHadron(miPDG,curString4M);
02588 theResult->push_back(sHad);
02589 #ifdef debug
02590 G4cout<<">>>G4QIonIonCollision::Breeder:S->H="<<miPDG<<curString4M<<G4endl;
02591 #endif
02592 }
02593 else
02594 {
02595 G4QChipolino QCh(miQC);
02596 G4QPDGCode h1QPDG=QCh.GetQPDG1();
02597 G4double h1M =h1QPDG.GetMass();
02598 G4QPDGCode h2QPDG=QCh.GetQPDG2();
02599 G4double h2M =h2QPDG.GetMass();
02600 G4double pt1 =h1M/(h1M+h2M);
02601 G4LorentzVector h14M=pt1*curString4M;
02602 G4LorentzVector h24M=curString4M-h14M;
02603 G4QHadron* h1H = new G4QHadron(h1QPDG.GetPDGCode(),h14M);
02604 theResult->push_back(h1H);
02605 #ifdef debug
02606 G4LorentzVector f4M=h1H->Get4Momentum();
02607 G4int fPD=h1H->GetPDGCode();
02608 G4int fCg=h1H->GetCharge();
02609 G4int fBN=h1H->GetBaryonNumber();
02610 G4cout<<"-EMC->>G4QIonIonCollision::Breed:Str=2HadrAR Prod-F is filled, f4M="
02611 <<f4M<<", fPDG="<<fPD<<", fCg="<<fCg<<", fBN="<<fBN<<G4endl;
02612 #endif
02613 G4QHadron* h2H = new G4QHadron(h2QPDG.GetPDGCode(),h24M);
02614 theResult->push_back(h2H);
02615 #ifdef debug
02616 G4LorentzVector s4M=h2H->Get4Momentum();
02617 G4int sPD=h2H->GetPDGCode();
02618 G4int sCg=h2H->GetCharge();
02619 G4int sBN=h2H->GetBaryonNumber();
02620 G4cout<<"-EMC->>G4QIonIonCollision::Breed:Str=2HadrAR Prod-S is filled, s4M="
02621 <<s4M<<", sPDG="<<sPD<<", sCg="<<sCg<<", sBN="<<sBN<<G4endl;
02622 #endif
02623 }
02624 continue;
02625 }
02626 else
02627 {
02628 G4ThreeVector cP=curString4M.vect();
02629 G4double cE=curString4M.e();
02630 G4ThreeVector curV=cP/cE;
02631 G4double miM2=miM*miM;
02632 G4int restr=0;
02633 G4int fustr=0;
02634 G4double selX=0.;
02635 G4double maD=-DBL_MAX;
02636 G4double Vmin=DBL_MAX;
02637 G4LorentzVector s4M(0.,0.,0.,0.);
02638 #ifdef debug
02639 G4cout<<"G4QIonIonCollision::Breeder: TryRecover, cV="<<curV<<G4endl;
02640 #endif
02641 nOfStr=strings.size();
02642 for(restr=next; restr < nOfStr; ++restr) if(restr != astring)
02643 {
02644 G4QString* pString=strings[restr];
02645 G4LorentzVector p4M=pString->Get4Momentum();
02646 G4ThreeVector pP=p4M.vect();
02647 G4double pE=p4M.e();
02648 G4double D2=cE*pE-cP.dot(pP);
02649 G4double pM2=p4M.m2();
02650 G4double dM4=pM2*(miM2-cM2);
02651 G4double D=D2*D2+dM4;
02652 G4double x=-1.;
02653 if(D >= 0. && pM2>.01) x=(std::sqrt(D)-D2)/pM2;
02654 #ifdef debug
02655 else G4cout<<"G4QIonIonCollis::Breed:D="<<D<<",D2="<<D2<<",dM="<<dM4<<G4endl;
02656 G4cout<<"G4QIonIonCollision::Breed: pM2="<<pM2<<",D2="<<D2<<",x="<<x<<G4endl;
02657 #endif
02658 if(x > 0. && x < 1.)
02659 {
02660 G4QContent pQC=pString->GetQC();
02661 G4int pPDG=pQC.GetSPDGCode();
02662 G4double pM=0.;
02663 if(pPDG==10)
02664 {
02665 G4QChipolino QCh(pQC);
02666 pM=QCh.GetQPDG1().GetMass()+QCh.GetQPDG2().GetMass();
02667 }
02668 else pM=G4QPDGCode(pPDG).GetMass();
02669 G4double rM=std::sqrt(pM2);
02670 G4double delta=(1.-x)*rM-pM;
02671 if(delta > 0. && delta > maD)
02672 {
02673 maD=delta;
02674 #ifdef debug
02675 G4cout<<"G4QIonIonCollision::Breed: Subtr,S#"<<restr<<",d="<<maD<<G4endl;
02676 #endif
02677 fustr=restr;
02678 selX=x;
02679 s4M=p4M;
02680 }
02681 }
02682 else if(x <= 0.)
02683 {
02684 G4ThreeVector pV=pP/pE;
02685 G4double dV=(curV-pV).mag2();
02686 if(dV < Vmin)
02687 {
02688 #ifdef debug
02689 G4cout<<"G4QIonIonCollision::Breed: FreeAdd,S#"<<restr<<",x="<<x<<G4endl;
02690 #endif
02691 Vmin=dV;
02692 fustr=restr;
02693 selX=x;
02694 s4M=p4M;
02695 }
02696 }
02697 #ifdef debug
02698 G4cout<<"G4QIonIonCollision::Breed:EndOfLOOP r="<<restr<<"<"<<nOfStr<<G4endl;
02699 #endif
02700 }
02701 #ifdef debug
02702 G4cout<<"G4QIonIonCollision::Breeder: AfterLOOP fustr="<<fustr<<G4endl;
02703 #endif
02704 if(fustr)
02705 {
02706 #ifdef edebug
02707 G4LorentzVector sum4M=s4M+curString4M;
02708 G4cout<<"G4QIonIonCollision::Breeder: Found Sum4M="<<sum4M<<G4endl;
02709 #endif
02710 G4QString* pString=strings[fustr];
02711 curString4M+=selX*s4M;
02712 if(std::abs(miPDG)%10 > 2)
02713 {
02714 G4Quasmon Quasm;
02715 G4QHadron* sHad = new G4QHadron(miPDG,curString4M);
02716 G4QHadronVector* tmpQHadVec=Quasm.DecayQHadron(sHad);
02717 #ifdef debug
02718 G4cout<<"G4QIonIonCollision::Breed:DecStH,nH="<<tmpQHadVec->size()<<G4endl;
02719 #endif
02720 for(unsigned aH=0; aH < tmpQHadVec->size(); aH++)
02721 {
02722 theResult->push_back((*tmpQHadVec)[aH]);
02723 #ifdef debug
02724 G4QHadron* prodH =(*tmpQHadVec)[aH];
02725 G4LorentzVector p4M=prodH->Get4Momentum();
02726 G4int PDG=prodH->GetPDGCode();
02727 G4int Chg=prodH->GetCharge();
02728 G4int BaN=prodH->GetBaryonNumber();
02729 G4cout<<"-EMC->>G4QIonIonCollision::Breed:St=Had,pH#"<<aH<<" filled, 4M="
02730 <<p4M<<", PDG="<<PDG<<", Chg="<<Chg<<", BaN="<<BaN<<G4endl;
02731 #endif
02732 }
02733 tmpQHadVec->clear();
02734 delete tmpQHadVec;
02735 }
02736 else if(miPDG == 10)
02737 {
02738 G4QChipolino QCh(miQC);
02739 G4QPDGCode h1QPDG=QCh.GetQPDG1();
02740 G4double h1M =h1QPDG.GetMass();
02741 G4QPDGCode h2QPDG=QCh.GetQPDG2();
02742 G4double h2M =h2QPDG.GetMass();
02743 G4double ttM =curString4M.m();
02744 if(h1M+h2M<ttM+eps)
02745 {
02746 G4LorentzVector h14M(0.,0.,0.,h1M);
02747 G4LorentzVector h24M(0.,0.,0.,h2M);
02748 if(std::fabs(ttM-h1M-h2M)<=eps)
02749 {
02750 G4double part1=h1M/(h1M+h2M);
02751 h14M=part1*curString4M;
02752 h24M=curString4M-h14M;
02753 }
02754 else
02755 {
02756 if(!G4QHadron(curString4M).DecayIn2(h14M,h24M))
02757 {
02758 G4cerr<<"***G4QIonIonCollision::Breeder: tM="<<ttM<<"->h1="<<h1QPDG
02759 <<"("<<h1M<<")+h2="<<h1QPDG<<"("<<h2M<<")="<<h1M+h2M<<G4endl;
02760 G4Exception("G4QIonIonCollision::Breeder:","72",FatalException,"CD");
02761 }
02762 }
02763 G4QHadron* h1H = new G4QHadron(h1QPDG.GetPDGCode(),h14M);
02764 theResult->push_back(h1H);
02765 #ifdef debug
02766 G4LorentzVector f4M=h1H->Get4Momentum();
02767 G4int fPD=h1H->GetPDGCode();
02768 G4int fCg=h1H->GetCharge();
02769 G4int fBN=h1H->GetBaryonNumber();
02770 G4cout<<"-EMC->>>G4QIonIonCollision::Breed:Str=Hadr Prod-F's filled,f4M="
02771 <<f4M<<", fPDG="<<fPD<<", fCg="<<fCg<<", fBN="<<fBN<<G4endl;
02772 #endif
02773 G4QHadron* h2H = new G4QHadron(h2QPDG.GetPDGCode(),h24M);
02774 theResult->push_back(h2H);
02775 #ifdef debug
02776 G4LorentzVector s4M=h2H->Get4Momentum();
02777 G4int sPD=h2H->GetPDGCode();
02778 G4int sCg=h2H->GetCharge();
02779 G4int sBN=h2H->GetBaryonNumber();
02780 G4cout<<"-EMC->>>G4QIonIonCollision::Breed:Str=Hadr Prod-S's filled,s4M="
02781 <<s4M<<", sPDG="<<sPD<<", sCg="<<sCg<<", sBN="<<sBN<<G4endl;
02782 #endif
02783 #ifdef edebug
02784 G4cout<<"-EMC-Chipo.G4QIonIonCollision::Breed:DecCHECK,c4M="<<curString4M
02785 <<", ChQC="<<miQC<<", d4M="<<curString4M-h14M-h24M<<G4endl;
02786 #endif
02787 }
02788 else
02789 {
02790 G4cerr<<"***G4QInel::Breeder: tM="<<ttM<<miQC<<"->h1="<<h1QPDG<<"(" <<h1M
02791 <<")+h2="<<h1QPDG<<"("<<h2M<<") = "<<h1M+h2M<<G4endl;
02792 G4Exception("G4QIonIonCollision::Breeder:","72",FatalException,"ChiDec");
02793 }
02794 }
02795 else
02796 {
02797 G4QHadron* sHad = new G4QHadron(miPDG,curString4M);
02798 theResult->push_back(sHad);
02799 #ifdef debug
02800 G4cout<<"-EMC->>>G4QIonIonCollision::Breeder:Str=Hadr Filled, 4M="
02801 <<curString4M<<", PDG="<<miPDG<<G4endl;
02802 #endif
02803 }
02804 G4double corF=1-selX;
02805 G4QParton* Left=pString->GetLeftParton();
02806 G4QParton* Right=pString->GetRightParton();
02807 Left->Set4Momentum(corF*Left->Get4Momentum());
02808 Right->Set4Momentum(corF*Right->Get4Momentum());
02809 #ifdef edebug
02810 G4cout<<"-EMC-...Cor...G4QIonIonCollision::Breeder:CorCHECK Sum="<<sum4M
02811 <<" =? "<<curString4M+pString->Get4Momentum()<<", M="<<miM<<" ?= "
02812 <<curString4M.m()<<G4endl;
02813 #endif
02814 #ifdef debug
02815 G4cout<<">>>G4QIonIonCollision::Breeder:*Corrected* String->Hadr="<<miPDG
02816 <<curString4M<<" by String #"<<fustr<<G4endl;
02817 #endif
02818 continue;
02819 }
02820 }
02821 }
02822 }
02823 #ifdef debug
02824 else G4cerr<<"***G4QIonIonCollision::Breed:**No SSCorrection**, next="<<next<<G4endl;
02825 #endif
02826
02827 G4QParton* lpcS=curString->GetLeftParton();
02828 G4QParton* rpcS=curString->GetRightParton();
02829 G4int lPDGcS=lpcS->GetPDGCode();
02830 G4int rPDGcS=rpcS->GetPDGCode();
02831 if (lPDGcS==3 && rPDGcS==-3)
02832 {
02833 lpcS->SetPDGCode( 1);
02834 rpcS->SetPDGCode(-1);
02835 }
02836 else if(lPDGcS==-3 && rPDGcS==3)
02837 {
02838 lpcS->SetPDGCode(-1);
02839 rpcS->SetPDGCode( 1);
02840 }
02841
02842 G4int nofRH=theResult->size();
02843 #ifdef debug
02844 G4cout<<"G4QIonIonCollision::Breeder: theH="<<theHadrons<<", #OfH="<<nofRH<<G4endl;
02845 #endif
02846 if(!theHadrons && nofRH)
02847 {
02848 #ifdef debug
02849 G4cerr<<"!G4QIonIonCollision::Breeder:CanTrySHCor, nH="<<theResult->size()<<G4endl;
02850 #endif
02851
02852 G4QContent miQC=curString->GetQC();
02853 G4int miPDG=miQC.GetSPDGCode();
02854 G4double miM=0.;
02855 if(miPDG==10)
02856 {
02857 G4QChipolino QCh(miQC);
02858 miM=QCh.GetQPDG1().GetMass()+QCh.GetQPDG2().GetMass();
02859 }
02860 else miM=G4QPDGCode(miPDG).GetMass();
02861 G4double spM=0.;
02862 G4ThreeVector cP=curString4M.vect();
02863 G4double cE=curString4M.e();
02864 G4ThreeVector curV=cP/cE;
02865 G4int reha=0;
02866 G4int fuha=0;
02867 G4double dMmin=DBL_MAX;
02868 G4LorentzVector s4M(0.,0.,0.,0.);
02869 G4double sM=0.;
02870 for (reha=next; reha < nofRH; reha++)
02871 {
02872 G4QHadron* pHadron=(*theResult)[reha];
02873 G4LorentzVector p4M=pHadron->Get4Momentum();
02874 G4double pM=p4M.m();
02875 G4LorentzVector t4M=p4M+curString4M;
02876 G4double tM2=t4M.m2();
02877 if(tM2 >= sqr(pM+miM+eps))
02878 {
02879 G4double tM=std::sqrt(tM2);
02880 G4double dM=tM-pM-miM;
02881 if(dM < dMmin)
02882 {
02883 dMmin=dM;
02884 fuha=reha;
02885 spM=pM;
02886 s4M=t4M;
02887 sM=tM;
02888 }
02889 }
02890 }
02891 if(fuha)
02892 {
02893 G4QHadron* pHadron=(*theResult)[fuha];
02894 G4LorentzVector mi4M(0.,0.,0.,miM);
02895 if(miM+spM<sM+eps)
02896 {
02897 G4LorentzVector sp4M(0.,0.,0.,spM);
02898 if(std::fabs(sM-miM-spM)<=eps)
02899 {
02900 G4double part1=miM/(miM+spM);
02901 mi4M=part1*s4M;
02902 sp4M=s4M-mi4M;
02903 }
02904 else
02905 {
02906 if(!G4QHadron(s4M).DecayIn2(mi4M,sp4M))
02907 {
02908 G4cerr<<"***G4QIonIonCollision::Breeder: *SH*, tM="<<sM<<"->h1=("<<miPDG
02909 <<")"<<miM<<" + h2="<<spM<<" = "<<miM+spM<<G4endl;
02910 G4Exception("G4QIonIonCollision::Breeder:","72",FatalException,"SHChiDec");
02911 }
02912 }
02913 pHadron->Set4Momentum(sp4M);
02914 #ifdef debug
02915 G4cout<<"-EMC->...>G4QIonIonCollision::Breed: H# "<<fuha<<" is updated, new4M="
02916 <<sp4M<<G4endl;
02917 #endif
02918 }
02919 else
02920 {
02921 G4cerr<<"***G4QInel::Breeder: HS Failed, tM="<<sM<<"->h1M=("<<miPDG<<")"<<miM
02922 <<"+h2M="<<spM<<" = "<<miM+spM<<G4endl;
02923 G4Exception("G4QIonIonCollision::Breeder:","72",FatalException,"HSChipoliDec");
02924 }
02925 if(std::abs(miPDG)%10 > 2)
02926 {
02927 G4Quasmon Quasm;
02928 G4QHadron* sHad = new G4QHadron(miPDG,mi4M);
02929 G4QHadronVector* tmpQHadVec=Quasm.DecayQHadron(sHad);
02930 #ifdef debug
02931 G4cout<<"G4QInelast::Breeder: *HS* DecStrHad, nH="<<tmpQHadVec->size()<<G4endl;
02932 #endif
02933 for(unsigned aH=0; aH < tmpQHadVec->size(); aH++)
02934 {
02935 theResult->push_back((*tmpQHadVec)[aH]);
02936 #ifdef debug
02937 G4QHadron* prodH =(*tmpQHadVec)[aH];
02938 G4LorentzVector p4M=prodH->Get4Momentum();
02939 G4int PDG=prodH->GetPDGCode();
02940 G4int Chg=prodH->GetCharge();
02941 G4int BaN=prodH->GetBaryonNumber();
02942 G4cout<<"-EMC->>>G4QIonIonCollision::Breed:Str+Hadr PrH#"<<aH<<" filled, 4M="
02943 <<p4M<<", PDG="<<PDG<<", Chg="<<Chg<<", BaN="<<BaN<<G4endl;
02944 #endif
02945 }
02946 tmpQHadVec->clear();
02947 delete tmpQHadVec;
02948 }
02949 else if(miPDG == 10)
02950 {
02951 G4QChipolino QCh(miQC);
02952 G4QPDGCode h1QPDG=QCh.GetQPDG1();
02953 G4double h1M =h1QPDG.GetMass();
02954 G4QPDGCode h2QPDG=QCh.GetQPDG2();
02955 G4double h2M =h2QPDG.GetMass();
02956 G4double ttM =curString4M.m();
02957 if(h1M+h2M<miM+eps)
02958 {
02959 G4LorentzVector h14M(0.,0.,0.,h1M);
02960 G4LorentzVector h24M(0.,0.,0.,h2M);
02961 if(std::fabs(ttM-h1M-h2M)<=eps)
02962 {
02963 G4double part1=h1M/(h1M+h2M);
02964 h14M=part1*mi4M;
02965 h24M=mi4M-h14M;
02966 }
02967 else
02968 {
02969 if(!G4QHadron(mi4M).DecayIn2(h14M,h24M))
02970 {
02971 G4cerr<<"***G4QIonIonCollision::Breeder: HS tM="<<ttM<<"->h1="<<h1QPDG
02972 <<"("<<h1M<<")+h2="<<h1QPDG<<"("<<h2M<<")="<<h1M+h2M<<G4endl;
02973 G4Exception("G4QIonIonCollision::Breeder:","72",FatalException,"ChiDec");
02974 }
02975 }
02976 G4QHadron* h1H = new G4QHadron(h1QPDG.GetPDGCode(),h14M);
02977 theResult->push_back(h1H);
02978 #ifdef debug
02979 G4LorentzVector f4M=h1H->Get4Momentum();
02980 G4int fPD=h1H->GetPDGCode();
02981 G4int fCg=h1H->GetCharge();
02982 G4int fBN=h1H->GetBaryonNumber();
02983 G4cout<<"-EMC->>G4QIonIonCollision::Breed: CorStrHadr Prod-1 is filled, f4M="
02984 <<f4M<<", fPDG="<<fPD<<", fCg="<<fCg<<", fBN="<<fBN<<G4endl;
02985 #endif
02986 G4QHadron* h2H = new G4QHadron(h2QPDG.GetPDGCode(),h24M);
02987 theResult->push_back(h2H);
02988 #ifdef debug
02989 G4LorentzVector n4M=h2H->Get4Momentum();
02990 G4int nPD=h2H->GetPDGCode();
02991 G4int nCg=h2H->GetCharge();
02992 G4int nBN=h2H->GetBaryonNumber();
02993 G4cout<<"-EMC->>>G4QIonIonCollision::Breed:CorStrHadr Prod-2 is filled, n4M="
02994 <<n4M<<", nPDG="<<nPD<<", nCg="<<nCg<<", nBN="<<nBN<<G4endl;
02995 #endif
02996 #ifdef edebug
02997 G4cout<<"-EMC-...HS-Chipo...G4QIonIonCollision::Breeder:DecCHECK, Ch4M="
02998 <<mi4M<<", ChQC="<<miQC<<", d4M="<<mi4M-h14M-h24M<<G4endl;
02999 #endif
03000 }
03001 }
03002 else
03003 {
03004 G4QHadron* sHad = new G4QHadron(miPDG, mi4M);
03005 theResult->push_back(sHad);
03006 #ifdef debug
03007 G4cout<<">>..>>G4QIonIonCollision::Breeder: CorStr=Hadr is Filled, 4M="
03008 <<curString4M<<", StPDG="<<miPDG<<G4endl;
03009 #endif
03010 }
03011 #ifdef edebug
03012 G4cout<<"-EMC-...Cor...G4QIonIonCollision::Breeder:StHadCor CHECK Sum="<<s4M
03013 <<" =? "<<mi4M+pHadron->Get4Momentum()<<G4endl;
03014 #endif
03015 #ifdef debug
03016 G4cout<<">>>G4QIonIonCollision::Breeder:*Corrected* String+Hadr="<<miPDG
03017 <<mi4M<<" by Hadron #"<<reha<<G4endl;
03018 #endif
03019 continue;
03020 }
03021 else
03022 {
03023 #ifdef debug
03024 G4cout<<"-EMC->>>G4QIonIonCollision::Breeder: Str+Hadr Failed, 4M="<<curString4M
03025 <<", PDG="<<miPDG<<G4endl;
03026 #endif
03027 }
03028
03029 G4QContent curStringQC=curString->GetQC();
03030 G4Quasmon* stringQuasmon = new G4Quasmon(curStringQC, curString4M);
03031 theQuasmons.push_back(stringQuasmon);
03032 continue;
03033 }
03034 }
03035 G4Quasmon tmpQ;
03036 G4int nHfin=0;
03037 if(theHadrons) nHfin=theHadrons->size();
03038 else
03039 {
03040 G4LorentzVector ss4M(0.,0.,0.,0.);
03041 G4QContent ssQC(0,0,0,0,0,0);
03042 G4int tnSt=strings.size();
03043 for(G4int i=astring; i < tnSt; ++i)
03044 {
03045 G4LorentzVector pS4M=strings[i]->Get4Momentum();
03046 ss4M+=pS4M;
03047 G4QContent pSQC=strings[i]->GetQC();
03048 ssQC+=pSQC;
03049 #ifdef debug
03050 G4cout<<"=--=>G4QIonIonCollision::Breed:S#"<<i<<",4M="<<pS4M<<",QC="<<pSQC<<G4endl;
03051 #endif
03052 }
03053 #ifdef debug
03054 G4cout<<"==>G4QIonIonCollision::Breed:AllStrings are summed up in a Quasmon"<<G4endl;
03055 #endif
03056 G4Quasmon* stringQuasmon = new G4Quasmon(ssQC, ss4M);
03057 theQuasmons.push_back(stringQuasmon);
03058 break;
03059 }
03060 #ifdef debug
03061 G4cout<<"G4QIonIonCollision::Breeder: Trying to decay hadrons #ofHRes="<<nHfin<<G4endl;
03062 #endif
03063 for(G4int aTrack=0; aTrack<nHfin; aTrack++)
03064 {
03065 G4QHadron* curHadron=(*theHadrons)[aTrack];
03066 G4int hPDG=curHadron->GetPDGCode();
03067 #ifdef edebug
03068 G4LorentzVector curH4M=curHadron->Get4Momentum();
03069 G4int curHCh=curHadron->GetCharge();
03070 G4int curHBN=curHadron->GetBaryonNumber();
03071 #endif
03072 #ifdef debug
03073 G4cout<<">>...>>G4QIonIonCollision::Breeder:S#"<<astring<<",H#"<<aTrack<<",PDG="
03074 <<hPDG<<",4M="<<curHadron->Get4Momentum()<<G4endl;
03075 #endif
03076 if(std::abs(hPDG)%10 > 2)
03077 {
03078 G4QHadronVector* tmpQHadVec=tmpQ.DecayQHadron(curHadron);
03079 #ifdef debug
03080 G4cout<<"G4QIonIonCollision::Breed:-DECAY'S DONE-,nH="<<tmpQHadVec->size()<<G4endl;
03081 #endif
03082
03083
03084
03085 for(unsigned aH=0; aH < tmpQHadVec->size(); aH++)
03086 {
03087 theResult->push_back((*tmpQHadVec)[aH]);
03088 #ifdef edebug
03089 G4QHadron* prodH =(*tmpQHadVec)[aH];
03090 G4LorentzVector p4M=prodH->Get4Momentum();
03091 G4int PDG=prodH->GetPDGCode();
03092 G4int Chg=prodH->GetCharge();
03093 G4int BaN=prodH->GetBaryonNumber();
03094 curH4M-=p4M;
03095 curString4M-=p4M;
03096 curStrChg-=Chg;
03097 curStrBaN-=BaN;
03098 curHCh-=Chg;
03099 curHBN-=BaN;
03100 G4cout<<"-EMC->.>>G4QIonIonCollision::Breed:Str*Filled, 4M="<<p4M<<", PDG="<<PDG
03101 <<", Chg="<<Chg<<", BaN="<<BaN<<G4endl;
03102 #endif
03103 }
03104 #ifdef edebug
03105 G4cout<<"-EMC-.G4QIn::Br:Dec,r4M="<<curH4M<<",rC="<<curHCh<<",rB="<<curHBN<<G4endl;
03106 #endif
03107 tmpQHadVec->clear();
03108 delete tmpQHadVec;
03109 }
03110 else
03111 {
03112 theResult->push_back(curHadron);
03113 #ifdef edebug
03114 curString4M-=curH4M;
03115 G4int curCh=curHadron->GetCharge();
03116 G4int curBN=curHadron->GetBaryonNumber();
03117 curStrChg-=curCh;
03118 curStrBaN-=curBN;
03119 G4cout<<"-EMC->>..>>G4QIonIonCollision::Breeder: curH filled 4M="<<curH4M<<",PDG="
03120 <<curHadron->GetPDGCode()<<", Chg="<<curCh<<", BaN="<<curBN<<G4endl;
03121 #endif
03122 }
03123 }
03124
03125 if(theHadrons) delete theHadrons;
03126 #ifdef edebug
03127 G4cout<<"-EMC-.......G4QIonIonCollision::Breeder: StringDecay CHECK, r4M="<<curString4M
03128 <<", rChg="<<curStrChg<<", rBaN="<<curStrBaN<<G4endl;
03129 #endif
03130 }
03131 G4LorentzVector rp4M=theProjNucleus.Get4Momentum();
03132 G4int rpPDG=theProjNucleus.GetPDG();
03133 G4QHadron* resPNuc = new G4QHadron(rpPDG,rp4M);
03134 theResult->push_back(resPNuc);
03135 G4LorentzVector rt4M=theTargNucleus.Get4Momentum();
03136 G4int rtPDG=theTargNucleus.GetPDG();
03137 G4QHadron* resTNuc = new G4QHadron(rtPDG,rt4M);
03138 theResult->push_back(resTNuc);
03139 #ifdef edebug
03140 G4LorentzVector s4M(0.,0.,0.,0.);
03141 G4int rCh=totChg;
03142 G4int rBN=totBaN;
03143 G4int nHadr=theResult->size();
03144 G4int nQuasm=theQuasmons.size();
03145 G4cout<<"-EMCLS-G4QInel::Breeder: #ofHadr="<<nHadr<<", #OfQ="<<nQuasm<<", rPA="<<rp4M.m()
03146 <<"="<<G4QNucleus(rpPDG).GetGSMass()<<", rTA="<<rt4M.m()<<"="
03147 <<G4QNucleus(rtPDG).GetGSMass()<<G4endl;
03148 for(G4int i=0; i<nHadr; i++)
03149 {
03150 G4LorentzVector hI4M=(*theResult)[i]->Get4Momentum();
03151 s4M+=hI4M;
03152 G4int hChg=(*theResult)[i]->GetCharge();
03153 rCh-=hChg;
03154 G4int hBaN=(*theResult)[i]->GetBaryonNumber();
03155 rBN-=hBaN;
03156 G4cout<<"-EMCLS-G4QIonIonCollision::Breeder: Hadron#"<<i<<", 4M="<<hI4M<<", PDG="
03157 <<(*theResult)[i]->GetPDGCode()<<", C="<<hChg<<", B="<<hBaN<<G4endl;
03158 }
03159 for(G4int i=0; i<nQuasm; i++)
03160 {
03161 G4LorentzVector hI4M=theQuasmons[i]->Get4Momentum();
03162 s4M+=hI4M;
03163 G4int hChg=theQuasmons[i]->GetCharge();
03164 rCh-=hChg;
03165 G4int hBaN=theQuasmons[i]->GetBaryonNumber();
03166 rBN-=hBaN;
03167 G4cout<<"-EMCLS-G4QIonIonCollision::Breeder: Quasmon#"<<i<<", 4M="<<hI4M<<", C="<<hChg
03168 <<", B="<<hBaN<<G4endl;
03169 }
03170 G4cout<<"-EMCLS-G4QInel::Breed: LS r4M="<<s4M-totLS4M<<",rC="<<rCh<<",rB="<<rBN<<G4endl;
03171 #endif
03172
03173 G4int nRes=theResult->size();
03174 if(nRes>2)
03175 {
03176 G4QHadronVector::iterator ih;
03177 G4QHadronVector::iterator nih;
03178 G4QHadronVector::iterator mih;
03179 G4QHadronVector::iterator lst=theResult->end();
03180 lst--;
03181 G4double minMesEn=DBL_MAX;
03182 G4double minBarEn=DBL_MAX;
03183 G4bool nfound=false;
03184 G4bool mfound=false;
03185 for(ih = theResult->begin(); ih < theResult->end(); ++ih) if(ih != lst)
03186 {
03187 G4LorentzVector h4M=(*ih)->Get4Momentum();
03188 G4int hPDG=(*ih)->GetPDGCode();
03189 #ifdef debug
03190 G4cout<<"%->G4QIonIonCollision::Breeder: TRY hPDG="<<hPDG<<", h4M="<<h4M<<G4endl;
03191 #endif
03192 if(hPDG>1111 && hPDG<3333)
03193 {
03194 G4double bE=h4M.e()-(*ih)->GetMass();
03195 if(bE < minBarEn)
03196 {
03197 minBarEn=bE;
03198 nih=ih;
03199 nfound=true;
03200 }
03201 }
03202 else if(hPDG>-1111)
03203 {
03204 G4double mE=h4M.e();
03205 if(mE < minMesEn)
03206 {
03207 minMesEn=mE;
03208 mih=ih;
03209 mfound=true;
03210 }
03211 }
03212 }
03213 if(nfound && mfound && minMesEn+minBarEn < maxEn)
03214 {
03215 G4QHadron* resNuc = (*theResult)[nRes-1];
03216 theResult->pop_back();
03217 G4LorentzVector q4M=(*nih)->Get4Momentum()+(*mih)->Get4Momentum();
03218 G4QContent qQC=(*nih)->GetQC()+(*mih)->GetQC();
03219 maxEn -= minBarEn+minMesEn;
03220 #ifdef debug
03221 G4cout<<"%->G4QIonIonCollision::Breeder:Exclude,4M="<<(*nih)->Get4Momentum()
03222 <<", & 4m="<<(*mih)->Get4Momentum()<<G4endl;
03223 #endif
03224 delete (*nih);
03225 delete (*mih);
03226 if(nih > mih)
03227 {
03228 theResult->erase(nih);
03229 theResult->erase(mih);
03230 }
03231 else
03232 {
03233 theResult->erase(mih);
03234 theResult->erase(nih);
03235 }
03236 #ifdef debug
03237 G4cout<<"%->G4QI::Breed: BeforeLOOP, dE="<<maxEn<<", nR="<<theResult->size()<<G4endl;
03238 #endif
03239 if(maxEn > 0.)
03240 {
03241 for(ih = theResult->begin(); ih < theResult->end(); ih++)
03242 {
03243 G4LorentzVector h4M=(*ih)->Get4Momentum();
03244 G4int hPDG=(*ih)->GetPDGCode();
03245 G4double dE=0.;
03246 if (hPDG> 1111 && hPDG<3333) dE=h4M.e()-(*ih)->GetMass();
03247 else if(hPDG>-1111) dE=h4M.e();
03248 else continue;
03249 if(dE < maxEn)
03250 {
03251 maxEn-=dE;
03252 q4M+=h4M;
03253 qQC+=(*ih)->GetQC();
03254 #ifdef debug
03255 G4cout<<"%->G4QIonIonCollision::Breed:Exclude,4M="<<h4M<<",dE="<<maxEn<<G4endl;
03256 #endif
03257 delete (*ih);
03258 theResult->erase(ih);
03259 --ih;
03260 }
03261 }
03262 }
03263 G4Quasmon* softQuasmon = new G4Quasmon(qQC, q4M);
03264 #ifdef debug
03265 G4cout<<"%->G4QIonIonCollision::Breed:QuasmonIsFilled,4M="<<q4M<<",QC="<<qQC<<G4endl;
03266 #endif
03267 theQuasmons.push_back(softQuasmon);
03268 theResult->push_back(resNuc);
03269 }
03270 #ifdef edebug
03271 G4LorentzVector f4M(0.,0.,0.,0.);
03272 G4int fCh=totChg;
03273 G4int fBN=totBaN;
03274 G4int nHd=theResult->size();
03275 G4int nQm=theQuasmons.size();
03276 G4cout<<"-EMCLS-G4QIonIonCollision::Breeder:#ofHadr="<<nHd<<", #OfQ="<<nQm
03277 <<",rPA="<<rt4M.m()<<"="<<G4QNucleus(rtPDG).GetGSMass()
03278 <<",rTA="<<rt4M.m()<<"="<<G4QNucleus(rtPDG).GetGSMass()<<G4endl;
03279 for(G4int i=0; i<nHd; i++)
03280 {
03281 G4LorentzVector hI4M=(*theResult)[i]->Get4Momentum();
03282 f4M+=hI4M;
03283 G4int hChg=(*theResult)[i]->GetCharge();
03284 fCh-=hChg;
03285 G4int hBaN=(*theResult)[i]->GetBaryonNumber();
03286 fBN-=hBaN;
03287 G4cout<<"-EMCLS-G4QIonIonCollision::Breeder: Hadron#"<<i<<", 4M="<<hI4M<<", PDG="
03288 <<(*theResult)[i]->GetPDGCode()<<", C="<<hChg<<", B="<<hBaN<<G4endl;
03289 }
03290 for(G4int i=0; i<nQm; i++)
03291 {
03292 G4LorentzVector hI4M=theQuasmons[i]->Get4Momentum();
03293 f4M+=hI4M;
03294 G4int hChg=theQuasmons[i]->GetCharge();
03295 fCh-=hChg;
03296 G4int hBaN=theQuasmons[i]->GetBaryonNumber();
03297 fBN-=hBaN;
03298 G4cout<<"-EMCLS-G4QIonIonCollision::Breeder: Quasmon#"<<i<<", 4M="<<hI4M<<", C="
03299 <<hChg<<", B="<<hBaN<<G4endl;
03300 }
03301 G4cout<<"-EMCLS-G4QInel::Breed:, r4M="<<f4M-totLS4M<<", rC="<<fCh<<",rB="<<fBN<<G4endl;
03302 #endif
03303 }
03304 return;
03305 }
03306
03307
03308 G4bool G4QIonIonCollision::ExciteDiffParticipants(G4QHadron* projectile,
03309 G4QHadron* target) const
03310 {
03311 G4LorentzVector Pprojectile=projectile->Get4Momentum();
03312 G4double Mprojectile=projectile->GetMass();
03313 G4double Mprojectile2=Mprojectile*Mprojectile;
03314 G4LorentzVector Ptarget=target->Get4Momentum();
03315 G4double Mtarget=target->GetMass();
03316 G4double Mtarget2=Mtarget*Mtarget;
03317 #ifdef debug
03318 G4cout<<"G4QInel::ExciteDiffPartici: Ep="<<Pprojectile.e()<<", Et="<<Ptarget.e()<<G4endl;
03319 #endif
03320
03321 G4LorentzVector Psum=Pprojectile+Ptarget;
03322 G4LorentzRotation toCms(-Psum.boostVector());
03323 G4LorentzVector Ptmp=toCms*Pprojectile;
03324 if(Ptmp.pz()<=0.)
03325 {
03326 #ifdef debug
03327 G4cout<<"G4QIonIonCollision::ExciteDiffParticipants:*1* abort Collision!! *1*"<<G4endl;
03328 #endif
03329 return false;
03330 }
03331 toCms.rotateZ(-Ptmp.phi());
03332 toCms.rotateY(-Ptmp.theta());
03333 #ifdef debug
03334 G4cout<<"G4QIonIonCollision::ExciteDiffParticipantts:BeforBoost Pproj="<<Pprojectile
03335 <<",Ptarg="<<Ptarget<<G4endl;
03336 #endif
03337 G4LorentzRotation toLab(toCms.inverse());
03338 Pprojectile.transform(toCms);
03339 Ptarget.transform(toCms);
03340 #ifdef debug
03341 G4cout<< "G4QInelast::ExciteDiffParticipantts: AfterBoost Pproj="<<Pprojectile<<",Ptarg="
03342 <<Ptarget<<", cms4M="<<Pprojectile+Ptarget<<G4endl;
03343 G4cout<<"G4QIonIonCollis::ExciteDiffParticipants:ProjX+="<<Pprojectile.plus()<<",ProjX-="
03344 <<Pprojectile.minus()<<", tX+="<< Ptarget.plus()<<",tX-="<<Ptarget.minus()<<G4endl;
03345 #endif
03346 G4LorentzVector Qmomentum(0.,0.,0.,0.);
03347 G4int whilecount=0;
03348 #ifdef debug
03349 G4cout<<"G4QIonIonCollision::ExciteDiffParticipants: Before DO"<<G4endl;
03350 #endif
03351 do
03352 {
03353
03354 G4double maxPtSquare=sqr(Ptarget.pz());
03355 #ifdef debug
03356 G4cout<<"G4QIonIonCollision::ExciteDiffParticipants: maxPtSq="<<maxPtSquare<<G4endl;
03357 if(whilecount++>=500 && whilecount%100==0)
03358 G4cout<<"G4QIonIonCollision::ExciteDiffParticipantts: can loop, loopCount="<<whilecount
03359 <<", maxPtSquare="<<maxPtSquare<<G4endl;
03360 #endif
03361 if(whilecount>1000)
03362 {
03363 #ifdef debug
03364 G4cout<<"G4QIonIonCollision::ExciteDiffParticipants: *2* abort Loop!! *2*"<<G4endl;
03365 #endif
03366 return false;
03367 }
03368 Qmomentum=G4LorentzVector(GaussianPt(widthOfPtSquare,maxPtSquare),0);
03369 #ifdef debug
03370 G4cout<<"G4QIonIonCollision::ExciteDiffParticipants: generated Pt="<<Qmomentum
03371 <<", ProjPt="<<Pprojectile+Qmomentum<<", TargPt="<<Ptarget-Qmomentum<<G4endl;
03372 #endif
03373
03374 G4double Xmin=0.;
03375 G4double Xmax=1.;
03376 G4double Xplus =ChooseX(Xmin,Xmax);
03377 G4double Xminus=ChooseX(Xmin,Xmax);
03378 #ifdef debug
03379 G4cout<<"G4QIonIonCollision::ExciteDiffParticipant:X+="<<Xplus<<",X-="<<Xminus<<G4endl;
03380 #endif
03381 G4double pt2=Qmomentum.vect().mag2();
03382 G4double Qplus =-pt2/Xminus/Ptarget.minus();
03383 G4double Qminus= pt2/Xplus /Pprojectile.plus();
03384 Qmomentum.setPz((Qplus-Qminus)/2);
03385 Qmomentum.setE( (Qplus+Qminus)/2);
03386 #ifdef debug
03387 G4cout<<"G4QInelast::ExciteDiffParticip: Qplus="<<Qplus<<", Qminus="<<Qminus<<", pt2="
03388 <<pt2<<", Qmomentum="<<Qmomentum<<", ProjM="<<(Pprojectile+Qmomentum).mag()
03389 <<", TargM="<<(Ptarget-Qmomentum).mag()<<G4endl;
03390 #endif
03391 } while((Pprojectile+Qmomentum).mag2()<=Mprojectile2 ||
03392 (Ptarget-Qmomentum).mag2()<=Mtarget2);
03393 Pprojectile += Qmomentum;
03394 Ptarget -= Qmomentum;
03395 #ifdef debug
03396 G4cout<<"G4QInelast::ExciteDiffParticipan: Proj(Q)="<<Pprojectile<<", Targ(Q)="<<Ptarget
03397 <<", Proj(back)="<<toLab*Pprojectile<<", Targ(bac)="<< toLab*Ptarget << G4endl;
03398 #endif
03399
03400 Pprojectile.transform(toLab);
03401 Ptarget.transform(toLab);
03402 #ifdef debug
03403 G4cout<< "G4QIonIonCollision::ExciteDiffParticipants:TargetMass="<<Ptarget.mag()<<G4endl;
03404 #endif
03405 target->Set4Momentum(Ptarget);
03406 #ifdef debug
03407 G4cout<<"G4QIonIonCollision::ExciteDiffParticipant:ProjMass="<<Pprojectile.mag()<<G4endl;
03408 #endif
03409 projectile->Set4Momentum(Pprojectile);
03410 return true;
03411 }
03412
03413
03414
03415 G4bool G4QIonIonCollision::ExciteSingDiffParticipants(G4QHadron* projectile,
03416 G4QHadron* target) const
03417 {
03418 G4LorentzVector Pprojectile=projectile->Get4Momentum();
03419 G4double Mprojectile=projectile->GetMass();
03420 G4double Mprojectile2=Mprojectile*Mprojectile;
03421 G4LorentzVector Ptarget=target->Get4Momentum();
03422 G4double Mtarget=target->GetMass();
03423 G4double Mtarget2=Mtarget*Mtarget;
03424 #ifdef debug
03425 G4cout<<"G4QInel::ExSingDiffPartici: Ep="<<Pprojectile.e()<<", Et="<<Ptarget.e()<<G4endl;
03426 #endif
03427 G4bool KeepProjectile= G4UniformRand() > 0.5;
03428
03429 if(KeepProjectile )
03430 {
03431 #ifdef debug
03432 G4cout<<"-1/2-G4QIonIonCollision::ExSingDiffParticipants: Projectile is fixed"<<G4endl;
03433 #endif
03434 Mprojectile2 = projectile->GetMass2()*(1.-perCent);
03435 }
03436 else
03437 {
03438 #ifdef debug
03439 G4cout<<"---1/2---G4QIonIonCollision::ExSingDiffParticipants: Target is fixed"<<G4endl;
03440 #endif
03441 Mtarget2 = target->GetMass2()*(1.-perCent);
03442 }
03443
03444
03445 G4LorentzVector Psum=Pprojectile+Ptarget;
03446 G4LorentzRotation toCms(-Psum.boostVector());
03447 G4LorentzVector Ptmp=toCms*Pprojectile;
03448 if(Ptmp.pz()<=0.)
03449 {
03450 #ifdef debug
03451 G4cout<<"G4QIonIonCollision::ExciteSingDiffParticipants: *1*abortCollision*1*"<<G4endl;
03452 #endif
03453 return false;
03454 }
03455 toCms.rotateZ(-Ptmp.phi());
03456 toCms.rotateY(-Ptmp.theta());
03457 #ifdef debug
03458 G4cout<<"G4QInel::ExciteSingDiffParticipantts: Be4Boost Pproj="<<Pprojectile<<", Ptarg="
03459 <<Ptarget<<G4endl;
03460 #endif
03461 G4LorentzRotation toLab(toCms.inverse());
03462 Pprojectile.transform(toCms);
03463 Ptarget.transform(toCms);
03464 #ifdef debug
03465 G4cout<< "G4QInelast::ExciteDiffParticipantts: AfterBoost Pproj="<<Pprojectile<<"Ptarg="
03466 <<Ptarget<<", cms4M="<<Pprojectile+Ptarget<<G4endl;
03467
03468 G4cout<<"G4QInelast::ExciteDiffParticipantts: ProjX+="<<Pprojectile.plus()<<", ProjX-="
03469 <<Pprojectile.minus()<<", TargX+="<< Ptarget.plus()<<", TargX-="<<Ptarget.minus()
03470 <<G4endl;
03471 #endif
03472 G4LorentzVector Qmomentum(0.,0.,0.,0.);
03473 G4int whilecount=0;
03474 do
03475 {
03476
03477 G4double maxPtSquare=sqr(Ptarget.pz());
03478 if(whilecount++>=500 && whilecount%100==0)
03479 #ifdef debug
03480 G4cout<<"G4QIonIonCollision::ExciteSingDiffParticipantts: can loop, loopCount="
03481 <<whilecount<<", maxPtSquare="<<maxPtSquare<<G4endl;
03482 #endif
03483 if(whilecount>1000)
03484 {
03485 #ifdef debug
03486 G4cout<<"G4QIonIonCollision::ExciteSingDiffParticipants:*2* abortLoop!! *2*"<<G4endl;
03487 #endif
03488 return false;
03489 }
03490 Qmomentum=G4LorentzVector(GaussianPt(widthOfPtSquare,maxPtSquare),0);
03491 #ifdef debug
03492 G4cout<<"G4QInelast::ExciteSingDiffParticipants: generated Pt="<<Qmomentum
03493 <<", ProjPt="<<Pprojectile+Qmomentum<<", TargPt="<<Ptarget-Qmomentum<<G4endl;
03494 #endif
03495
03496 G4double Xmin=0.;
03497 G4double Xmax=1.;
03498 G4double Xplus =ChooseX(Xmin,Xmax);
03499 G4double Xminus=ChooseX(Xmin,Xmax);
03500 #ifdef debug
03501 G4cout<<"G4QInel::ExciteSingDiffPartici: X-plus="<<Xplus<<", X-minus="<<Xminus<<G4endl;
03502 #endif
03503 G4double pt2=G4ThreeVector(Qmomentum.vect()).mag2();
03504 G4double Qplus =-pt2/Xminus/Ptarget.minus();
03505 G4double Qminus= pt2/Xplus /Pprojectile.plus();
03506 if (KeepProjectile)
03507 Qminus=(projectile->GetMass2()+pt2)/(Pprojectile.plus()+Qplus) - Pprojectile.minus();
03508 else Qplus=Ptarget.plus() - (target->GetMass2()+pt2)/(Ptarget.minus()-Qminus);
03509 Qmomentum.setPz((Qplus-Qminus)/2);
03510 Qmomentum.setE( (Qplus+Qminus)/2);
03511 #ifdef debug
03512 G4cout<<"G4QInel::ExciteDiffParticip: Qplus="<<Qplus<<", Qminus="<<Qminus<<", pt2="
03513 <<pt2<<", Qmomentum="<<Qmomentum<<", ProjM="<<(Pprojectile+Qmomentum).mag()
03514 <<", TargM="<<(Ptarget-Qmomentum).mag()<<G4endl;
03515 #endif
03516
03517
03518
03519 } while((Ptarget-Qmomentum).mag2()<=Mtarget2 ||
03520 (Pprojectile+Qmomentum).mag2()<=Mprojectile2 ||
03521 (Ptarget-Qmomentum).e() < 0. || (Pprojectile+Qmomentum).e() < 0.);
03522 Pprojectile += Qmomentum;
03523 Ptarget -= Qmomentum;
03524 #ifdef debug
03525 G4cout<<"G4QIonIonCollision::ExciteSingDiffParticipan: Proj(Q)="<<Pprojectile<<"(E="
03526 <<Pprojectile.e()<<"), Targ(Q)="<<Ptarget<<"(E="<<Ptarget.e()
03527 <<"), Proj(back)="<<toLab*Pprojectile<<", Targ(bac)="<< toLab*Ptarget << G4endl;
03528 #endif
03529
03530 Pprojectile.transform(toLab);
03531 Ptarget.transform(toLab);
03532 #ifdef debug
03533 G4cout<< "G4QIonIonCollision::ExciteSingleDiffParticipants: TgM="<<Ptarget.mag()<<G4endl;
03534 #endif
03535 target->Set4Momentum(Ptarget);
03536 #ifdef debug
03537 G4cout<<"G4QIonIonCollision::ExciteSingleParticipants:ProjM="<<Pprojectile.mag()<<G4endl;
03538 #endif
03539 projectile->Set4Momentum(Pprojectile);
03540 return true;
03541 }
03542
03543 void G4QIonIonCollision::SetParameters(G4int nC,G4double stT, G4double tbD, G4double SigPt)
03544 {
03545 nCutMax = nC;
03546 stringTension = stT;
03547 tubeDensity = tbD;
03548 widthOfPtSquare = -2*SigPt*SigPt;
03549 }
03550
03551 G4double G4QIonIonCollision::ChooseX(G4double Xmin, G4double Xmax) const
03552 {
03553
03554
03555 if(Xmax == Xmin) return Xmin;
03556 if( Xmin < 0. || Xmax < Xmin)
03557 {
03558 G4cerr<<"***G4QIonIonCollision::ChooseX: Xmin="<<Xmin<<", Xmax="<<Xmax<< G4endl;
03559 G4Exception("G4QIonIonCollision::ChooseX:","72",FatalException,"Bad X or X-Range");
03560 }
03561
03562
03563 G4double sxi=std::sqrt(Xmin);
03564 G4double x=sqr(sxi+G4UniformRand()*(std::sqrt(Xmax)-sxi));
03565 #ifdef debug
03566 G4cout<<"G4QIonIonCollision::ChooseX: DiffractiveX="<<x<<G4endl;
03567 #endif
03568 return x;
03569 }
03570
03571
03572 G4ThreeVector G4QIonIonCollision::GaussianPt(G4double widthSq, G4double maxPtSquare) const
03573 {
03574 #ifdef debug
03575 G4cout<<"G4QIonIonCollision::GaussianPt: w^2="<<widthSq<<",maxPt2="<<maxPtSquare<<G4endl;
03576 #endif
03577 G4double pt2=0.;
03578 G4double rm=maxPtSquare/widthSq;
03579 if(rm>-.01) pt2=widthSq*(std::sqrt(1.-G4UniformRand()*(1.-sqr(1.+rm)))-1.);
03580 else pt2=widthSq*std::log(1.-G4UniformRand()*(1.-std::exp(rm)));
03581 pt2=std::sqrt(pt2);
03582 G4double phi=G4UniformRand()*twopi;
03583 return G4ThreeVector(pt2*std::cos(phi),pt2*std::sin(phi),0.);
03584 }
03585
03586 G4int G4QIonIonCollision::SumPartonPDG(G4int PDG1, G4int PDG2) const
03587 {
03588 if (PDG1 < 7 && PDG1 > 0 && PDG2 < 7 && PDG2 > 0)
03589 {
03590 if(PDG1 > PDG2) return PDG1*1000+PDG2*100+1;
03591 else return PDG2*1000+PDG1*100+1;
03592 }
03593 else if (PDG1 >-7 && PDG1 < 0 && PDG2 >-7 && PDG2 < 0)
03594 {
03595 if(-PDG1 > -PDG2) return PDG1*1000+PDG2*100-1;
03596 else return PDG2*1000+PDG1*100-1;
03597 }
03598 else if (PDG1 <-99 && PDG2 < 7 && PDG2 > 0)
03599 {
03600 G4int PDG=-PDG1/100;
03601 if(PDG2==PDG/10) return -PDG%10;
03602 if(PDG2==PDG%10) return -PDG/10;
03603 else
03604 {
03605 G4cerr<<"***4QFragmentation::SumPartonPDG: PDG1="<<PDG1<<", PDG2="<<PDG2<<G4endl;
03606 G4Exception("G4QIonIonCollision::SumPartonPDG:","72",FatalException,"Q&ADiQnoMatch");
03607 }
03608 }
03609 else if (PDG2 <-99 && PDG1 < 7 && PDG1 > 0)
03610 {
03611 G4int PDG=-PDG2/100;
03612 if(PDG1==PDG/10) return -PDG%10;
03613 if(PDG1==PDG%10) return -PDG/10;
03614 else
03615 {
03616 G4cerr<<"***4QFragmentation::SumPartonPDG: PDG1="<<PDG1<<", PDG2="<<PDG2<<G4endl;
03617 G4Exception("G4QIonIonCollision::SumPartonPDG:","72",FatalException,"ADiQ&QnoMatch");
03618 }
03619 }
03620 else if (PDG1 > 99 && PDG2 >-7 && PDG2 < 0)
03621 {
03622 G4int PDG=PDG1/100;
03623 if(PDG2==-PDG/10) return PDG%10;
03624 if(PDG2==-PDG%10) return PDG/10;
03625 else
03626 {
03627 G4cerr<<"***4QFragmentation::SumPartonPDG: PDG1="<<PDG1<<", PDG2="<<PDG2<<G4endl;
03628 G4Exception("G4QIonIonCollision::SumPartonPDG:","72",FatalException,"DiQ&AQnoMatch");
03629 }
03630 }
03631 else if (PDG2 > 99 && PDG1 >-7 && PDG1 < 0)
03632 {
03633 G4int PDG=PDG2/100;
03634 if(PDG1==-PDG/10) return PDG%10;
03635 if(PDG1==-PDG%10) return PDG/10;
03636 else
03637 {
03638 G4cerr<<"***G4QIonIonCollision::SumPartonPDG: PDG1="<<PDG1<<", PDG2="<<PDG2<<G4endl;
03639 G4Exception("G4QIonIonCollision::SumPartonPDG:","72",FatalException,"AQ&DiQnoMatch");
03640 }
03641 }
03642 else
03643 {
03644 G4cerr<<"***G4QIonIonCollision::SumPartonPDG: PDG1="<<PDG1<<", PDG2="<<PDG2<<G4endl;
03645 G4Exception("G4QIonIonCollision::SumPartonPDG:","72",FatalException,"Can'tSumUpParts");
03646 }
03647 return 0;
03648 }
03649
03650
03651 std::pair<G4int,G4int> G4QIonIonCollision::ReducePair(G4int P1, G4int P2) const
03652 {
03653 #ifdef debug
03654 G4cout<<"G4QIonIonCollision::ReducePair: **Called** P1="<<P1<<", P2="<<P2<<G4endl;
03655 #endif
03656 G4int P11=P1/10;
03657 G4int P12=P1%10;
03658 G4int P21=P2/10;
03659 G4int P22=P2%10;
03660 if (P11==P21) return std::make_pair(P12,P22);
03661 else if(P11==P22) return std::make_pair(P12,P21);
03662 else if(P12==P21) return std::make_pair(P11,P22);
03663 else if(P12==P22) return std::make_pair(P11,P21);
03664
03665 G4cout<<"-Warning-G4QIonIonCollision::ReducePair:**Failed**,P1="<<P1<<",P2="<<P2<<G4endl;
03666
03667 return std::make_pair(0,0);
03668 }
03669
03670
03671 G4int G4QIonIonCollision::AnnihilationOrder(G4int LS, G4int MPS, G4int uPDG, G4int mPDG,
03672 G4int sPDG, G4int nPDG)
03673 {
03674 G4int Ord=0;
03675
03676 if (LS==2 && MPS==2 )
03677 {
03678 #ifdef debug
03679 G4cout<<"G4QIonIonCollision::AnnihilationOrder:QaQ/QaQ->DiQ-aDiQ, uPDG="<<uPDG<<G4endl;
03680 #endif
03681 if ( (uPDG>0 && sPDG>0 && mPDG<0 && nPDG<0) ||
03682 (uPDG<0 && sPDG<0 && mPDG>0 && nPDG>0) ) Ord= 1;
03683 else if( (uPDG>0 && nPDG>0 && mPDG<0 && sPDG<0) ||
03684 (uPDG<0 && nPDG<0 && mPDG>0 && sPDG>0) ) Ord=-1;
03685 else G4cerr<<"-Warning-G4QIonIonCollision::AnnihilationOrder: Wrong 22 fusion, L="
03686 <<uPDG<<",R="<<mPDG<<",cL="<<sPDG<<",cR="<<nPDG<<G4endl;
03687 }
03688 else if ( LS==2 && MPS==3 )
03689 {
03690 if (uPDG > 7)
03691 {
03692 #ifdef debug
03693 G4cout<<"G4QIonIonCollision::AnnihOrder:pLDiQ, sPDG="<<sPDG<<", nPDG="<<nPDG<<G4endl;
03694 #endif
03695 if ( sPDG<0 && (-sPDG==uPDG/1000 || -sPDG==(uPDG/100)%10) ) Ord= 1;
03696 else if( nPDG<0 && (-nPDG==uPDG/1000 || -nPDG==(uPDG/100)%10) ) Ord=-1;
03697 else G4cerr<<"-Warning-G4QIonIonCollision::AnnihilationOrder: Wrong pLDiQ, L="<<uPDG
03698 <<",R="<<mPDG<<",cL="<<sPDG<<",cR="<<nPDG<<G4endl;
03699 }
03700 else if (mPDG > 7)
03701 {
03702 #ifdef debug
03703 G4cout<<"G4QIonIonCollision::AnnihOrder:pRDiQ, sPDG="<<sPDG<<", nPDG="<<nPDG<<G4endl;
03704 #endif
03705 if ( sPDG<0 && (-sPDG==mPDG/1000 || -sPDG==(mPDG/100)%10) ) Ord=-1;
03706 else if( nPDG<0 && (-nPDG==mPDG/1000 || -nPDG==(mPDG/100)%10) ) Ord= 1;
03707 else G4cerr<<"-Warning-G4QIonIonCollision::AnnihilationOrder: Wrong pRDiQ, L="<<uPDG
03708 <<",R="<<mPDG<<",cL="<<sPDG<<",cR="<<nPDG<<G4endl;
03709 }
03710 else if (uPDG <-7)
03711 {
03712 #ifdef debug
03713 G4cout<<"G4QIonIonCollision::AnnihOrder:pLaDiQ, sPDG="<<sPDG<<",nPDG="<<nPDG<<G4endl;
03714 #endif
03715 if ( sPDG>0 && (sPDG==(-uPDG)/1000 || sPDG==((-uPDG)/100)%10) ) Ord= 1;
03716 else if( nPDG>0 && (nPDG==(-uPDG)/1000 || nPDG==((-uPDG)/100)%10) ) Ord=-1;
03717 else G4cerr<<"-Warning-G4QIonIonCollision::AnnihilationOrder: Wrong pLaDiQ, L="<<uPDG
03718 <<", R="<<mPDG<<", cL="<<sPDG<<", cR="<<nPDG<<G4endl;
03719 }
03720 else if (mPDG <-7)
03721 {
03722 #ifdef debug
03723 G4cout<<"G4QIonIonCollision::AnnihOrder:pRaDiQ, sPDG="<<sPDG<<",nPDG="<<nPDG<<G4endl;
03724 #endif
03725 if ( sPDG>0 && (sPDG==(-mPDG)/1000 || sPDG==((-mPDG)/100)%10) ) Ord=-1;
03726 else if( nPDG>0 && (nPDG==(-mPDG)/1000 || nPDG==((-mPDG)/100)%10) ) Ord= 1;
03727 else G4cerr<<"-Warning-G4QIonIonCollision::AnnihilationOrder: Wrong pRaDiQ, L="<<uPDG
03728 <<", R="<<mPDG<<", cL="<<sPDG<<", cR="<<nPDG<<G4endl;
03729 }
03730 else if( (sPDG<0 && (-sPDG==mPDG || -sPDG==uPDG) ) ||
03731 (nPDG<0 && (-nPDG==mPDG || -nPDG==uPDG) ) ) Ord= 2;
03732 #ifdef debug
03733 else G4cout<<"-Warning-G4QIonIonCollision::AnnihilatOrder:Wrong23StringFusion"<<G4endl;
03734 G4cout<<"G4QIonIonCollision::AnnihilationOrder: Ord="<<Ord<<",sPDG="<<sPDG<<",nPDG="
03735 <<nPDG<<", uPDG="<<uPDG<<",mPDG="<<mPDG<<G4endl;
03736 #endif
03737 }
03738 else if ( LS==3 && MPS==2 )
03739 {
03740 if (sPDG > 7)
03741 {
03742 #ifdef debug
03743 G4cout<<"G4QIonIonCollision::AnnihOrder:cLDiQ, uPDG="<<uPDG<<", mPDG="<<mPDG<<G4endl;
03744 #endif
03745 if ( uPDG<0 && (-uPDG==sPDG/1000 || -uPDG==(sPDG/100)%10) ) Ord= 1;
03746 else if( mPDG<0 && (-mPDG==sPDG/1000 || -mPDG==(sPDG/100)%10) ) Ord=-1;
03747 else G4cerr<<"-Warning-G4QIonIonCollision::AnnihilationOrder: Wrong cLDiQ, L="<<uPDG
03748 <<",R="<<mPDG<<",cL="<<sPDG<<",cR="<<nPDG<<G4endl;
03749 }
03750 else if (nPDG > 7)
03751 {
03752 #ifdef debug
03753 G4cout<<"G4QIonIonCollision::AnnihOrder:cRDiQ, uPDG="<<uPDG<<", mPDG="<<mPDG<<G4endl;
03754 #endif
03755 if ( uPDG<0 && (-uPDG==nPDG/1000 || -uPDG==(nPDG/100)%10) ) Ord=-1;
03756 else if( mPDG<0 && (-mPDG==nPDG/1000 || -mPDG==(nPDG/100)%10) ) Ord= 1;
03757 else G4cerr<<"-Warning-G4QIonIonCollision::AnnihilationOrder: Wrong cRDiQ, L="<<uPDG
03758 <<",R="<<mPDG<<",cL="<<sPDG<<",cR="<<nPDG<<G4endl;
03759 }
03760 else if (sPDG <-7)
03761 {
03762 #ifdef debug
03763 G4cout<<"G4QIonIonCollision::AnnihOrder:cLaDiQ, uPDG="<<uPDG<<",mPDG="<<mPDG<<G4endl;
03764 #endif
03765 if ( uPDG>0 && (uPDG==(-sPDG)/1000 || uPDG==((-sPDG)/100)%10) ) Ord= 1;
03766 else if( mPDG>0 && (mPDG==(-sPDG)/1000 || mPDG==((-sPDG)/100)%10) ) Ord=-1;
03767 else G4cerr<<"-Warning-G4QIonIonCollision::AnnihilationOrder: Wrong cLaDiQ, L="<<uPDG
03768 <<", R="<<mPDG<<", cL="<<sPDG<<", cR="<<nPDG<<G4endl;
03769 }
03770 else if (nPDG <-7)
03771 {
03772 #ifdef debug
03773 G4cout<<"G4QIonIonCollision::AnnihOrder:cRaDiQ, uPDG="<<uPDG<<",mPDG="<<mPDG<<G4endl;
03774 #endif
03775 if ( uPDG>0 && (uPDG==(-nPDG)/1000 || uPDG==((-nPDG)/100)%10) ) Ord=-1;
03776 else if( mPDG>0 && (mPDG==(-nPDG)/1000 || mPDG==((-nPDG)/100)%10) ) Ord= 1;
03777 else G4cerr<<"-Warning-G4QIonIonCollision::AnnihilationOrder: Wrong cRaDiQ, L="<<uPDG
03778 <<", R="<<mPDG<<", cL="<<sPDG<<", cR="<<nPDG<<G4endl;
03779 }
03780 else if( (uPDG<0 && (-uPDG==sPDG || -uPDG==nPDG) ) ||
03781 (mPDG<0 && (-mPDG==sPDG || -mPDG==nPDG) ) ) Ord=2;
03782 #ifdef debug
03783 else G4cout<<"-Warning-G4QIonIonCollision::AnnihilatOrder:Wrong32StringFusion"<<G4endl;
03784 G4cout<<"G4QIonIonCollision::AnnihilationOrder: Ord="<<Ord<<",sPDG="<<sPDG<<",nPDG="
03785 <<nPDG<<", uPDG="<<uPDG<<",mPDG="<<mPDG<<G4endl;
03786 #endif
03787 }
03788 else if ( (LS==2 && MPS==4) || (LS==4 && MPS==2) )
03789 {
03790 if (uPDG > 7)
03791 {
03792 #ifdef debug
03793 G4cout<<"G4QIonIonCollision::AnnihilOrder:pLDiQ,sPDG="<<sPDG<<",nPDG="<<nPDG<<G4endl;
03794 #endif
03795 if ( sPDG<0 && (-sPDG==uPDG/1000 || -sPDG==(uPDG/100)%10) &&
03796 (nPDG==(-mPDG)/1000 || nPDG==((-mPDG)/100)%10) ) Ord= 1;
03797 else if( nPDG<0 && (-nPDG==uPDG/1000 || -nPDG==(uPDG/100)%10) &&
03798 (sPDG==(-mPDG)/1000 || sPDG==((-mPDG)/100)%10) ) Ord=-1;
03799 else G4cerr<<"-Warning-G4QIonIonCollision::AnnihilationOrder: Wrong 24 pLDiQ, L="
03800 <<uPDG<<",R="<<mPDG<<",cL="<<sPDG<<",cR="<<nPDG<<G4endl;
03801 }
03802 else if (mPDG >7)
03803 {
03804 #ifdef debug
03805 G4cout<<"G4QIonIonCollision::AnnihilOrder:PRDiQ,sPDG="<<sPDG<<",nPDG="<<nPDG<<G4endl;
03806 #endif
03807 if ( sPDG<0 && (-sPDG==mPDG/1000 || -sPDG==(mPDG/100)%10) &&
03808 (nPDG==(-uPDG)/1000 || nPDG==((-uPDG)/100)%10) ) Ord=-1;
03809 else if( nPDG<0 && (-nPDG==mPDG/1000 || -nPDG==(mPDG/100)%10) &&
03810 (sPDG==(-uPDG)/1000 || sPDG==((-uPDG)/100)%10) ) Ord= 1;
03811 else G4cerr<<"-Warning-G4QIonIonCollision::AnnihilationOrder: Wrong 24 pLDiQ, L="
03812 <<uPDG<<",R="<<mPDG<<",cL="<<sPDG<<",cR="<<nPDG<<G4endl;
03813 }
03814 else if (sPDG > 7)
03815 {
03816 #ifdef debug
03817 G4cout<<"G4QIonIonCollision::AnnihilOrder:cLDiQ,uPDG="<<uPDG<<",mPDG="<<mPDG<<G4endl;
03818 #endif
03819 if ( uPDG<0 && (-uPDG==sPDG/1000 || -uPDG==(sPDG/100)%10) &&
03820 (mPDG==(-nPDG)/1000 || mPDG==((-nPDG)/100)%10) ) Ord= 1;
03821 else if( mPDG<0 && (-mPDG==sPDG/1000 || -mPDG==(sPDG/100)%10) &&
03822 (uPDG==(-nPDG)/1000 || uPDG==((-nPDG)/100)%10) ) Ord=-1;
03823 else G4cerr<<"-Warning-G4QIonIonCollision::AnnihilationOrder: Wrong 24 cLDiQ, L="
03824 <<uPDG<<",R="<<mPDG<<",cL="<<sPDG<<",cR="<<nPDG<<G4endl;
03825 }
03826 else if (nPDG > 7)
03827 {
03828 #ifdef debug
03829 G4cout<<"G4QIonIonCollision::AnnihilOrder:cRDiQ,uPDG="<<uPDG<<",mPDG="<<mPDG<<G4endl;
03830 #endif
03831 if ( uPDG<0 && (-uPDG==nPDG/1000 || -uPDG==(nPDG/100)%10) &&
03832 (mPDG==(-sPDG)/1000 || mPDG==((-sPDG)/100)%10) ) Ord=-1;
03833 else if( mPDG<0 && (-mPDG==nPDG/1000 || -mPDG==(nPDG/100)%10) &&
03834 (uPDG==(-sPDG)/1000 || uPDG==((-sPDG)/100)%10) ) Ord= 1;
03835 else G4cerr<<"-Warning-G4QIonIonCollision::AnnihilationOrder: Wrong 24 cRDiQ, L="
03836 <<uPDG<<",R="<<mPDG<<",cL="<<sPDG<<",cR="<<nPDG<<G4endl;
03837 }
03838 #ifdef debug
03839 else G4cout<<"-Warning-G4QIonIonCollision::AnnihilOrder: Wrong 24StringFusion"<<G4endl;
03840 G4cout<<"G4QIonIonCollision::AnnihilationOrder: Ord="<<Ord<<",sPDG="<<sPDG<<",nPDG="
03841 <<nPDG<<", uPDG="<<uPDG<<",mPDG="<<mPDG<<G4endl;
03842 #endif
03843 }
03844 else if ( LS==3 && MPS==3 )
03845 {
03846 if (uPDG > 7)
03847 {
03848 #ifdef debug
03849 G4cout<<"G4QIonIonCollision::AnnihilOrder:pLDiQ,sPDG="<<sPDG<<",nPDG="<<nPDG<<G4endl;
03850 #endif
03851 if ( sPDG<-7 && (-nPDG==uPDG/1000 || -nPDG==(uPDG/100)%10) &&
03852 (mPDG==(-sPDG)/1000 || mPDG==((-sPDG)/100)%10) ) Ord=-1;
03853 else if( nPDG<-7 && (-sPDG==uPDG/1000 || -sPDG==(uPDG/100)%10) &&
03854 (mPDG==(-nPDG)/1000 || mPDG==((-nPDG)/100)%10) ) Ord= 1;
03855 else G4cerr<<"-Warning-G4QIonIonCollision::AnnihilationOrder: Wrong 33 pLDiQ, L="
03856 <<uPDG<<",R="<<mPDG<<",cL="<<sPDG<<",cR="<<nPDG<<G4endl;
03857 }
03858 else if(mPDG > 7)
03859 {
03860 #ifdef debug
03861 G4cout<<"G4QIonIonCollision::AnnihilOrder:pRDiQ,sPDG="<<sPDG<<",nPDG="<<nPDG<<G4endl;
03862 #endif
03863 if ( sPDG<-7 && (-nPDG==mPDG/1000 || -nPDG==(mPDG/100)%10) &&
03864 (uPDG==(-sPDG)/1000 || uPDG==((-sPDG)/100)%10) ) Ord= 1;
03865 else if( nPDG<-7 && (-sPDG==mPDG/1000 || -sPDG==(mPDG/100)%10) &&
03866 (uPDG==(-nPDG)/1000 || uPDG==((-nPDG)/100)%10) ) Ord=-1;
03867 else G4cerr<<"-Warning-G4QIonIonCollision::AnnihilationOrder: Wrong33pRDiQ, L="<<uPDG
03868 <<",R="<<mPDG<<",cL="<<sPDG<<",cR="<<nPDG<<G4endl;
03869 }
03870 else if(sPDG > 7)
03871 {
03872 #ifdef debug
03873 G4cout<<"G4QIonIonCollision::AnnihilOrder:cLDiQ,uPDG="<<uPDG<<",mPDG="<<mPDG<<G4endl;
03874 #endif
03875 if ( uPDG<-7 && (-mPDG==sPDG/1000 || -mPDG==(sPDG/100)%10) &&
03876 (nPDG==(-uPDG)/1000 || nPDG==((-uPDG)/100)%10) ) Ord=-1;
03877 else if( mPDG<-7 && (-uPDG==sPDG/1000 || -uPDG==(sPDG/100)%10) &&
03878 (nPDG==(-mPDG)/1000 || nPDG==((-mPDG)/100)%10) ) Ord= 1;
03879 else G4cerr<<"-Warning-G4QIonIonCollision::AnnihilationOrder: Wrong 33 cLDiQ, L="
03880 <<uPDG<<",R="<<mPDG<<",cL="<<sPDG<<",cR="<<nPDG<<G4endl;
03881 }
03882 else if(nPDG > 7)
03883 {
03884 #ifdef debug
03885 G4cout<<"G4QIonIonCollision::AnnihilOrder:cRDiQ,uPDG="<<uPDG<<",mPDG="<<mPDG<<G4endl;
03886 #endif
03887 if ( uPDG<-7 && (-mPDG==nPDG/1000 || -mPDG==(nPDG/100)%10) &&
03888 (nPDG==(-uPDG)/1000 || nPDG==((-uPDG)/100)%10) ) Ord= 1;
03889 else if( mPDG<-7 && (-uPDG==nPDG/1000 || -sPDG==(nPDG/100)%10) &&
03890 (sPDG==(-mPDG)/1000 || sPDG==((-mPDG)/100)%10) ) Ord=-1;
03891 else G4cerr<<"-Warning-G4QIonIonCollision::AnnihilationOrder: Wrong 33 cRDiQ, L="
03892 <<uPDG<<",R="<<mPDG<<",cL="<<sPDG<<",cR="<<nPDG<<G4endl;
03893 }
03894 #ifdef debug
03895 else G4cout<<"-Warning-G4QIonIonCollision::AnnihilOrder: Wrong 33StringFusion"<<G4endl;
03896 G4cout<<"G4QIonIonCollision::AnnihilationOrder: Ord="<<Ord<<",sPDG="<<sPDG<<",nPDG="
03897 <<nPDG<<", uPDG="<<uPDG<<",mPDG="<<mPDG<<G4endl;
03898 #endif
03899 }
03900 return Ord;
03901 }
03902
03903 void G4QIonIonCollision::SwapPartons()
03904 {
03905 static const G4double baryM=800.;
03906 G4QStringVector::iterator ist;
03907 for(ist = strings.begin(); ist < strings.end(); ist++)
03908 {
03909 G4LorentzVector cS4M=(*ist)->Get4Momentum();
03910 G4double cSM2=cS4M.m2();
03911 if(cSM2<0.1)
03912 {
03913 G4QParton* cLeft=(*ist)->GetLeftParton();
03914 G4QParton* cRight=(*ist)->GetRightParton();
03915 G4int cLPDG=cLeft->GetPDGCode();
03916 G4int cRPDG=cRight->GetPDGCode();
03917 G4LorentzVector cL4M=cLeft->Get4Momentum();
03918 G4LorentzVector cR4M=cRight->Get4Momentum();
03919 G4int cLT=cLeft->GetType();
03920 G4int cRT=cRight->GetType();
03921 G4QStringVector::iterator sst;
03922 G4QStringVector::iterator pst;
03923 G4double maxM=-DBL_MAX;
03924 G4int sDir=0;
03925 for(pst = strings.begin(); pst < strings.end(); pst++) if(pst != ist)
03926 {
03927 G4QParton* pLeft=(*pst)->GetLeftParton();
03928 G4QParton* pRight=(*pst)->GetRightParton();
03929 G4int pLPDG=pLeft->GetPDGCode();
03930 G4int pRPDG=pRight->GetPDGCode();
03931 G4LorentzVector pL4M=pLeft->Get4Momentum();
03932 G4LorentzVector pR4M=pRight->Get4Momentum();
03933 G4int pLT=pLeft->GetType();
03934 G4int pRT=pRight->GetType();
03935 G4double LM=0.;
03936 G4double RM=0.;
03937 if(((cLPDG<-7 || (cLPDG>0 && cLPDG< 7) ) && (pLPDG<-7 || (pLPDG>0 && pLPDG< 7) ))||
03938 ((cLPDG> 7 || (cLPDG<0 && cLPDG>-7) ) && (pLPDG> 7 || (pLPDG<0 && cLPDG>-7) )))
03939 {
03940 G4double pLM2=(cL4M+pR4M).m2();
03941 G4double cLM2=(cR4M+pL4M).m2();
03942 if(pLM2>0. && cLM2>0.)
03943 {
03944 G4double pLM=std::sqrt(pLM2);
03945 if(cLT+pRT==3) pLM-=baryM;
03946 G4double cLM=std::sqrt(cLM2);
03947 if(cRT+pLT==3) cLM-=baryM;
03948 LM=std::min(pLM2,cLM2);
03949 }
03950 }
03951 if(((cRPDG<-7 || (cRPDG>0 && cRPDG< 7) ) && (pRPDG<-7 || (pRPDG>0 && pRPDG< 7) ))||
03952 ((cRPDG> 7 || (cRPDG<0 && cRPDG>-7) ) && (pRPDG> 7 || (pRPDG<0 && cRPDG>-7) )) )
03953 {
03954 G4double pRM2=(cR4M+pL4M).m2();
03955 G4double cRM2=(cL4M+pR4M).m2();
03956 if(pRM2>0. && cRM2>0.)
03957 {
03958 G4double pRM=std::sqrt(pRM2);
03959 if(cRT+pLT==3) pRM-=baryM;
03960 G4double cRM=std::sqrt(cRM2);
03961 if(cLT+pRT==3) cRM-=baryM;
03962 RM=std::min(pRM,cRM);
03963 }
03964 }
03965 G4int dir=0;
03966 G4double sM=0.;
03967 if( LM && LM > RM )
03968 {
03969 dir= 1;
03970 sM=LM;
03971 }
03972 else if(RM)
03973 {
03974 dir=-1;
03975 sM=RM;
03976 }
03977 if(sM > maxM)
03978 {
03979 sst=pst;
03980 maxM=sM;
03981 sDir=dir;
03982 }
03983 }
03984 if(sDir)
03985 {
03986 G4QParton* pLeft=(*sst)->GetLeftParton();
03987 G4QParton* pRight=(*sst)->GetRightParton();
03988 G4QParton* swap=pLeft;
03989 if(sDir>0)
03990 {
03991 (*sst)->SetLeftParton(cLeft);
03992 (*ist)->SetLeftParton(swap);
03993 }
03994 else
03995 {
03996 swap=pRight;
03997 (*sst)->SetRightParton(cRight);
03998 (*ist)->SetRightParton(swap);
03999 }
04000 }
04001 #ifdef debug
04002 else G4cout<<"***G4QIonIonCollision::SwapPartons:**Failed**,cLPDG="<<cLPDG<<",cRPDG="
04003 <<cRPDG<<",-->cM2="<<cSM2<<G4endl;
04004 #endif
04005
04006 }
04007 }
04008 }