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 "G4QFragmentation.hh"
00050 #include "G4PhysicalConstants.hh"
00051 #include "G4SystemOfUnits.hh"
00052
00053
00054
00055
00056 G4int G4QFragmentation::nCutMax=27;
00057
00058
00059 G4double G4QFragmentation::stringTension=2*GeV/fermi;
00060
00061
00062
00063 G4double G4QFragmentation::tubeDensity =.5/fermi;
00064
00065
00066 G4double G4QFragmentation::widthOfPtSquare=-GeV*GeV;
00067
00068 G4QFragmentation::G4QFragmentation(const G4QNucleus &aNucleus, const G4QHadron &aPrimary)
00069 {
00070 static const G4double mProt= G4QPDGCode(2212).GetMass();
00071 static const G4double mProt2= mProt*mProt;
00072 static const G4double mPi0= G4QPDGCode(111).GetMass();
00073 static const G4double thresh= 1000;
00074 theWorld= G4QCHIPSWorld::Get();
00075 theQuasiElastic=G4QuasiFreeRatios::GetPointer();
00076 theDiffraction=G4QDiffractionRatio::GetPointer();
00077 theResult = new G4QHadronVector;
00078 G4bool stringsInitted=false;
00079 G4QHadron aProjectile = aPrimary;
00080 G4LorentzRotation toZ;
00081 G4LorentzVector proj4M=aProjectile.Get4Momentum();
00082 #ifdef edebug
00083 G4LorentzVector targ4M=aProjectile.Get4Momentum();
00084 G4double tgMass=aNucleus.GetGSMass();
00085 G4double cM=0.;
00086 G4double cs=(proj4M+targ4M).mag2();
00087 if(cs > 0.) cM=std::sqrt(cs);
00088 G4cout<<"G4QFragmentation::Construct: *Called*, p4M="<<proj4M<<", A="<<aNucleus<<tgMass
00089 <<",M="<<cM<<",s="<<cs<<",t4M="<<targ4M<<G4endl;
00090 #endif
00091 G4int tZ=aNucleus.GetZ();
00092 G4int tN=aNucleus.GetN();
00093 G4int tPDG=90000000+tZ*1000+tN;
00094 toZ.rotateZ(-proj4M.phi());
00095 toZ.rotateY(-proj4M.theta());
00096 G4LorentzVector zProj4M=toZ*proj4M;
00097 aProjectile.Set4Momentum(zProj4M);
00098 #ifdef edebug
00099 G4int totChg=aProjectile.GetCharge()+tZ;
00100 G4int totBaN=aProjectile.GetBaryonNumber()+tZ+tN;
00101 G4LorentzVector tgLS4M(0.,0.,0.,tgMass);
00102 G4LorentzVector totLS4M=proj4M+tgLS4M;
00103 G4LorentzVector totZLS4M=zProj4M+tgLS4M;
00104 G4cout<<"-EMC-G4QFragmentation::Construct:tLS4M="<<totLS4M<<",tZLS4M="<<totZLS4M<<G4endl;
00105
00106 #endif
00107 G4LorentzRotation toLab(toZ.inverse());
00108 G4int pPDG=aProjectile.GetPDGCode();
00109 G4double projM=aProjectile.GetMass();
00110 G4QInteractionVector theInteractions;
00111 G4QPartonPairVector thePartonPairs;
00112 G4int ModelMode=SOFT;
00113 theNucleus=G4QNucleus(tZ,tN);
00114 theNucleus.InitByPDG(tPDG);
00115 #ifdef debug
00116 G4cout<<"G4QFragmentation::Construct: Nucleus4Mom="<<theNucleus.Get4Momentum()<<G4endl;
00117 #endif
00118 theNucleus.Init3D();
00119
00120 std::pair<G4double,G4double> ratios=std::make_pair(0.,0.);
00121 G4int apPDG=std::abs(pPDG);
00122 if(apPDG>99)
00123 {
00124 G4double pMom=proj4M.vect().mag();
00125 ratios = theQuasiElastic->GetRatios(pMom, pPDG, tZ, tN);
00126 G4double qeRat = ratios.first*ratios.second;
00127 G4double difRat= theDiffraction->GetRatio(pMom, pPDG, tZ, tN);
00128
00129 if(qeRat<1. && proj4M.e()>thresh) difRat /= 1.-qeRat;
00130 else difRat=0.;
00131 difRat += qeRat;
00132 G4double rnd=G4UniformRand();
00133 if(rnd < qeRat)
00134 {
00135 theNucleus.StartLoop();
00136 G4QHadron* pNucleon = theNucleus.GetNextNucleon();
00137 G4LorentzVector pN4M=pNucleon->Get4Momentum();
00138 G4int pNPDG=pNucleon->GetPDGCode();
00139 #ifdef debug
00140 G4cout<<">QE>G4QFragmentation::Construct:TryQE,R="<<ratios.first*ratios.second<<",N4M="
00141 <<pN4M<<",NPDG="<<pNPDG<<G4endl;
00142 #endif
00143 std::pair<G4LorentzVector,G4LorentzVector> QEout=theQuasiElastic->Scatter(pNPDG,pN4M,
00144 pPDG,proj4M);
00145 G4bool CoulB = true;
00146 G4double CB=theNucleus.CoulombBarrier(1, 1);
00147 if( (pNPDG==2212 && QEout.first.e()-mProt < CB) ||
00148 (pPDG==2212 && QEout.second.e()-mProt < CB) ) CoulB = false;
00149 if(QEout.first.e() > 0 && CoulB)
00150 {
00151 G4QHadron* qfN = new G4QHadron(pNPDG,QEout.first);
00152 theResult->push_back(qfN);
00153 G4QHadron* qfP = new G4QHadron(pPDG,QEout.second);
00154 theResult->push_back(qfP);
00155 theNucleus.SubtractNucleon(pNucleon);
00156 G4LorentzVector r4M=theNucleus.Get4Momentum();
00157 G4int rPDG=theNucleus.GetPDG();
00158 G4QHadron* resNuc = new G4QHadron(rPDG,r4M);
00159 theResult->push_back(resNuc);
00160 #ifdef debug
00161 G4cout<<"-->QE-->G4QFragmentation::Construct:QEDone, N4M="<<QEout.first<<", p4M="
00162 <<QEout.second<<G4endl;
00163 #endif
00164 return;
00165 }
00166 }
00167 else if(rnd < difRat)
00168 {
00169 #ifdef debug
00170 G4cout<<"-->Dif-->G4QFragmentation::Construct: qe="<<qeRat<<", dif="<<difRat-qeRat
00171 <<",P="<<proj4M.vect().mag()<<", tZ="<<tZ<<", tN="<<tN<<G4endl;
00172 #endif
00173 G4QHadronVector* out=0;
00174
00175
00176
00177 out=theDiffraction->ProjFragment(pPDG, proj4M, tZ, tN);
00178 G4int nSec=out->size();
00179 if(nSec>1) for(G4int i=0; i<nSec; i++) theResult->push_back((*out)[i]);
00180 else if(nSec>0)
00181 {
00182 #ifdef debug
00183 G4cout<<"-Warning-G4QFragmentation::Construct: 1 secondary in Diffractionp 4M="
00184 <<proj4M<<", s4M="<<(*out)[0]->Get4Momentum()<<G4endl;
00185 #endif
00186 delete (*out)[0];
00187 }
00188 out->clear();
00189 delete out;
00190 if(nSec>1) return;
00191 }
00192 }
00193 #ifdef edebug
00194 G4LorentzVector sum1=theNucleus.GetNucleons4Momentum();
00195 G4cout<<"-EMC-G4QFragmentation::Construct: Nuc4M="<<sum1<<G4endl;
00196 #endif
00197 G4ThreeVector theCurrentVelocity(0.,0.,0.);
00198
00199 G4double nCons = 1;
00200 G4int projAbsB=std::abs(aProjectile.GetBaryonNumber());
00201 if(projAbsB>1) nCons = projAbsB;
00202
00203 proj4M = aProjectile.Get4Momentum();
00204 G4double pz_per_projectile = proj4M.pz()/nCons;
00205
00206 G4double e_per_projectile = proj4M.e()/nCons+mProt;
00207 G4double vz = pz_per_projectile/e_per_projectile;
00208 #ifdef debug
00209 G4cout<<"G4QFragmentation::Construct: Projectile4M="<<proj4M<<", Vz="<<vz<<", nC="
00210 <<nCons<<", pE="<<e_per_projectile<<G4endl;
00211 #endif
00212 theCurrentVelocity.setZ(vz);
00213 theNucleus.DoLorentzBoost(-theCurrentVelocity);
00214 #ifdef edebug
00215 G4LorentzVector sum2=theNucleus.GetNucleons4Momentum();
00216 G4cout<<"-EMC-G4QFragmentation::Construct: AfterBoost, v="<<vz<<", Nuc4M="<<sum2<<G4endl;
00217 #endif
00218 G4LorentzVector cmProjMom = proj4M;
00219 cmProjMom.boost(-theCurrentVelocity);
00220 G4double kE=cmProjMom.e()-projM;
00221 #ifdef debug
00222 G4cout<<"G4QFragmentation::Construct: kE="<<kE<<G4endl;
00223 #endif
00224 G4int maxCt=1;
00225 if(kE > 720.) maxCt=static_cast<int>(std::log(kE/270.));
00226
00227
00228 G4int maxCuts=std::min( 7 , std::max(1, maxCt) );
00229 #ifdef debug
00230 G4cout<<"G4QFragmentation::Construct: Proj4MInCM="<<cmProjMom<<", pPDG="<<pPDG<<G4endl;
00231 #endif
00232
00233
00234
00235 G4QHadron* cmProjectile = new G4QHadron(pPDG,cmProjMom);
00236
00237 G4QProbability theProbability(pPDG);
00238 G4double outerRadius = theNucleus.GetOuterRadius();
00239 #ifdef debug
00240 G4cout<<"G4QFrag::Constr:OutR="<<outerRadius<<",mC="<<maxCuts<<",A="<<theNucleus<<G4endl;
00241 #endif
00242 G4QHadron* pNucleon=0;
00243
00244 G4int theNA=theNucleus.GetA();
00245 G4LorentzVector pNuc4M=theNucleus.Get4Momentum()/theNA;
00246 G4double s_value = (cmProjMom + pNuc4M).mag2();
00247 G4double ThresholdMass = projM + theNucleus.GetGSMass()/theNA;
00248 #ifdef debug
00249 G4cout<<"G4QFrag::Construc: p4M="<<cmProjMom<<", tgN4M="<<pNuc4M<<", s="<<s_value<<", ThreshM="
00250 <<ThresholdMass<<G4endl;
00251 #endif
00252 ModelMode = SOFT;
00253 if (s_value < 0.)
00254 {
00255 G4cerr<<"***G4QFragmentation::Construct: s="<<s_value<<", pN4M="<<pNuc4M<<G4endl;
00256 G4Exception("G4QFragmentation::Construct:","72",FatalException,"LowEnergy(NegativeS)");
00257 }
00258 else if(s_value < mProt2)
00259 {
00260 theNucleus.StartLoop();
00261 G4QHadron* aTarget=0;
00262 G4QHadron* pProjectile=0;
00263 G4QHadron* bNucleon=0;
00264 G4double maxS=0.;
00265 while((bNucleon=theNucleus.GetNextNucleon()))
00266 {
00267 G4LorentzVector cp4M=bNucleon->Get4Momentum();
00268 G4double cs=(cmProjMom + cp4M).mag2();
00269 if(cs > maxS)
00270 {
00271 maxS=cs;
00272 pNucleon=bNucleon;
00273 }
00274 }
00275 aTarget = new G4QHadron(*pNucleon);
00276 pProjectile =cmProjectile;
00277 theNucleus.DoLorentzBoost(theCurrentVelocity);
00278 theNucleus.SubtractNucleon(pNucleon);
00279 theNucleus.DoLorentzBoost(-theCurrentVelocity);
00280 G4QContent QQC=aTarget->GetQC()+pProjectile->GetQC();
00281 G4LorentzVector Q4M=aTarget->Get4Momentum()+pProjectile->Get4Momentum();
00282 delete aTarget;
00283 delete pProjectile;
00284
00285
00286
00287
00288
00289
00290
00291
00292
00293
00294
00295
00296
00297 Q4M.boost(theCurrentVelocity);
00298 Q4M=toLab*Q4M;
00299 G4Quasmon* stringQuasmon = new G4Quasmon(QQC, Q4M);
00300 theQuasmons.push_back(stringQuasmon);
00301 theNucleus.DoLorentzBoost(theCurrentVelocity);
00302 theNucleus.DoLorentzRotation(toLab);
00303 return;
00304 }
00305 if (s_value < sqr(ThresholdMass))
00306 {
00307 #ifdef debug
00308 G4cout<<"G4QFragmentation::Construct:*OnlyDiffraction*ThM="<<ThresholdMass<<">sqrt(s)="
00309 <<std::sqrt(s_value)<<" -> only Diffraction is possible"<<G4endl;
00310 #endif
00311 ModelMode = DIFFRACTIVE;
00312 }
00313
00314 #ifdef debug
00315 G4cout<<"G4QFragmentation::Construct: theIntSize="<<theInteractions.size()<<G4endl;
00316 #endif
00317 std::for_each(theInteractions.begin(),theInteractions.end(), DeleteQInteraction());
00318 theInteractions.clear();
00319 G4int totalCuts = 0;
00320 G4int attCnt=0;
00321
00322 G4int maxAtt=27;
00323 G4double prEn=proj4M.e();
00324 G4int proB=aProjectile.GetBaryonNumber();
00325 if (proB>0) prEn-=aProjectile.GetMass();
00326 else if(proB<0) prEn+=mProt;
00327 #ifdef debug
00328 G4double impactUsed = 0.;
00329 G4cout<<"G4QFragmentation::Construct: estimated energy, prEn="<<prEn<<G4endl;
00330 #endif
00331 while(!theInteractions.size() && ++attCnt < maxAtt)
00332 {
00333 #ifdef debug
00334 G4cout<<"G4QFragmentation::Construct: *EnterTheInteractionLOOP*, att#"<<attCnt<<G4endl;
00335 #endif
00336
00337 std::pair<G4double, G4double> theImpactParameter;
00338 theImpactParameter = theNucleus.ChooseImpactXandY(outerRadius);
00339 G4double impactX = theImpactParameter.first;
00340 G4double impactY = theImpactParameter.second;
00341 #ifdef debug
00342 G4cout<<"G4QFragmentation::Construct: Impact Par X="<<impactX<<", Y="<<impactY<<G4endl;
00343 #endif
00344 G4double impar=std::sqrt(impactX*impactX+impactY*impactY);
00345 G4int nA=theNucleus.GetA();
00346 G4double eflen=theNucleus.GetThickness(impar);
00347 maxEn=eflen*stringTension;
00348 maxNuc=static_cast<int>(eflen*tubeDensity+0.5);
00349 #ifdef debug
00350 G4cout<<"G4QFragment::Construct: pE="<<prEn<<" <? mE="<<maxEn<<", mN="<<maxNuc<<G4endl;
00351 #endif
00352 if(prEn < maxEn)
00353 {
00354 theNucleus.StartLoop();
00355 pNucleon=theNucleus.GetNextNucleon();
00356 G4QHadron* aTarget = new G4QHadron(*pNucleon);
00357 G4QInteraction* anInteraction = new G4QInteraction(cmProjectile);
00358 anInteraction->SetTarget(aTarget);
00359 anInteraction->SetNumberOfDINRCollisions(1);
00360 theInteractions.push_back(anInteraction);
00361 theNucleus.DoLorentzBoost(theCurrentVelocity);
00362 theNucleus.SubtractNucleon(pNucleon);
00363 theNucleus.DoLorentzBoost(-theCurrentVelocity);
00364 #ifdef debug
00365 G4cout<<"G4QFragmentation::Construct: DINR interaction is created"<<G4endl;
00366 #endif
00367 break;
00368 }
00369
00370 theNucleus.StartLoop();
00371 G4int nucleonCount = 0;
00372 #ifdef debug
00373 G4cout<<"G4QFragment::Construct:BeforeWhileOveNuc, A="<<nA<<",p4M="<<cmProjMom<<G4endl;
00374 #endif
00375 while( (pNucleon=theNucleus.GetNextNucleon()) && nucleonCount<nA && totalCuts<maxCuts)
00376 {
00377 ++nucleonCount;
00378
00379 s_value = (cmProjMom + pNucleon->Get4Momentum()).mag2();
00380 #ifdef debug
00381 G4cout<<"G4QFrag::Constr:N# "<<nucleonCount<<", s="<<s_value<<", tgN4M="
00382 <<pNucleon->Get4Momentum()<<G4endl;
00383 #endif
00384 if(s_value<=10000.)
00385 {
00386 #ifdef debug
00387 G4cout<<"G4QFragmentation::Construct: SKIP, s<.01 GeV^2, p4M="<<cmProjMom
00388 <<",t4M="<<pNucleon->Get4Momentum()<<G4endl;
00389 #endif
00390 continue;
00391 }
00392 #ifdef sdebug
00393 G4cout<<"G4QFragmentation::Construct:LOOPovNuc,nC="<<nucleonCount<<", s="<<s_value<<G4endl;
00394 G4cout<<"G4QFragmentation::Construct:LOOPovNuc, R="<<pNucleon->GetPosition()<<G4endl;
00395 #endif
00396 G4double Distance2 = sqr(impactX - pNucleon->GetPosition().x()) +
00397 sqr(impactY - pNucleon->GetPosition().y());
00398 #ifdef sdebug
00399 G4cout<<"G4QFragmentation::Construct: s="<<s_value<<", D2="<<Distance2<<G4endl;
00400 #endif
00401 G4double Probability = theProbability.GetPomInelProbability(s_value, Distance2);
00402
00403 #ifdef sdebug
00404 G4cout<<"G4QFragmentation::Construct: Probubility="<<Probability<<G4endl;
00405 #endif
00406 G4double rndNumber = G4UniformRand();
00407
00408 #ifdef sdebug
00409 G4cout<<"G4QFragmentation::Construct: NLOOP prob="<<Probability<<", rndm="<<rndNumber
00410 <<", d="<<std::sqrt(Distance2)<<G4endl;
00411 #endif
00412 if (Probability > rndNumber)
00413 {
00414 G4QHadron* aTarget = new G4QHadron(*pNucleon);
00415 #ifdef edebug
00416 G4cout<<"--->EMC-->G4QFragmentation::Construct: Target Nucleon is filled, 4M/PDG="
00417 <<aTarget->Get4Momentum()<<aTarget->GetPDGCode()<<G4endl;
00418 #endif
00419
00420 theNucleus.DoLorentzBoost(theCurrentVelocity);
00421 theNucleus.SubtractNucleon(pNucleon);
00422 theNucleus.DoLorentzBoost(-theCurrentVelocity);
00423 if((theProbability.GetPomDiffProbability(s_value,Distance2)/Probability >
00424 G4UniformRand() && ModelMode==SOFT ) || ModelMode==DIFFRACTIVE)
00425 {
00426
00427 if(IsSingleDiffractive()) ExciteSingDiffParticipants(cmProjectile, aTarget);
00428 else ExciteDiffParticipants(cmProjectile, aTarget);
00429 G4QInteraction* anInteraction = new G4QInteraction(cmProjectile);
00430 anInteraction->SetTarget(aTarget);
00431 anInteraction->SetNumberOfDiffractiveCollisions(1);
00432 theInteractions.push_back(anInteraction);
00433
00434 totalCuts++;
00435 #ifdef debug
00436 G4cout<<"G4QFragmentation::Construct:NLOOP DiffInteract, tC="<<totalCuts<<G4endl;
00437 #endif
00438 }
00439 else
00440 {
00441
00442
00443 G4int nCut;
00444 G4double* running = new G4double[nCutMax];
00445 for(nCut = 0; nCut < nCutMax; nCut++)
00446 {
00447 running[nCut]= theProbability.GetCutPomProbability(s_value, Distance2, nCut+1);
00448 if(nCut) running[nCut] += running[nCut-1];
00449 }
00450 G4double random = running[nCutMax-1]*G4UniformRand();
00451 for(nCut = 0; nCut < nCutMax; nCut++) if(running[nCut] > random) break;
00452 delete [] running;
00453 #ifdef debug
00454 G4cout<<"G4QFragmentation::Construct: NLOOP-Soft Chosen nCut="<<nCut<<G4endl;
00455 #endif
00456
00457
00458
00459 aTarget->IncrementCollisionCount(nCut+1);
00460 cmProjectile->IncrementCollisionCount(nCut+1);
00461 G4QInteraction* anInteraction = new G4QInteraction(cmProjectile);
00462 anInteraction->SetTarget(aTarget);
00463 anInteraction->SetNumberOfSoftCollisions(nCut+1);
00464 theInteractions.push_back(anInteraction);
00465 totalCuts += nCut+1;
00466 #ifdef debug
00467 G4cout<<"G4QFragmentation::Construct:NLOOP SoftInteract, tC="<<totalCuts<<G4endl;
00468 impactUsed=Distance2;
00469 #endif
00470 }
00471 }
00472 }
00473
00474 #ifdef debug
00475 G4cout<<"G4QFragmentation::Construct: NUCLEONCOUNT="<<nucleonCount<<G4endl;
00476 #endif
00477 }
00478 G4int nInt=theInteractions.size();
00479 #ifdef debug
00480 G4cout<<"G4QFrag::Con:CUT="<<totalCuts<<",ImpPar="<<impactUsed<<",#ofInt="<<nInt<<G4endl;
00481 #endif
00482
00483 if(!nInt || (nInt==1 && theInteractions[0]->GetNumberOfDINRCollisions()==1))
00484 {
00485 G4QHadron* aTarget=0;
00486 G4QHadron* pProjectile=0;
00487 if(nInt)
00488 {
00489 aTarget=theInteractions[0]->GetTarget();
00490 pProjectile=theInteractions[0]->GetProjectile();
00491 delete theInteractions[0];
00492 theInteractions.clear();
00493 }
00494 else
00495 {
00496 theNucleus.StartLoop();
00497 pNucleon=theNucleus.GetNextNucleon();
00498 aTarget = new G4QHadron(*pNucleon);
00499 pProjectile =cmProjectile;
00500 theNucleus.DoLorentzBoost(theCurrentVelocity);
00501 theNucleus.SubtractNucleon(pNucleon);
00502 theNucleus.DoLorentzBoost(-theCurrentVelocity);
00503 }
00504 G4QContent QQC=aTarget->GetQC()+pProjectile->GetQC();
00505 G4LorentzVector Q4M=aTarget->Get4Momentum()+pProjectile->Get4Momentum();
00506 delete aTarget;
00507 delete pProjectile;
00508
00509
00510
00511
00512
00513
00514
00515
00516
00517
00518
00519
00520
00521 Q4M.boost(theCurrentVelocity);
00522 Q4M=toLab*Q4M;
00523 G4Quasmon* stringQuasmon = new G4Quasmon(QQC, Q4M);
00524 theQuasmons.push_back(stringQuasmon);
00525 theNucleus.DoLorentzBoost(theCurrentVelocity);
00526 theNucleus.DoLorentzRotation(toLab);
00527 return;
00528 }
00529
00530
00531
00532 #ifdef debug
00533 G4cout<<"G4QFragmentation::Construct: Before PartPairCreation nInt="<<nInt<<G4endl;
00534 #endif
00535 for(G4int i=0; i<nInt; i++)
00536 {
00537 theInteractions[i]->SplitHadrons();
00538 #ifdef edebug
00539 G4QHadron* projH=theInteractions[i]->GetProjectile();
00540 G4QHadron* targH=theInteractions[i]->GetTarget();
00541 G4LorentzVector pSP(0.,0.,0.,0.);
00542 G4LorentzVector tSP(0.,0.,0.,0.);
00543 std::list<G4QParton*> projCP=projH->GetColor();
00544 std::list<G4QParton*> projAC=projH->GetAntiColor();
00545 std::list<G4QParton*> targCP=targH->GetColor();
00546 std::list<G4QParton*> targAC=targH->GetAntiColor();
00547 std::list<G4QParton*>::iterator picp = projCP.begin();
00548 std::list<G4QParton*>::iterator pecp = projCP.end();
00549 std::list<G4QParton*>::iterator piac = projAC.begin();
00550 std::list<G4QParton*>::iterator peac = projAC.end();
00551 std::list<G4QParton*>::iterator ticp = targCP.begin();
00552 std::list<G4QParton*>::iterator tecp = targCP.end();
00553 std::list<G4QParton*>::iterator tiac = targAC.begin();
00554 std::list<G4QParton*>::iterator teac = targAC.end();
00555 for(; picp!=pecp&& piac!=peac&& ticp!=tecp&& tiac!=teac; ++picp,++piac,++ticp,++tiac)
00556 {
00557 pSP+=(*picp)->Get4Momentum();
00558 pSP+=(*piac)->Get4Momentum();
00559 tSP+=(*ticp)->Get4Momentum();
00560 tSP+=(*tiac)->Get4Momentum();
00561 }
00562 G4cout<<"-EMC-G4QFragmentation::Construct: Interaction#"<<i<<",dP4M="
00563 <<projH->Get4Momentum()-pSP<<",dT4M="<<targH->Get4Momentum()-tSP<<G4endl;
00564 #endif
00565 }
00566
00567
00568
00569 G4QInteractionVector::iterator it;
00570 #ifdef debug
00571 G4cout<<"G4QFragmentation::Construct: Creation ofSoftCollisionPartonPair STARTS"<<G4endl;
00572 #endif
00573 G4bool rep=true;
00574 while(rep && theInteractions.size())
00575 {
00576 for(it = theInteractions.begin(); it != theInteractions.end(); ++it)
00577 {
00578 G4QInteraction* anIniteraction = *it;
00579 G4QPartonPair* aPair=0;
00580 G4int nSoftCollisions = anIniteraction->GetNumberOfSoftCollisions();
00581 #ifdef debug
00582 G4cout<<"G4QFragmentation::Construct: #0f SOFT collisions ="<<nSoftCollisions<<G4endl;
00583 #endif
00584 if (nSoftCollisions)
00585 {
00586 G4QHadron* pProjectile = anIniteraction->GetProjectile();
00587 G4QHadron* pTarget = anIniteraction->GetTarget();
00588 for (G4int j = 0; j < nSoftCollisions; j++)
00589 {
00590 aPair = new G4QPartonPair(pTarget->GetNextParton(),
00591 pProjectile->GetNextAntiParton(),
00592 G4QPartonPair::SOFT, G4QPartonPair::TARGET);
00593 thePartonPairs.push_back(aPair);
00594 aPair = new G4QPartonPair(pProjectile->GetNextParton(),
00595 pTarget->GetNextAntiParton(),
00596 G4QPartonPair::SOFT, G4QPartonPair::PROJECTILE);
00597 thePartonPairs.push_back(aPair);
00598 #ifdef debug
00599 G4cout<<"--->G4QFragmentation::Construct: SOFT, 2 parton pairs are filled"<<G4endl;
00600 #endif
00601 }
00602 delete *it;
00603 it=theInteractions.erase(it);
00604 if( it != theInteractions.begin() )
00605 {
00606 it--;
00607 rep=false;
00608 #ifdef debug
00609 G4cout<<"G4QFragmentation::Construct: *** Decremented ***"<<G4endl;
00610 #endif
00611 }
00612 else
00613 {
00614 rep=true;
00615 #ifdef debug
00616 G4cout<<"G4QFragmentation::Construct: *** IT Begin ***"<<G4endl;
00617 #endif
00618 break;
00619 }
00620 }
00621 else rep=false;
00622 #ifdef debug
00623 G4cout<<"G4QFragmentation::Construct: #0fSC="<<nSoftCollisions<<", r="<<rep<<G4endl;
00624 #endif
00625 }
00626 #ifdef debug
00627 G4cout<<"G4QFragmentation::Construct: *** IT While *** , r="<<rep<<G4endl;
00628 #endif
00629 }
00630 #ifdef debug
00631 G4cout<<"G4QFragmentation::Construct: -> Parton pairs for SOFT strings are made"<<G4endl;
00632 #endif
00633
00634
00635
00636 for(unsigned i = 0; i < theInteractions.size(); i++)
00637 {
00638
00639 G4QInteraction* anIniteraction = theInteractions[i];
00640 G4QPartonPair* aPartonPair;
00641 #ifdef debug
00642 G4cout<<"G4QFragmentation::Construct: CreationOfDiffractivePartonPairs, i="<<i<<G4endl;
00643 #endif
00644
00645 G4QHadron* pProjectile = anIniteraction->GetProjectile();
00646 G4QParton* aParton = pProjectile->GetNextParton();
00647 if (aParton)
00648 {
00649 aPartonPair = new G4QPartonPair(aParton, pProjectile->GetNextAntiParton(),
00650 G4QPartonPair::DIFFRACTIVE,
00651 G4QPartonPair::PROJECTILE);
00652 thePartonPairs.push_back(aPartonPair);
00653 #ifdef debug
00654 G4cout<<"G4QFragmentation::Construct: proj Diffractive PartonPair is filled"<<G4endl;
00655 #endif
00656 }
00657
00658 G4QHadron* aTarget = anIniteraction->GetTarget();
00659 aParton = aTarget->GetNextParton();
00660 if (aParton)
00661 {
00662 aPartonPair = new G4QPartonPair(aParton, aTarget->GetNextAntiParton(),
00663 G4QPartonPair::DIFFRACTIVE, G4QPartonPair::TARGET);
00664 thePartonPairs.push_back(aPartonPair);
00665 #ifdef debug
00666 G4cout<<"G4QFragmentation::Constr: target Diffractive PartonPair is filled"<<G4endl;
00667 #endif
00668 }
00669 }
00670 #ifdef debug
00671 G4cout<<"G4QFragmentation::Construct: DiffractivePartonPairs are created"<<G4endl;
00672 #endif
00673
00674
00675
00676 std::for_each(theInteractions.begin(),theInteractions.end(), DeleteQInteraction());
00677 theInteractions.clear();
00678 delete cmProjectile;
00679 #ifdef debug
00680 G4cout<<"G4QFragmentation::Construct: Temporary objects are cleaned up"<<G4endl;
00681 #endif
00682
00683 theNucleus.DoLorentzBoost(theCurrentVelocity);
00684
00685 #ifdef debug
00686 G4cout<<"--------->>G4QFragmentation::Construct: ------->> Strings are created "<<G4endl;
00687 #endif
00688 G4QPartonPair* aPair;
00689 G4QString* aString=0;
00690 while(thePartonPairs.size())
00691 {
00692 aPair = thePartonPairs.back();
00693 thePartonPairs.pop_back();
00694 #ifdef debug
00695 G4cout<<"G4QFragmentation::Construct: StringType="<<aPair->GetCollisionType()<<G4endl;
00696 #endif
00697 aString= new G4QString(aPair);
00698 #ifdef debug
00699 G4cout<<"G4QFragmentation::Construct:NewString4M="<<aString->Get4Momentum()<<G4endl;
00700 #endif
00701 aString->Boost(theCurrentVelocity);
00702 strings.push_back(aString);
00703 stringsInitted=true;
00704 delete aPair;
00705 }
00706 #ifdef edebug
00707 G4LorentzVector sum=theNucleus.Get4Momentum();
00708 G4int rChg=totChg-theNucleus.GetZ();
00709 G4int rBaN=totBaN-theNucleus.GetA();
00710 G4int nStrings=strings.size();
00711 G4cout<<"-EMC-G4QFragmentation::Construct:#ofString="<<nStrings<<",tNuc4M="<<sum<<G4endl;
00712 for(G4int i=0; i<nStrings; i++)
00713 {
00714 G4QString* prString=strings[i];
00715 G4LorentzVector strI4M=prString->Get4Momentum();
00716 sum+=strI4M;
00717 G4int sChg=prString->GetCharge();
00718 G4int sBaN=prString->GetBaryonNumber();
00719 G4int LPDG=prString->GetLeftParton()->GetPDGCode();
00720 G4int RPDG=prString->GetRightParton()->GetPDGCode();
00721 G4QContent LQC =prString->GetLeftParton()->GetQC();
00722 G4QContent RQC =prString->GetRightParton()->GetQC();
00723 rChg-=sChg;
00724 rBaN-=sBaN;
00725 G4cout<<"-EMC-G4QFragmentation::Construct: String#"<<i<<", 4M="<<strI4M<<",LPDG="<<LPDG
00726 <<LQC<<",RPDG="<<RPDG<<RQC<<", Ch="<<sChg<<", BN="<<sBaN<<G4endl;
00727 }
00728 G4cout<<"-EMC-G4QFragm::Constr: r4M="<<sum-totZLS4M<<",rC="<<rChg<<",rB="<<rBaN<<G4endl;
00729 #endif
00730 if(!stringsInitted)
00731 {
00732 G4cerr<<"******G4QFragmentation::Construct:***** No strings are created *****"<<G4endl;
00733 G4Exception("G4QFragmentation::Construct:","72",FatalException,"NoStrings're created");
00734 }
00735 #ifdef debug
00736 G4cout<<"G4QFragmentation::Constr: BeforeRotation, #0fStrings="<<strings.size()<<G4endl;
00737 #endif
00738
00739
00740
00741 for(unsigned astring=0; astring < strings.size(); astring++)
00742 strings[astring]->LorentzRotate(toLab);
00743 theNucleus.DoLorentzRotation(toLab);
00744
00745 #ifdef edebug
00746 G4LorentzVector sm=theNucleus.Get4Momentum();
00747 G4int rCg=totChg-theNucleus.GetZ();
00748 G4int rBC=totBaN-theNucleus.GetA();
00749 G4int nStrs=strings.size();
00750 G4cout<<"-EMCLS-G4QFragmentation::Constr: #ofS="<<nStrings<<",tNuc4M(E=M)="<<sum<<G4endl;
00751 for(G4int i=0; i<nStrs; i++)
00752 {
00753 G4LorentzVector strI4M=strings[i]->Get4Momentum();
00754 sm+=strI4M;
00755 G4int sChg=strings[i]->GetCharge();
00756 rCg-=sChg;
00757 G4int sBaN=strings[i]->GetBaryonNumber();
00758 rBC-=sBaN;
00759 G4cout<<"-EMCLS-G4QFragm::Construct:String#"<<i<<",4M="<<strI4M<<strI4M.m()<<",Charge="
00760 <<sChg<<",BaryN="<<sBaN<<G4endl;
00761 }
00762 G4cout<<"-EMCLS-...G4QFragm::Constr:r4M="<<sm-totLS4M<<",rC="<<rCg<<",rB="<<rBC<<G4endl;
00763 #endif
00764
00765
00766
00767 SwapPartons();
00768 #ifdef edebug
00769 sm=theNucleus.Get4Momentum();
00770 rCg=totChg-theNucleus.GetZ();
00771 rBC=totBaN-theNucleus.GetA();
00772 nStrs=strings.size();
00773 G4cout<<"-EMCLS-G4QFrag::Constr:AfterSwap #ofS="<<nStrings<<",tNuc4M(E=M)="<<sum<<G4endl;
00774 for(G4int i=0; i<nStrs; i++)
00775 {
00776 G4LorentzVector strI4M=strings[i]->Get4Momentum();
00777 sm+=strI4M;
00778 G4int sChg=strings[i]->GetCharge();
00779 rCg-=sChg;
00780 G4int sBaN=strings[i]->GetBaryonNumber();
00781 rBC-=sBaN;
00782 G4cout<<"-EMCLS-G4QFragm::Construct:String#"<<i<<",4M="<<strI4M<<strI4M.m()<<",Charge="
00783 <<sChg<<",BaryN="<<sBaN<<G4endl;
00784 }
00785 G4cout<<"-EMCLS-...G4QFragm::Constr:r4M="<<sm-totLS4M<<",rC="<<rCg<<",rB="<<rBC<<G4endl;
00786 #endif
00787
00788
00789
00790 G4int problem=0;
00791 G4QStringVector::iterator ist;
00792 G4bool con=true;
00793 while(con && strings.size())
00794 {
00795 for(ist = strings.begin(); ist < strings.end(); ++ist)
00796 {
00797 G4bool bad=true;
00798 G4LorentzVector cS4M=(*ist)->Get4Momentum();
00799 G4double cSM2=cS4M.m2();
00800 G4QParton* cLeft=(*ist)->GetLeftParton();
00801 G4QParton* cRight=(*ist)->GetRightParton();
00802 G4int cLT=cLeft->GetType();
00803 G4int cRT=cRight->GetType();
00804 G4int cLPDG=cLeft->GetPDGCode();
00805 G4int cRPDG=cRight->GetPDGCode();
00806 G4int aLPDG=0;
00807 G4int aRPDG=0;
00808 if (cLPDG > 7) aLPDG= cLPDG/100;
00809 else if(cLPDG <-7) aLPDG=(-cLPDG)/100;
00810 if (cRPDG > 7) aRPDG= cRPDG/100;
00811 else if(cRPDG <-7) aRPDG=(-cRPDG)/100;
00812 G4int L1=0;
00813 G4int L2=0;
00814 if(aLPDG)
00815 {
00816 L1=aLPDG/10;
00817 L2=aLPDG%10;
00818 }
00819 G4int R1=0;
00820 G4int R2=0;
00821 if(aRPDG)
00822 {
00823 R1=aRPDG/10;
00824 R2=aRPDG%10;
00825 }
00826 G4double cSM=cSM2;
00827 if(cSM2>0.) cSM=std::sqrt(cSM2);
00828 #ifdef debug
00829 G4cout<<"G4QFrag::Constr:NeedsFusion? cLPDG="<<cLPDG<<",cRPDG="<<cRPDG<<",cM(cM2If<0)="
00830 <<cSM<<",c4M"<<cS4M<<G4endl;
00831 #endif
00832 if(cSM>0.)
00833 {
00834 G4bool single=true;
00835 G4double miM=0.;
00836 if(cLT==2 && cRT==2)
00837 {
00838 if(L1!=R1 && L1!=R2 && L2!=R1 && L2!=R2)
00839 {
00840 single=false;
00841 G4QPDGCode tmp;
00842 std::pair<G4int,G4int> paB=tmp.MakeTwoBaryons(L1, L2, R1, R2);
00843 miM=G4QPDGCode(paB.first).GetMass()+G4QPDGCode(paB.second).GetMass();
00844 }
00845 }
00846 if(single) miM=G4QPDGCode((*ist)->GetQC().GetSPDGCode()).GetMass() + mPi0;
00847
00848 #ifdef debug
00849 G4cout<<"G4QFrag::Const:*IsItGood? realM="<<std::sqrt(cSM2)<<" > GSM="<<miM<<G4endl;
00850 #endif
00851 if(std::sqrt(cSM2) > miM) bad=false;
00852 }
00853 if(bad)
00854 {
00855 #ifdef debug
00856 G4cout<<"G4QFrag::Const:TryFuse,L1="<<L1<<",L2="<<L2<<",R1="<<R1<<",R2="<<R2<<G4endl;
00857 #endif
00858 G4int cST=cLT+cRT;
00859 G4double excess=-DBL_MAX;
00860 G4double maxiM2=-DBL_MAX;
00861 G4QStringVector::iterator sst;
00862 G4QStringVector::iterator pst;
00863 G4int sLPDG=0;
00864 G4int sRPDG=0;
00865 G4int sOrd=0;
00866 G4bool minC=true;
00867 if(cSM2>0.) minC=false;
00868 for(pst = strings.begin(); pst < strings.end(); pst++) if(pst != ist)
00869 {
00870 G4LorentzVector sS4M=(*pst)->Get4Momentum();
00871 G4LorentzVector pS4M=sS4M+cS4M;
00872 G4int nLPDG=0;
00873 G4int nRPDG=0;
00874 G4double pExcess=-DBL_MAX;
00875 G4double pSM2=pS4M.m2();
00876 #ifdef debug
00877 G4cout<<"->G4QFragm::Construct: sum4M="<<pS4M<<",M2="<<pSM2<<",p4M="<<sS4M<<G4endl;
00878 #endif
00879
00880
00881 G4QParton* pLeft=(*pst)->GetLeftParton();
00882 G4QParton* pRight=(*pst)->GetRightParton();
00883 G4int pLT=pLeft->GetType();
00884 G4int pRT=pRight->GetType();
00885 G4int pLPDG=pLeft->GetPDGCode();
00886 G4int pRPDG=pRight->GetPDGCode();
00887 G4int LPDG=0;
00888 G4int RPDG=0;
00889 if (pLPDG > 7) LPDG= pLPDG/100;
00890 else if(pLPDG <-7) LPDG=(-pLPDG)/100;
00891 if (pRPDG > 7) RPDG= pRPDG/100;
00892 else if(pRPDG <-7) RPDG=(-pRPDG)/100;
00893 G4int pL1=0;
00894 G4int pL2=0;
00895 if(LPDG)
00896 {
00897 pL1=LPDG/10;
00898 pL2=LPDG%10;
00899 }
00900 G4int pR1=0;
00901 G4int pR2=0;
00902 if(RPDG)
00903 {
00904 pR1=RPDG/10;
00905 pR2=RPDG%10;
00906 }
00907 G4int pST=pLT+pRT;
00908 #ifdef debug
00909 G4cout<<"G4QFragm::Construct: Partner/w pLPDG="<<pLPDG<<", pRPDG="<<pRPDG<<", pM2="
00910 <<pSM2<<G4endl;
00911 #endif
00912
00913 G4int tf=0;
00914 G4int af=0;
00915 if (cST==2 && pST==2)
00916 {
00917 tf=1;
00918 af=1;
00919 }
00920 else if(cST==2 && pST==3)
00921 {
00922 tf=2;
00923 if (pLPDG > 7 &&
00924 ( (cLPDG<0 && (-cLPDG==pL1 || -cLPDG==pL2 || -cLPDG==pRPDG) ) ||
00925 (cRPDG<0 && (-cRPDG==pL1 || -cRPDG==pL2 || -cRPDG==pRPDG) )
00926 )
00927 ) af=1;
00928 else if(pRPDG > 7 &&
00929 ( (cLPDG<0 && (-cLPDG==pR1 || -cLPDG==pR2 || -cLPDG==pLPDG) ) ||
00930 (cRPDG<0 && (-cRPDG==pR1 || -cRPDG==pR2 || -cRPDG==pLPDG) )
00931 )
00932 ) af=2;
00933 else if(pLPDG <-7 &&
00934 ( (cLPDG>0 && ( cLPDG==pL1 || cLPDG==pL2 || cLPDG==-pRPDG) ) ||
00935 (cRPDG>0 && ( cRPDG==pL1 || cRPDG==pL2 || cRPDG==-pRPDG) )
00936 )
00937 ) af=3;
00938 else if(pRPDG <-7 &&
00939 ( (cLPDG>0 && ( cLPDG==pR1 || cLPDG==pR2 || cLPDG==-pLPDG) ) ||
00940 (cRPDG>0 && ( cRPDG==pR1 || cRPDG==pR2 || cRPDG==-pLPDG) )
00941 )
00942 ) af=4;
00943 #ifdef debug
00944 else G4cout<<"G4QFragmentation::Construct:2(QaQ+QDiQ/aQaDiQ) Can't fuse"<<G4endl;
00945 #endif
00946 }
00947 else if(cST==3 && pST==2)
00948 {
00949 tf=3;
00950 if (cLPDG > 7 &&
00951 ( (pLPDG<0 && (-pLPDG==L1 || -pLPDG==L2 || -pLPDG==cRPDG) ) ||
00952 (pRPDG<0 && (-pRPDG==L1 || -pRPDG==L2 || -pRPDG==cRPDG) )
00953 )
00954 ) af=1;
00955 else if(cRPDG > 7 &&
00956 ( (pLPDG<0 && (-pLPDG==R1 || -pLPDG==R2 || -pLPDG==cLPDG) ) ||
00957 (pRPDG<0 && (-pRPDG==R1 || -pRPDG==R2 || -pRPDG==cLPDG) )
00958 )
00959 ) af=2;
00960 else if(cLPDG <-7 &&
00961 ( (pLPDG>0 && ( pLPDG==L1 || pLPDG==L2 || pLPDG==-cRPDG) ) ||
00962 (pRPDG>0 && ( pRPDG==L1 || pRPDG==L2 || pRPDG==-cRPDG) )
00963 )
00964 ) af=3;
00965 else if(cRPDG <-7 &&
00966 ( (pLPDG>0 && ( pLPDG==R1 || pLPDG==R2 || pLPDG==-cLPDG) ) ||
00967 (pRPDG>0 && ( pRPDG==R1 || pRPDG==R2 || pRPDG==-cLPDG) )
00968 )
00969 ) af=4;
00970 #ifdef debug
00971 else G4cout<<"G4QFragmentation::Construct:3(QDiQ/aQaDiQ+QaQ) Can't fuse"<<G4endl;
00972 #endif
00973 }
00974 else if(cST==2 && pST==4)
00975 {
00976 tf=4;
00977 if (pLPDG > 7)
00978 {
00979 if ( (-cLPDG==pL1 || -cLPDG==pL2) && (cRPDG==pR1 || cRPDG==pR2) ) af=1;
00980 else if( (-cRPDG==pL1 || -cRPDG==pL2) && (cLPDG==pR1 || cLPDG==pR2) ) af=2;
00981 }
00982 else if(pRPDG > 7)
00983 {
00984 if ( (-cRPDG==pR1 || -cRPDG==pR2) && (cLPDG==pL1 || cLPDG==pL2) ) af=3;
00985 else if( (-cLPDG==pR1 || -cLPDG==pR2) && (cRPDG==pL1 || cRPDG==pL2) ) af=4;
00986 }
00987 #ifdef debug
00988 else G4cout<<"-G4QFragmentation::Construct: 4 (QaQ+aQDiQDiQ) Can't fuse"<<G4endl;
00989 #endif
00990 }
00991 else if(cST==4 && pST==2)
00992 {
00993 tf=5;
00994 if (cLPDG > 7)
00995 {
00996 if ( (-pLPDG==L1 || -pLPDG==L2) && (pRPDG==R1 || pRPDG==R2) ) af=1;
00997 else if( (-pRPDG==L1 || -pRPDG==L2) && (pLPDG==R1 || pLPDG==R2) ) af=2;
00998 }
00999 else if(cRPDG > 7)
01000 {
01001 if ( (-pRPDG==R1 || -pRPDG==R2) && (pLPDG==L1 || pLPDG==L2) ) af=3;
01002 else if( (-pLPDG==R1 || -pLPDG==R2) && (pRPDG==L1 || pRPDG==L2) ) af=4;
01003 }
01004 #ifdef debug
01005 else G4cout<<"-G4QFragmentation::Construct: 5 (aQDiQDiQ+QaQ) Can't fuse"<<G4endl;
01006 #endif
01007 }
01008 else if(cST==3 && pST==3)
01009 {
01010 tf=6;
01011 if(pLPDG > 7)
01012 {
01013 if (cLPDG<-7 && (-cRPDG==pL1 || -cRPDG==pL2) && (pRPDG==L1 || pRPDG==L2))
01014 af=1;
01015 else if(cRPDG<-7 && (-cLPDG==pL1 || -cLPDG==pL2) && (pRPDG==R1 || pRPDG==R2))
01016 af=2;
01017 }
01018 else if(pRPDG > 7)
01019 {
01020 if (cLPDG<-7 && (-cRPDG==pR1 || -cRPDG==pR2) && (pLPDG==L1 || pLPDG==L2))
01021 af=3;
01022 else if(cRPDG<-7 && (-cLPDG==pR1 || -cLPDG==pR2) && (pLPDG==R1 || pLPDG==R2))
01023 af=4;
01024 }
01025 else if(cLPDG > 7)
01026 {
01027 if (pLPDG<-7 && (-pRPDG==L1 || -pRPDG==L2) && (cRPDG==pL1 || cRPDG==pL2))
01028 af=5;
01029 else if(pRPDG<-7 && (-pLPDG==L1 || -pLPDG==L2) && (cRPDG==pR1 || cRPDG==pR2))
01030 af=6;
01031 }
01032 else if(cRPDG > 7)
01033 {
01034 if (pLPDG<-7 && (-pRPDG==R1 || -pRPDG==R2) && (cLPDG==pL1 || cLPDG==pL2))
01035 af=7;
01036 else if(pRPDG<-7 && (-pLPDG==R1 || -pLPDG==R2) && (cLPDG==pR1 || cLPDG==pR2))
01037 af=8;
01038 }
01039 #ifdef debug
01040 else G4cout<<"-G4QFragmentation::Construct: 6 (QDiQ+aQaDiQ) Can't fuse"<<G4endl;
01041 #endif
01042 }
01043 #ifdef debug
01044 G4cout<<"G4QFragmentation::Const: ***Possibility***, tf="<<tf<<", af="<<af<<G4endl;
01045 #endif
01046 if(tf && af)
01047 {
01048
01049 G4int order=0;
01050 switch (tf)
01051 {
01052 case 1:
01053 if (cLPDG > 0 && pLPDG > 0)
01054 {
01055 order= 1;
01056 if (cLPDG > pLPDG) nLPDG=cLPDG*1000+pLPDG*100+1;
01057 else if(cLPDG < pLPDG) nLPDG=pLPDG*1000+cLPDG*100+1;
01058 else nLPDG=pLPDG*1000+cLPDG*100+3;
01059 if (cRPDG < pRPDG) nRPDG=cRPDG*1000+pRPDG*100-1;
01060 else if(cRPDG > pRPDG) nRPDG=pRPDG*1000+cRPDG*100-1;
01061 else nRPDG=pRPDG*1000+cRPDG*100-3;
01062 }
01063 else if(cLPDG < 0 && pLPDG < 0)
01064 {
01065 order= 1;
01066 if (cRPDG > pRPDG) nRPDG=cRPDG*1000+pRPDG*100+1;
01067 else if(cRPDG < pRPDG) nRPDG=pRPDG*1000+cRPDG*100+1;
01068 else nRPDG=pRPDG*1000+cRPDG*100+3;
01069 if (cLPDG < pLPDG) nLPDG=cLPDG*1000+pLPDG*100-1;
01070 else if(cLPDG > pLPDG) nLPDG=pLPDG*1000+cLPDG*100-1;
01071 else nLPDG=pLPDG*1000+cLPDG*100-3;
01072 }
01073 else if(cRPDG > 0 && pLPDG > 0)
01074 {
01075 order=-1;
01076 if (cRPDG > pLPDG) nLPDG=cRPDG*1000+pLPDG*100+1;
01077 else if(cRPDG < pLPDG) nLPDG=pLPDG*1000+cRPDG*100+1;
01078 else nLPDG=pLPDG*1000+cRPDG*100+3;
01079 if (cLPDG < pRPDG) nRPDG=cLPDG*1000+pRPDG*100-1;
01080 else if(cLPDG > pRPDG) nRPDG=pRPDG*1000+cLPDG*100-1;
01081 else nRPDG=pRPDG*1000+cLPDG*100-3;
01082 }
01083 else if(cRPDG < 0 && pLPDG < 0)
01084 {
01085 order=-1;
01086 if (cLPDG > pRPDG) nRPDG=cLPDG*1000+pRPDG*100+1;
01087 else if(cLPDG < pRPDG) nRPDG=pRPDG*1000+cLPDG*100+1;
01088 else nRPDG=pRPDG*1000+cLPDG*100+3;
01089 if (cRPDG < pLPDG) nLPDG=cRPDG*1000+pLPDG*100-1;
01090 else if(cRPDG > pLPDG) nLPDG=pLPDG*1000+cRPDG*100-1;
01091 else nLPDG=pLPDG*1000+cRPDG*100-3;
01092 }
01093 break;
01094 case 2:
01095 switch (af)
01096 {
01097 case 1:
01098 if(cLPDG < 0)
01099 {
01100 order= 1;
01101 if(-cLPDG==pRPDG)
01102 {
01103 nLPDG=pLPDG;
01104 nRPDG=cRPDG;
01105 }
01106 else
01107 {
01108 if (cRPDG > pRPDG) nRPDG=cRPDG*1000+pRPDG*100+1;
01109 else if(cRPDG < pRPDG) nRPDG=pRPDG*1000+cRPDG*100+1;
01110 else nRPDG=pRPDG*1000+cRPDG*100+3;
01111 if (-cLPDG == pL1) nLPDG=pL2;
01112 else nLPDG=pL1;
01113 }
01114 }
01115 else
01116 {
01117 order=-1;
01118 if(-cRPDG==pRPDG)
01119 {
01120 nLPDG=pLPDG;
01121 nRPDG=cLPDG;
01122 }
01123 else
01124 {
01125 if (cLPDG > pRPDG) nRPDG=cLPDG*1000+pRPDG*100+1;
01126 else if(cLPDG < pRPDG) nRPDG=pRPDG*1000+cLPDG*100+1;
01127 else nRPDG=pRPDG*1000+cLPDG*100+3;
01128 if (-cRPDG == pL1) nLPDG=pL2;
01129 else nLPDG=pL1;
01130 }
01131 }
01132 break;
01133 case 2:
01134 if(cLPDG < 0)
01135 {
01136 order=-1;
01137 if(-cLPDG==pLPDG)
01138 {
01139 nLPDG=cRPDG;
01140 nRPDG=pRPDG;
01141 }
01142 else
01143 {
01144 if (cRPDG > pLPDG) nLPDG=cRPDG*1000+pLPDG*100+1;
01145 else if(cRPDG < pLPDG) nLPDG=pLPDG*1000+cRPDG*100+1;
01146 else nLPDG=pLPDG*1000+cRPDG*100+3;
01147 if (-cLPDG == pR1) nRPDG=pR2;
01148 else nRPDG=pR1;
01149 }
01150 }
01151 else
01152 {
01153 order= 1;
01154 if(-cRPDG==pLPDG)
01155 {
01156 nLPDG=cLPDG;
01157 nRPDG=pRPDG;
01158 }
01159 else
01160 {
01161 if (cLPDG > pLPDG) nLPDG=cLPDG*1000+pLPDG*100+1;
01162 else if(cLPDG < pLPDG) nLPDG=pLPDG*1000+cLPDG*100+1;
01163 else nLPDG=pLPDG*1000+cLPDG*100+3;
01164 if (-cRPDG == pR1) nRPDG=pR2;
01165 else nRPDG=pR1;
01166 }
01167 }
01168 break;
01169 case 3:
01170 if(cLPDG > 0)
01171 {
01172 order= 1;
01173 if(cLPDG==-pRPDG)
01174 {
01175 nLPDG=pLPDG;
01176 nRPDG=cRPDG;
01177 }
01178 else
01179 {
01180 if (cRPDG < pRPDG) nRPDG=cRPDG*1000+pRPDG*100-1;
01181 else if(cRPDG > pRPDG) nRPDG=pRPDG*1000+cRPDG*100-1;
01182 else nRPDG=pRPDG*1000+cRPDG*100-3;
01183 if ( cLPDG == pL1) nLPDG=-pL2;
01184 else nLPDG=-pL1;
01185 }
01186 }
01187 else
01188 {
01189 order=-1;
01190 if(cRPDG==-pRPDG)
01191 {
01192 nLPDG=pLPDG;
01193 nRPDG=cLPDG;
01194 }
01195 else
01196 {
01197 if (cLPDG < pRPDG) nRPDG=cLPDG*1000+pRPDG*100-1;
01198 else if(cLPDG > pRPDG) nRPDG=pRPDG*1000+cLPDG*100-1;
01199 else nRPDG=pRPDG*1000+cLPDG*100-3;
01200 if ( cRPDG == pL1) nLPDG=-pL2;
01201 else nLPDG=-pL1;
01202 }
01203 }
01204 break;
01205 case 4:
01206 if(cLPDG > 0)
01207 {
01208 order=-1;
01209 if(cLPDG==-pLPDG)
01210 {
01211 nLPDG=cRPDG;
01212 nRPDG=pRPDG;
01213 }
01214 else
01215 {
01216 if (cRPDG < pLPDG) nLPDG=cRPDG*1000+pLPDG*100-1;
01217 else if(cRPDG > pLPDG) nLPDG=pLPDG*1000+cRPDG*100-1;
01218 else nLPDG=pLPDG*1000+cRPDG*100-3;
01219 if ( cLPDG == pR1) nRPDG=-pR2;
01220 else nRPDG=-pR1;
01221 }
01222 }
01223 else
01224 {
01225 order= 1;
01226 if(cRPDG==-pLPDG)
01227 {
01228 nLPDG=cLPDG;
01229 nRPDG=pRPDG;
01230 }
01231 else
01232 {
01233 if (cLPDG < pLPDG) nLPDG=cLPDG*1000+pLPDG*100-1;
01234 else if(cLPDG > pLPDG) nLPDG=pLPDG*1000+cLPDG*100-1;
01235 else nLPDG=pLPDG*1000+cLPDG*100-3;
01236 if ( cRPDG == pR1) nRPDG=-pR2;
01237 else nRPDG=-pR1;
01238 }
01239 }
01240 break;
01241 }
01242 break;
01243 case 3:
01244 switch (af)
01245 {
01246 case 1:
01247 if(pLPDG < 0)
01248 {
01249 order= 1;
01250 if(-pLPDG==cRPDG)
01251 {
01252 nLPDG=cLPDG;
01253 nRPDG=pRPDG;
01254 }
01255 else
01256 {
01257 if (pRPDG > cRPDG) nRPDG=pRPDG*1000+cRPDG*100+1;
01258 else if(pRPDG < cRPDG) nRPDG=cRPDG*1000+pRPDG*100+1;
01259 else nRPDG=cRPDG*1000+pRPDG*100+3;
01260 if (-pLPDG == L1) nLPDG=L2;
01261 else nLPDG=L1;
01262 }
01263 }
01264 else
01265 {
01266 order=-1;
01267 if(-pRPDG==cRPDG)
01268 {
01269 nLPDG=cLPDG;
01270 nRPDG=pLPDG;
01271 }
01272 else
01273 {
01274 if (pLPDG > cRPDG) nLPDG=pLPDG*1000+cRPDG*100+1;
01275 else if(pLPDG < cRPDG) nLPDG=cRPDG*1000+pLPDG*100+1;
01276 else nLPDG=cRPDG*1000+pLPDG*100+3;
01277 if (-pRPDG == L1) nRPDG=L2;
01278 else nRPDG=L1;
01279 }
01280 }
01281 break;
01282 case 2:
01283 if(pLPDG < 0)
01284 {
01285 order=-1;
01286 if(-pLPDG==cLPDG)
01287 {
01288 nLPDG=pRPDG;
01289 nRPDG=cRPDG;
01290 }
01291 else
01292 {
01293 if (pRPDG > cLPDG) nRPDG=pRPDG*1000+cLPDG*100+1;
01294 else if(pRPDG < cLPDG) nRPDG=cLPDG*1000+pRPDG*100+1;
01295 else nRPDG=cLPDG*1000+pRPDG*100+3;
01296 if (-pLPDG == R1) nLPDG=R2;
01297 else nLPDG=R1;
01298 }
01299 }
01300 else
01301 {
01302 order= 1;
01303 if(-pRPDG==cLPDG)
01304 {
01305 nLPDG=pLPDG;
01306 nRPDG=cRPDG;
01307 }
01308 else
01309 {
01310 if (pLPDG > cLPDG) nLPDG=pLPDG*1000+cLPDG*100+1;
01311 else if(pLPDG < cLPDG) nLPDG=cLPDG*1000+pLPDG*100+1;
01312 else nLPDG=cLPDG*1000+pLPDG*100+3;
01313 if (-pRPDG == R1) nRPDG=R2;
01314 else nRPDG=R1;
01315 }
01316 }
01317 break;
01318 case 3:
01319 if(pLPDG > 0)
01320 {
01321 order= 1;
01322 if(pLPDG==-cRPDG)
01323 {
01324 nLPDG=cLPDG;
01325 nRPDG=pRPDG;
01326 }
01327 else
01328 {
01329 if (pRPDG < cRPDG) nRPDG=pRPDG*1000+cRPDG*100-1;
01330 else if(pRPDG > cRPDG) nRPDG=cRPDG*1000+pRPDG*100-1;
01331 else nRPDG=cRPDG*1000+pRPDG*100-3;
01332 if ( pLPDG == L1) nLPDG=-L2;
01333 else nLPDG=-L1;
01334 }
01335 }
01336 else
01337 {
01338 order=-1;
01339 if(pRPDG==-cRPDG)
01340 {
01341 nLPDG=cLPDG;
01342 nRPDG=pLPDG;
01343 }
01344 else
01345 {
01346 if (pLPDG < cRPDG) nLPDG=pLPDG*1000+cRPDG*100-1;
01347 else if(pLPDG > cRPDG) nLPDG=cRPDG*1000+pLPDG*100-1;
01348 else nLPDG=cRPDG*1000+pLPDG*100-3;
01349 if ( pRPDG == L1) nRPDG=-L2;
01350 else nRPDG=-L1;
01351 }
01352 }
01353 break;
01354 case 4:
01355 if(pLPDG > 0)
01356 {
01357 order=-1;
01358 if(pLPDG==-cLPDG)
01359 {
01360 nLPDG=pRPDG;
01361 nRPDG=cRPDG;
01362 }
01363 else
01364 {
01365 if (pRPDG < cLPDG) nRPDG=pRPDG*1000+cLPDG*100-1;
01366 else if(pRPDG > cLPDG) nRPDG=cLPDG*1000+pRPDG*100-1;
01367 else nRPDG=cLPDG*1000+pRPDG*100-3;
01368 if ( pLPDG == R1) nLPDG=-R2;
01369 else nLPDG=-R1;
01370 }
01371 }
01372 else
01373 {
01374 order= 1;
01375 if(pRPDG==-cLPDG)
01376 {
01377 nLPDG=pLPDG;
01378 nRPDG=cRPDG;
01379 }
01380 else
01381 {
01382 if (pLPDG < cLPDG) nLPDG=pLPDG*1000+cLPDG*100-1;
01383 else if(pLPDG > cLPDG) nLPDG=cLPDG*1000+pLPDG*100-1;
01384 else nLPDG=cLPDG*1000+pLPDG*100-3;
01385 if ( pRPDG == R1) nRPDG=-R2;
01386 else nRPDG=-R1;
01387 }
01388 }
01389 break;
01390 }
01391 break;
01392 case 4:
01393 switch (af)
01394 {
01395 case 1:
01396 order= 1;
01397 if(-cLPDG == pL1) nLPDG= pL2;
01398 else nLPDG= pL1;
01399 if( cRPDG == pR1) nRPDG=-pR2;
01400 else nRPDG=-pR1;
01401 break;
01402 case 2:
01403 order=-1;
01404 if(-cRPDG == pL1) nLPDG= pL2;
01405 else nLPDG= pL1;
01406 if( cLPDG == pR1) nRPDG=-pR2;
01407 else nRPDG=-pR1;
01408 break;
01409 case 3:
01410 order= 1;
01411 if( cLPDG == pL1) nLPDG=-pL2;
01412 else nLPDG=-pL1;
01413 if(-cRPDG == pR1) nRPDG= pR2;
01414 else nRPDG= pR1;
01415 break;
01416 case 4:
01417 order=-1;
01418 if( cRPDG == pL1) nLPDG=-pL2;
01419 else nLPDG=-pL1;
01420 if(-cLPDG == pR1) nRPDG= pR2;
01421 else nRPDG= pR1;
01422 break;
01423 }
01424 break;
01425 case 5:
01426 switch (af)
01427 {
01428 case 1:
01429 order= 1;
01430 if(-pLPDG == L1) nLPDG= L2;
01431 else nLPDG= L1;
01432 if( pRPDG == R1) nRPDG=-R2;
01433 else nRPDG=-R1;
01434 break;
01435 case 2:
01436 order=-1;
01437 if(-pRPDG == L1) nRPDG= L2;
01438 else nRPDG= L1;
01439 if( pLPDG == R1) nLPDG=-R2;
01440 else nLPDG=-R1;
01441 break;
01442 case 3:
01443 order= 1;
01444 if( pLPDG == L1) nLPDG=-L2;
01445 else nLPDG=-L1;
01446 if(-pRPDG == R1) nRPDG= R2;
01447 else nRPDG= R1;
01448 break;
01449 case 4:
01450 order=-1;
01451 if( pRPDG == L1) nRPDG=-L2;
01452 else nRPDG=-L1;
01453 if(-pLPDG == R1) nLPDG= R2;
01454 else nLPDG= R1;
01455 break;
01456 }
01457 break;
01458 case 6:
01459 switch (af)
01460 {
01461 case 1:
01462 order=-1;
01463 if(-cRPDG == pL1) nLPDG= pL2;
01464 else nLPDG= pL1;
01465 if( pRPDG == L1) nRPDG= -L2;
01466 else nRPDG= -L1;
01467 break;
01468 case 2:
01469 order= 1;
01470 if(-cLPDG == pL1) nLPDG= pL2;
01471 else nLPDG= pL1;
01472 if( pRPDG == R1) nRPDG= -R2;
01473 else nRPDG= -R1;
01474 break;
01475 case 3:
01476 order= 1;
01477 if(-cRPDG == pR1) nRPDG= pR2;
01478 else nRPDG= pR1;
01479 if( pLPDG == L1) nLPDG= -L2;
01480 else nLPDG= -L1;
01481 break;
01482 case 4:
01483 order=-1;
01484 if(-cLPDG == pR1) nRPDG= pR2;
01485 else nRPDG= pR1;
01486 if( pLPDG == R1) nLPDG= -R2;
01487 else nLPDG= -R1;
01488 break;
01489 case 5:
01490 order=-1;
01491 if(-pRPDG == L1) nRPDG= L2;
01492 else nRPDG= L1;
01493 if( cRPDG == pL1) nLPDG=-pL2;
01494 else nLPDG=-pL1;
01495 break;
01496 case 6:
01497 order= 1;
01498 if(-pLPDG == L1) nLPDG= L2;
01499 else nLPDG= L1;
01500 if( cRPDG == pR1) nRPDG=-pR2;
01501 else nRPDG=-pR1;
01502 break;
01503 case 7:
01504 order= 1;
01505 if(-pRPDG == R1) nRPDG= R2;
01506 else nRPDG= R1;
01507 if( cLPDG == pL1) nLPDG=-pL2;
01508 else nLPDG=-pL1;
01509 break;
01510 case 8:
01511 order=-1;
01512 if(-pLPDG == R1) nLPDG= R2;
01513 else nLPDG= R1;
01514 if( cLPDG == pR1) nRPDG=-pR2;
01515 else nRPDG=-pR1;
01516 break;
01517 }
01518 break;
01519 }
01520 if(!order) G4cerr<<"-Warning-G4QFrag::Constr: t="<<tf<<", a="<<af<<", cL="<<cLPDG
01521 <<", cR="<<cRPDG<<", pL="<<pLPDG<<", pR="<<pRPDG<<G4endl;
01522 else
01523 {
01524
01525 G4int LT=1;
01526 if(std::abs(nLPDG) > 7) ++LT;
01527 G4int RT=1;
01528 if(std::abs(nRPDG) > 7) ++RT;
01529 G4double minM=0.;
01530 G4bool sing=true;
01531 if(cLT==2 && cRT==2)
01532 {
01533 aLPDG=0;
01534 aRPDG=0;
01535 if(cLPDG>0)
01536 {
01537 aLPDG=nLPDG/100;
01538 aRPDG=(-nRPDG)/100;
01539 }
01540 else
01541 {
01542 aRPDG=nRPDG/100;
01543 aLPDG=(-nLPDG)/100;
01544 }
01545 G4int nL1=aLPDG/10;
01546 G4int nL2=aLPDG%10;
01547 G4int nR1=aRPDG/10;
01548 G4int nR2=aRPDG%10;
01549 if(nL1!=nR1 && nL1!=nR2 && nL2!=nR1 && nL2!=nR2)
01550 {
01551 #ifdef debug
01552 G4cout<<"G4QFragmentation::Const:aLPDG="<<aLPDG<<", aRPDG="<<aRPDG<<G4endl;
01553 #endif
01554 sing=false;
01555 G4QPDGCode tmp;
01556 std::pair<G4int,G4int> pB=tmp.MakeTwoBaryons(nL1, nL2, nR1, nR2);
01557 minM=G4QPDGCode(pB.first).GetMass()+G4QPDGCode(pB.second).GetMass();
01558 }
01559 }
01560 if(sing)
01561 {
01562 std::pair<G4int,G4int> newPair = std::make_pair(nLPDG,nRPDG);
01563 G4QContent newStQC(newPair);
01564 #ifdef debug
01565 G4cout<<"G4QFr::Con: LPDG="<<nLPDG<<",RPDG="<<nRPDG<<",QC="<<newStQC<<G4endl;
01566 #endif
01567 G4int minPDG=newStQC.GetSPDGCode();
01568 minM=G4QPDGCode(minPDG).GetMass() + mPi0;
01569 }
01570
01571 G4bool win=false;
01572 G4double pSM=0.;
01573 if(pSM2>0.) pSM=std::sqrt(pSM2);
01574 if(minC && pSM2 > maxiM2)
01575 {
01576 maxiM2=pSM2;
01577 win=true;
01578 }
01579 else if(!minC || pSM > minM)
01580 {
01581 pExcess=pSM-minM;
01582 if(minC || pExcess > excess)
01583 {
01584 minC=false;
01585 excess=pExcess;
01586 win=true;
01587 }
01588 }
01589 if(win)
01590 {
01591 sst=pst;
01592 sLPDG=nLPDG;
01593 sRPDG=nRPDG;
01594 sOrd=order;
01595 }
01596 }
01597 }
01598
01599 }
01600 if(sOrd)
01601 {
01602 G4LorentzVector cL4M=cLeft->Get4Momentum();
01603 G4LorentzVector cR4M=cRight->Get4Momentum();
01604 G4QParton* pLeft=(*sst)->GetLeftParton();
01605 G4QParton* pRight=(*sst)->GetRightParton();
01606 G4LorentzVector pL4M=pLeft->Get4Momentum();
01607 G4LorentzVector pR4M=pRight->Get4Momentum();
01608 #ifdef debug
01609 G4cout<<"G4QFragmentation::Const:cS4M="<<cS4M<<" fused/w pS4M="<<pL4M+pR4M<<G4endl;
01610 #endif
01611 if(sOrd>0)
01612 {
01613 pL4M+=cL4M;
01614 pR4M+=cR4M;
01615 }
01616 else
01617 {
01618 pL4M+=cR4M;
01619 pR4M+=cL4M;
01620 }
01621 pLeft->SetPDGCode(sLPDG);
01622 pLeft->Set4Momentum(pL4M);
01623 pRight->SetPDGCode(sRPDG);
01624 pRight->Set4Momentum(pR4M);
01625 delete (*ist);
01626 strings.erase(ist);
01627 #ifdef debug
01628 G4LorentzVector ss4M=pL4M+pR4M;
01629 G4cout<<"G4QFragmentation::Construct:Created,4M="<<ss4M<<",m2="<<ss4M.m2()<<G4endl;
01630 #endif
01631 if( ist != strings.begin() )
01632 {
01633 ist--;
01634 con=false;
01635 #ifdef debug
01636 G4cout<<"G4QFragmentation::Construct: *** IST Decremented ***"<<G4endl;
01637 #endif
01638 }
01639 else
01640 {
01641 con=true;
01642 #ifdef debug
01643 G4cout<<"G4QFragmentation::Construct: *** IST Begin ***"<<G4endl;
01644 #endif
01645 break;
01646 }
01647 }
01648 else
01649 {
01650 #ifdef debug
01651 G4cout<<"-Warning-G4QFragm::Const:S4M="<<cS4M<<",M2="<<cSM2<<" Leave ASIS"<<G4endl;
01652 #endif
01653 ++problem;
01654 con=false;
01655 }
01656 }
01657 else con=false;
01658 }
01659 #ifdef debug
01660 G4cout<<"G4QFragmentation::Construct: *** IST While *** , con="<<con<<G4endl;
01661 #endif
01662 }
01663 #ifdef edebug
01664
01665 G4LorentzVector t4M=theNucleus.Get4Momentum();
01666 G4int rC=totChg-theNucleus.GetZ();
01667 G4int rB=totBaN-theNucleus.GetA();
01668 G4int nStr=strings.size();
01669 G4cout<<"-EMCLS-G4QFr::Const: AfterSUPPRESION #ofS="<<nStr<<",tNuc4M(E=M)="<<sum<<G4endl;
01670 for(G4int i=0; i<nStr; i++)
01671 {
01672 G4LorentzVector strI4M=strings[i]->Get4Momentum();
01673 t4M+=strI4M;
01674 G4int sChg=strings[i]->GetCharge();
01675 rC-=sChg;
01676 G4int sBaN=strings[i]->GetBaryonNumber();
01677 rB-=sBaN;
01678 G4cout<<"-EMCLS-G4QFragm::Construct: St#"<<i<<", 4M="<<strI4M<<", M="<<strI4M.m()
01679 <<", C="<<sChg<<", B="<<sBaN<<G4endl;
01680 }
01681 G4cout<<"-EMCLS-G4QFragm::Construct:r4M="<<t4M-totLS4M<<",rC="<<rC<<",rB="<<rB<<G4endl;
01682 #endif
01683
01684
01685
01686 #ifdef debug
01687 G4cout<<"G4QFragmentation::Construct: problem="<<problem<<G4endl;
01688 #endif
01689 if(problem)
01690 {
01691 G4int nOfStr=strings.size();
01692 #ifdef debug
01693 G4cout<<"G4QFragmentation::Construct:SecurityDiQaDiQReduction,#OfStr="<<nOfStr<<G4endl;
01694 #endif
01695 for (G4int astring=0; astring < nOfStr; astring++)
01696 {
01697 G4QString* curString=strings[astring];
01698 G4QParton* cLeft=curString->GetLeftParton();
01699 G4QParton* cRight=curString->GetRightParton();
01700 G4int LT=cLeft->GetType();
01701 G4int RT=cRight->GetType();
01702 G4int sPDG=cLeft->GetPDGCode();
01703 G4int nPDG=cRight->GetPDGCode();
01704 if(LT==2 && RT==2)
01705 {
01706 #ifdef debug
01707 G4cout<<"G4QFragmentation::Constr:TrySelfReduString,L="<<sPDG<<",R="<<nPDG<<G4endl;
01708 #endif
01709 if( cLeft->ReduceDiQADiQ(cLeft, cRight) )
01710 {
01711 sPDG=cLeft->GetPDGCode();
01712 nPDG=cRight->GetPDGCode();
01713 #ifdef debug
01714 G4cout<<"+G4QFragm::Const:#"<<astring<<" Reduced, L="<<sPDG<<",R="<<nPDG<<G4endl;
01715 #endif
01716 }
01717 #ifdef debug
01718 else G4cout<<"--*--G4QFragm::Const:#"<<astring<<" DQ-aDQ reduction Failed"<<G4endl;
01719 #endif
01720 }
01721 else if(sPDG==3 && nPDG==-3)
01722 {
01723 sPDG= 1;
01724 nPDG=-1;
01725 cLeft->SetPDGCode(sPDG);
01726 cRight->SetPDGCode(nPDG);
01727 }
01728 else if(sPDG==-3 && nPDG==3)
01729 {
01730 sPDG=-1;
01731 nPDG= 1;
01732 cLeft->SetPDGCode(sPDG);
01733 cRight->SetPDGCode(nPDG);
01734 }
01735 }
01736 SwapPartons();
01737 }
01738 #ifdef edebug
01739 G4LorentzVector u4M=theNucleus.Get4Momentum();
01740 G4int rCh=totChg-theNucleus.GetZ();
01741 G4int rBa=totBaN-theNucleus.GetA();
01742 G4int nStri=strings.size();
01743 G4cout<<"-EMCLS-G4QFr::Const: FinalConstruct, #ofSt="<<nStri<<",tN4M(E=M)="<<t4M<<G4endl;
01744 for(G4int i=0; i<nStri; i++)
01745 {
01746 G4LorentzVector strI4M=strings[i]->Get4Momentum();
01747 u4M+=strI4M;
01748 G4int sChg=strings[i]->GetCharge();
01749 rCh-=sChg;
01750 G4int sBaN=strings[i]->GetBaryonNumber();
01751 rBa-=sBaN;
01752 G4cout<<"-EMCLS-G4QFragm::Construct: St#"<<i<<", 4M="<<strI4M<<", M="<<strI4M.m()
01753 <<", C="<<sChg<<", B="<<sBaN<<G4endl;
01754 }
01755 G4cout<<"-EMCLS-G4QFragm::Construct:r4M="<<u4M-totLS4M<<",rC="<<rCh<<",rB="<<rBa<<G4endl;
01756 #endif
01757 }
01758
01759 G4QFragmentation::~G4QFragmentation()
01760 {
01761 std::for_each(strings.begin(), strings.end(), DeleteQString() );
01762 }
01763
01764 G4QHadronVector* G4QFragmentation::Fragment()
01765 {
01766 static const G4QPDGCode nQPDG(2112);
01767 static const G4double mProt = G4QPDGCode(2212).GetMass();
01768 static const G4double mNeut = G4QPDGCode(2112).GetMass();
01769 static const G4double mPiCh = G4QPDGCode(211).GetMass();
01770 static const G4double mPiZr = G4QPDGCode(111).GetMass();
01771
01772 static const G4LorentzVector nul4M(0.,0.,0.,0.);
01773
01774 #ifdef debug
01775 G4cout<<"*******>G4QFragmentation::Fragment: ***Called***, Res="<<theResult<<G4endl;
01776 #endif
01777 G4int striNum=strings.size();
01778 G4int hadrNum=theResult->size();
01779 #ifdef edebug
01780 G4int nQm=theQuasmons.size();
01781 G4LorentzVector totLS4M=theNucleus.Get4Momentum();
01782 G4int totChg=theNucleus.GetZ();
01783 G4int totBaN=theNucleus.GetA();
01784 G4cout<<"-EMCLS-G4QF::Fragment: CHECKRecovery, #ofS="<<striNum<<", #Nuc4M(E=M)="<<totLS4M
01785 <<",#Q="<<nQm<<",#H="<<hadrNum<<G4endl;
01786 for(G4int i=0; i < striNum; i++)
01787 {
01788 G4LorentzVector strI4M=strings[i]->Get4Momentum();
01789 totLS4M+=strI4M;
01790 G4int sChg=strings[i]->GetCharge();
01791 totChg+=sChg;
01792 G4int sBaN=strings[i]->GetBaryonNumber();
01793 totBaN+=sBaN;
01794 G4cout<<"-EMCLS-G4QFragm::Fragment: String#"<<i<<", 4M="<<strI4M<<", M="<<strI4M.m()
01795 <<", C="<<sChg<<", B="<<sBaN<<G4endl;
01796 }
01797 for(G4int i=0; i < nQm; i++)
01798 {
01799 G4LorentzVector hI4M=theQuasmons[i]->Get4Momentum();
01800 totLS4M+=hI4M;
01801 G4int hChg=theQuasmons[i]->GetCharge();
01802 totChg+=hChg;
01803 G4int hBaN=theQuasmons[i]->GetBaryonNumber();
01804 totBaN+=hBaN;
01805 G4cout<<"-EMCLS-G4QFragmentation::Fragment: Quasmon#"<<i<<", 4M="<<hI4M<<", C="<<hChg
01806 <<", B="<<hBaN<<G4endl;
01807 }
01808 for(G4int i=0; i < hadrNum; i++)
01809 {
01810 G4LorentzVector hI4M=(*theResult)[i]->Get4Momentum();
01811 totLS4M+=hI4M;
01812 G4int hChg=(*theResult)[i]->GetCharge();
01813 totChg+=hChg;
01814 G4int hBaN=(*theResult)[i]->GetBaryonNumber();
01815 totBaN+=hBaN;
01816 G4cout<<"-EMCLS-G4QFr::Fragment:H#"<<i<<",4M="<<hI4M<<",C="<<hChg<<",B="<<hBaN<<G4endl;
01817 }
01818 #endif
01819 #ifdef debug
01820 G4cout<<"***>G4QFragmentation::Fragment: #OfStr="<<striNum<<", #OfRes="<<hadrNum<<G4endl;
01821 #endif
01822 if(!striNum && hadrNum)
01823 {
01824 #ifdef debug
01825 G4cout<<"***>G4QFragmentation::Fragment:**Quasi-Elastic**,#OfResult="<<hadrNum<<G4endl;
01826 #endif
01827 return theResult;
01828 }
01829 else if(striNum) Breeder();
01830 else
01831 {
01832 if(hadrNum)
01833 {
01834 for(G4int ih=0; ih<hadrNum; ih++) delete (*theResult)[ih];
01835 theResult->clear();
01836 }
01837 G4LorentzVector r4M=theNucleus.Get4Momentum();
01838 G4int rPDG=theNucleus.GetPDG();
01839 G4QHadron* resNuc = new G4QHadron(rPDG,r4M);
01840 theResult->push_back(resNuc);
01841 }
01842 G4int nQuas=theQuasmons.size();
01843 G4int theRS=theResult->size();
01844 #ifdef debug
01845 G4cout<<"***>G4QFragmentation::Fragment:beforeEnv,#OfQ="<<nQuas<<",#OfR="<<theRS<<G4endl;
01846 #endif
01847 if(nQuas && theRS)
01848 {
01849 G4QHadron* resNuc = (*theResult)[theRS-1];
01850 G4LorentzVector resNuc4M = resNuc->Get4Momentum();
01851 G4int resNucPDG= resNuc->GetPDGCode();
01852 if(resNucPDG==90000000 || resNuc4M.m2()<800000.)
01853 {
01854 resNuc4M=G4LorentzVector(0.,0.,0.,0.);
01855 if(resNucPDG == 90000000) resNuc->Set4Momentum(resNuc4M);
01856 }
01857 #ifdef edebug
01858 G4int rnChg=resNuc->GetCharge();
01859 G4int rnBaN=resNuc->GetBaryonNumber();
01860 #endif
01861 G4QNucleus theEnv(resNucPDG);
01862 delete resNuc;
01863 theResult->pop_back();
01864 --theRS;
01865 #ifdef debug
01866 G4cout<<"G4QFragmentation::Fragment:#OfRemainingHadron="<<theRS<<",A="<<theEnv<<G4endl;
01867 #endif
01868
01869 for(G4int j=theRS-1; j>-2; --j)
01870 {
01871 G4LorentzVector qsum4M=resNuc4M;
01872 G4QContent qsumQC=theEnv.GetQCZNS();
01873 #ifdef debug
01874 G4cout<<"G4QFragmentation::Fragm:rN4M"<<qsum4M<<qsum4M.m()<<",rNQC="<<qsumQC<<G4endl;
01875 #endif
01876 G4Quasmon* firstQ=0;
01877 G4LorentzVector first4M;
01878 G4QContent firstQC;
01879 for(G4int i=0; i<nQuas; ++i)
01880 {
01881 G4Quasmon* curQuasm=theQuasmons[i];
01882 G4LorentzVector cur4M=curQuasm->Get4Momentum();
01883 G4QContent curQC=curQuasm->GetQC();
01884 qsum4M+=cur4M;
01885 qsumQC+=curQC;
01886 #ifdef debug
01887 G4cout<<"G4QFr::Fr:Q#"<<i<<",Q4M="<<cur4M<<",QQC="<<curQC<<",sQC="<<qsumQC<<G4endl;
01888 #endif
01889 if(!i)
01890 {
01891 firstQ =curQuasm;
01892 first4M=cur4M;
01893 firstQC=curQC;
01894 }
01895 }
01896 G4int miPDG=qsumQC.GetSPDGCode();
01897 G4double gsM=0.;
01898 #ifdef debug
01899 G4cout<<"G4QFragmentation::Fragment: QC="<<qsumQC<<",PDG="<<miPDG<<G4endl;
01900 #endif
01901 if(miPDG == 10)
01902 {
01903 G4QChipolino QCh(qsumQC);
01904 gsM=QCh.GetQPDG1().GetMass()+QCh.GetQPDG2().GetMass();
01905
01906
01907 }
01908
01909
01910
01911
01912
01913
01914 else if(miPDG < 80000000 && std::abs(miPDG)%10 > 2)
01915 gsM=theWorld->GetQParticle(G4QPDGCode(miPDG))->MinMassOfFragm();
01916 else gsM=G4QPDGCode(miPDG).GetMass();
01917 G4double reM=qsum4M.m();
01918 #ifdef debug
01919 G4cout<<"G4QFragmentation::Fragment: PDG="<<miPDG<<",rM="<<reM<<",GSM="<<gsM<<G4endl;
01920 #endif
01921 if(reM > gsM) break;
01922 if(j > -1)
01923 {
01924 G4QHadron* cH = (*theResult)[j];
01925 G4LorentzVector h4M = cH->Get4Momentum();
01926 G4QContent hQC = cH->GetQC();
01927 firstQ->Set4Momentum(first4M+h4M);
01928 firstQ->SetQC(firstQC+hQC);
01929 delete cH;
01930 theResult->pop_back();
01931 #ifdef debug
01932 G4cout<<"G4QFragm::Fragm: H#"<<j<<",hQC="<<hQC<<",hPDG="<<cH->GetPDGCode()<<G4endl;
01933 #endif
01934 }
01935 else
01936 {
01937 G4cerr<<"*G4QFr::Fr:PDG="<<miPDG<<",M="<<reM<<",GSM="<<gsM<<",QC="<<qsumQC<<G4endl;
01938 G4Exception("G4QFragmentation::Fragment:","27",FatalException,"Can't recover GSM");
01939 }
01940 }
01941 G4double nucE=resNuc4M.e();
01942 if(nucE < 1.E-12) nucE=0.;
01943 else if(resNucPDG==22 && nQuas==1)
01944 {
01945 G4Quasmon* aQuasm=theQuasmons[0];
01946 aQuasm->Set4Momentum(aQuasm->Get4Momentum()+resNuc4M);
01947 nucE=0.;
01948 }
01949 G4ThreeVector nucVel(0.,0.,0.);
01950 G4QHadronVector* output=0;
01951 G4QEnvironment* pan= new G4QEnvironment(theEnv);
01952 #ifdef debug
01953 G4cout<<"G4QFragm::Fragm: nucE="<<nucE<<",nQ="<<nQuas<<G4endl;
01954 #endif
01955 if(nucE) nucVel=resNuc4M.vect()/nucE;
01956 #ifdef edebug
01957 G4LorentzVector sq4M=resNuc4M-totLS4M;
01958 G4int sqCg=rnChg-totChg;
01959 G4int sqBN=rnBaN-totBaN;
01960 #endif
01961 for(G4int i=0; i<nQuas; ++i)
01962 {
01963 G4Quasmon* curQuasm=theQuasmons[i];
01964 #ifdef debug
01965 if(nucE) G4cout<<"G4QFr::Fr:V="<<nucVel<<",Q="<<curQuasm->Get4Momentum()<<",R="
01966 <<resNuc4M<<resNucPDG<<G4endl;
01967 #endif
01968 if(nucE) curQuasm->Boost(-nucVel);
01969 pan->AddQuasmon(curQuasm);
01970 #ifdef edebug
01971 G4LorentzVector cQ4M=curQuasm->Get4Momentum();
01972 G4cout<<"G4QFragmentation::Fragment: Quasmon# "<<i<<" added, 4M="<<cQ4M<<G4endl;
01973 sq4M+=cQ4M;
01974 sqCg+=curQuasm->GetCharge();
01975 sqBN+=curQuasm->GetBaryonNumber();
01976 #endif
01977 }
01978 #ifdef edebug
01979 G4cout<<"-EMCLS-G4QFrag::Fragm: r4M="<<sq4M<<", rC="<<sqCg<<", rB="<<sqBN<<G4endl;
01980 #endif
01981 try
01982 {
01983 #ifdef debug
01984 G4cout<<"G4QFrag::Fragm: *** Before Del Output ***"<<G4endl;
01985 #endif
01986 delete output;
01987 #ifdef debug
01988 G4cout<<"G4QFrag::Fragm: *** After Del Output ***"<<G4endl;
01989 #endif
01990 output = pan->Fragment();
01991 }
01992 catch (G4QException& error)
01993 {
01994 G4cerr<<"***G4QFragmentation::Fragment: G4QE Exception is catched"<<G4endl;
01995
01996 G4Exception("G4QFragmentation::Fragment()", "HAD_CHPS_0027",
01997 FatalException, "CHIPSCrash");
01998 }
01999 #ifdef debug
02000 G4cout<<"G4QFrag::Fragm: *** Before Del Pan ***"<<G4endl;
02001 #endif
02002 delete pan;
02003 #ifdef debug
02004 G4cout<<"G4QFrag::Fragm: *** After Del Pan ***"<<G4endl;
02005 #endif
02006 if(output)
02007 {
02008 G4int nOut=output->size();
02009 for(G4int j=0; j<nOut; j++)
02010 {
02011 G4QHadron* curHadr=(*output)[j];
02012 if(nucE) curHadr->Boost(nucVel);
02013 G4int hPDG=curHadr->GetPDGCode();
02014 G4LorentzVector h4M=curHadr->Get4Momentum();
02015 if((hPDG!=90000000 && hPDG!=22) || h4M!=nul4M) theResult->push_back(curHadr);
02016 else delete curHadr;
02017 }
02018 delete output;
02019 }
02020 }
02021 else if(!striNum) G4cout<<"-Warning-G4QFragmentation::Fragment:Nothing was done"<<G4endl;
02022 #ifdef debug
02023 G4cout<<"=--=>G4QFragmentation::Fragment: Final #OfResult="<<theResult->size()<<G4endl;
02024 #endif
02025 G4int nQ =theQuasmons.size();
02026 if(nQ) theQuasmons.clear();
02027 G4int nHd=theResult->size();
02028 #ifdef edebug
02029 G4LorentzVector f4M(0.,0.,0.,0.);
02030 G4int fCh=totChg;
02031 G4int fBN=totBaN;
02032 G4cout<<"-EMCLS-G4QFragmentation::Fragment: #ofHadr="<<nHd<<", #OfQuasm="<<nQ<<G4endl;
02033 for(G4int i=0; i<nHd; i++)
02034 {
02035 G4LorentzVector hI4M=(*theResult)[i]->Get4Momentum();
02036 f4M+=hI4M;
02037 G4int hChg=(*theResult)[i]->GetCharge();
02038 fCh-=hChg;
02039 G4int hBaN=(*theResult)[i]->GetBaryonNumber();
02040 fBN-=hBaN;
02041 G4cout<<"-EMCLS-G4QFragmentation::Fragment: Hadron#"<<i<<", 4M="<<hI4M<<", PDG="
02042 <<(*theResult)[i]->GetPDGCode()<<", C="<<hChg<<", B="<<hBaN<<G4endl;
02043 }
02044 G4cout<<"-EMCLS-G4QFrag::Fragm: r4M="<<f4M-totLS4M<<", rC="<<fCh<<", rB="<<fBN<<G4endl;
02045 #endif
02046
02047 G4QHadron* resNuc = (*theResult)[nHd-1];
02048 G4int rnBn = resNuc->GetBaryonNumber();
02049 G4int rnCg = resNuc->GetCharge();
02050 if(rnBn==1 && (rnCg==-2 || rnCg==3 || rnCg==-1 || rnCg==2))
02051 {
02052 G4LorentzVector tot4M=resNuc->Get4Momentum();
02053 G4int nPDG=2212;
02054 G4int mPDG=211;
02055 G4double nM=mProt;
02056 if(rnCg<0)
02057 {
02058 nPDG=2112;
02059 mPDG=-211;
02060 nM=mNeut;
02061 }
02062 G4LorentzVector m14M(0.,0.,0.,mPiCh);
02063 G4LorentzVector n4M(0.,0.,0.,nM);
02064 if(rnCg==-2 || rnCg==3)
02065 {
02066 G4LorentzVector m24M(0.,0.,0.,mPiCh);
02067 if(!G4QHadron(tot4M).DecayIn3(m14M,m24M,n4M))
02068 {
02069 G4cerr<<"***G4QFrag::Frag: tM="<<tot4M.m()<<" -> m1="<<mPiCh<<" + m2="<<mPiCh
02070 <<" + nM="<<nM<<" = "<<2*mPiCh+nM<<G4endl;
02071 G4Exception("G4QFragmentation::Breeder:","72",FatalException,"ImpossibleDecayIn3");
02072 }
02073 theResult->pop_back();
02074 delete resNuc;
02075 G4QHadron* m1H = new G4QHadron(mPDG,m14M);
02076 theResult->push_back(m1H);
02077 #ifdef debug
02078 G4cout<<"G4QFragment::Fragment:DecayIn3, M1="<<mPDG<<m14M<<G4endl;
02079 #endif
02080 G4QHadron* m2H = new G4QHadron(mPDG,m24M);
02081 theResult->push_back(m2H);
02082 #ifdef debug
02083 G4cout<<"G4QFragment::Fragment:DecayIn3, M2="<<mPDG<<m24M<<G4endl;
02084 #endif
02085 G4QHadron* nH = new G4QHadron(nPDG,n4M);
02086 theResult->push_back(nH);
02087 #ifdef debug
02088 G4cout<<"G4QFragment::Fragment:DecayIn3, Nucleon="<<nPDG<<n4M<<G4endl;
02089 #endif
02090 }
02091 else
02092 {
02093 if(!G4QHadron(tot4M).DecayIn2(m14M,n4M))
02094 {
02095 G4cerr<<"***G4QFrag::Frag: tM="<<tot4M.m()<<" -> m1="<<mPiCh
02096 <<" + nM="<<nM<<" = "<<mPiCh+nM<<G4endl;
02097 G4Exception("G4QFragmentation::Breeder:","72",FatalException,"ImpossibleDecayIn2");
02098 }
02099 theResult->pop_back();
02100 delete resNuc;
02101 G4QHadron* m1H = new G4QHadron(mPDG,m14M);
02102 theResult->push_back(m1H);
02103 #ifdef debug
02104 G4cout<<"G4QFragment::Fragment:DecayIn2, M1="<<mPDG<<m14M<<G4endl;
02105 #endif
02106 G4QHadron* nH = new G4QHadron(nPDG,n4M);
02107 theResult->push_back(nH);
02108 #ifdef debug
02109 G4cout<<"G4QFragment::Fragment:DecayIn2, Nucleon="<<nPDG<<n4M<<G4endl;
02110 #endif
02111 }
02112 }
02113 if(rnBn==2)
02114 {
02115 if(!rnCg)
02116 {
02117 G4LorentzVector tot4M=resNuc->Get4Momentum();
02118 G4LorentzVector n14M(0.,0.,0.,mNeut);
02119 G4LorentzVector n24M(0.,0.,0.,mNeut);
02120 if(!G4QHadron(tot4M).DecayIn2(n14M,n24M))
02121 {
02122 G4cerr<<"***G4QFrag::Frag: tM="<<tot4M.m()<<" -> n*2="<<2*mNeut<<G4endl;
02123 G4Exception("G4QFragmentation::Breeder:","72",FatalException,"ImpossibleDecay-2n");
02124 }
02125 theResult->pop_back();
02126 delete resNuc;
02127 G4QHadron* n1H = new G4QHadron(2112,n14M);
02128 theResult->push_back(n1H);
02129 #ifdef debug
02130 G4cout<<"G4QFragment::Fragment:DecayIn2, Neutron1="<<n14M<<G4endl;
02131 #endif
02132 G4QHadron* n2H = new G4QHadron(2112,n24M);
02133 theResult->push_back(n2H);
02134 #ifdef debug
02135 G4cout<<"G4QFragment::Fragment:DecayIn2, Neutron2="<<n24M<<G4endl;
02136 #endif
02137 }
02138 else if(rnCg==2)
02139 {
02140 G4LorentzVector tot4M=resNuc->Get4Momentum();
02141 G4LorentzVector n14M(0.,0.,0.,mProt);
02142 G4LorentzVector n24M(0.,0.,0.,mProt);
02143 if(!G4QHadron(tot4M).DecayIn2(n14M,n24M))
02144 {
02145 G4cerr<<"***G4QFrag::Frag: tM="<<tot4M.m()<<" -> n*2="<<2*mProt<<G4endl;
02146 G4Exception("G4QFragmentation::Breeder:","72",FatalException,"ImpossibleDecay-2p");
02147 }
02148 theResult->pop_back();
02149 delete resNuc;
02150 G4QHadron* n1H = new G4QHadron(2212,n14M);
02151 theResult->push_back(n1H);
02152 #ifdef debug
02153 G4cout<<"G4QFragment::Fragment:DecayIn2, Proton1="<<n14M<<G4endl;
02154 #endif
02155 G4QHadron* n2H = new G4QHadron(2212,n24M);
02156 theResult->push_back(n2H);
02157 #ifdef debug
02158 G4cout<<"G4QFragment::Fragment:DecayIn2, Proton2="<<n24M<<G4endl;
02159 #endif
02160 }
02161 }
02162
02163 nHd=theResult->size();
02164 G4int maxChg=0;
02165 #ifdef debug
02166 G4int maxBN=0;
02167 #endif
02168 G4QContent maxQC(0,0,0,0,0,0);
02169 for(G4int i=0; i<nHd; ++i)
02170 {
02171 G4int found=0;
02172 G4QHadron* cHadr = (*theResult)[i];
02173 G4int hPDG= cHadr->GetPDGCode();
02174 if(hPDG==90000000 || hPDG==22)
02175 {
02176 G4QHadron* curHadr=(*theResult)[i];
02177 G4LorentzVector curh4M=curHadr->Get4Momentum();
02178 if ( curh4M.e() > 0.) curHadr->SetPDGCode(22);
02179 else if( curh4M == nul4M)
02180 {
02181 G4QHadron* theLast = (*theResult)[nHd-1];
02182 if(theLast != curHadr)
02183 {
02184 curHadr->Set4Momentum(theLast->Get4Momentum());
02185 G4QPDGCode lQP=theLast->GetQPDG();
02186 if(lQP.GetPDGCode()!=10) curHadr->SetQPDG(lQP);
02187 else curHadr->SetQC(theLast->GetQC());
02188 }
02189 theResult->pop_back();
02190 --nHd;
02191 delete theLast;
02192 if(i == nHd-1) break;
02193 }
02194 }
02195
02196 else if(2 > 3)
02197 {
02198 for(G4int j=i+1; j<nHd; ++j)
02199 {
02200 G4int pPDG=(*theResult)[j]->GetPDGCode();
02201 if(hPDG==pPDG)
02202 {
02203 G4LorentzVector h4M=(*theResult)[i]->Get4Momentum();
02204 G4LorentzVector p4M=(*theResult)[j]->Get4Momentum();
02205 G4LorentzVector d4M=h4M+p4M;
02206 G4double E=d4M.m();
02207 if(hPDG==2212) E -= mProt+mProt;
02208 else E -= mNeut+mNeut;
02209 if(E < 140. && G4UniformRand() < .6)
02210 {
02211 G4int piPDG= 211;
02212 if(hPDG==2212) piPDG=-211;
02213 for(G4int k=0; k<nHd; ++k)
02214 {
02215 G4int mPDG=(*theResult)[k]->GetPDGCode();
02216
02217
02218 if(mPDG==111 || mPDG==piPDG)
02219 {
02220 G4LorentzVector m4M=(*theResult)[k]->Get4Momentum();
02221 G4double mN=mProt;
02222 G4int nPDG=2212;
02223 G4int tPDG=-211;
02224 if(hPDG==2212)
02225 {
02226 mN=mNeut;
02227 nPDG=2112;
02228 tPDG= 211;
02229 }
02230 G4double mPi=mPiZr;
02231 G4int sPDG=111;
02232 if(mPDG==111)
02233 {
02234 mPi=mPiCh;
02235 sPDG=tPDG;
02236 }
02237
02238
02239 G4double D=mPi+mN;
02240 G4LorentzVector t4M=m4M+h4M;
02241 G4LorentzVector n4M=h4M;
02242 G4double D2=D*D;
02243 G4double S=t4M.m2();
02244 if(S > D2) found= 1;
02245 else
02246 {
02247 t4M=m4M+p4M;
02248 n4M=p4M;
02249 S=t4M.m2();
02250 if(S > D2) found=-1;
02251 }
02252 if(found)
02253 {
02254 G4ThreeVector tV=t4M.vect()/t4M.e();
02255
02256 m4M.boost(-tV);
02257
02258 n4M.boost(-tV);
02259 G4double mPi2=mPi*mPi;
02260 G4double mN2=mN*mN;
02261 G4double C=S-mPi2-mN2;
02262 G4double p2=(C*C/4.-mPi2*mN2)/S;
02263 if(p2 < 0.) G4cout<<"-Warning-G4QFragment::Fragment: P2="<<p2<<G4endl;
02264 G4double pc2=m4M.vect().mag2();
02265
02266 G4double r=1.;
02267 if(pc2 < .00000000000001)
02268 G4cout<<"-Warning-G4QFragment::Fragment: PC2="<<pc2<<m4M<<G4endl;
02269 else r=std::sqrt(p2/pc2);
02270 m4M.setV(r*m4M.vect());
02271 m4M.setE(std::sqrt(mPi2+p2));
02272
02273
02274 n4M.setV(r*n4M.vect());
02275 n4M.setE(std::sqrt(mN2+p2));
02276 m4M.boost(tV);
02277 n4M.boost(tV);
02278 (*theResult)[k]->SetPDGCode(sPDG);
02279 (*theResult)[k]->Set4Momentum(m4M);
02280 if(found > 0)
02281 {
02282 (*theResult)[i]->SetPDGCode(nPDG);
02283 (*theResult)[i]->Set4Momentum(n4M);
02284 }
02285 else
02286 {
02287 (*theResult)[j]->SetPDGCode(nPDG);
02288 (*theResult)[j]->Set4Momentum(n4M);
02289 }
02290 break;
02291 }
02292 }
02293 }
02294 if(found) break;
02295 }
02296 }
02297 }
02298 }
02299
02300 G4int hChg=cHadr->GetCharge();
02301 if(hChg > maxChg)
02302 {
02303 maxChg = hChg;
02304 maxQC = cHadr->GetQC();
02305 #ifdef debug
02306 maxBN = cHadr->GetBaryonNumber();
02307 #endif
02308 }
02309 }
02310 G4QNucleus ResNucEnv(maxQC);
02311 #ifdef debug
02312 G4cout<<"G4QFragmentation::Fra: ResNucEnv with A="<<maxBN<<", Z="<<maxChg<<G4endl;
02313 #endif
02314
02315 G4LorentzVector sum(0.,0.,0.,0.);
02316 G4QContent sumQC(0,0,0,0,0,0);
02317 G4int sumCount=0;
02318 G4int nHadr=theResult->size();
02319 G4bool frag=false;
02320 if(nHadr>2) for(unsigned f=0; f<theResult->size(); f++)
02321 {
02322 G4int fBN=(*theResult)[f]->GetBaryonNumber();
02323 #ifdef debug
02324 G4int fPDG=(*theResult)[f]->GetPDGCode();
02325 G4LorentzVector fLV=(*theResult)[f]->Get4Momentum();
02326 G4cout<<"G4QFragmentation::Fra:"<<f<<",PDG="<<fPDG<<",fBN="<<fBN<<",f4M="<<fLV<<G4endl;
02327 #endif
02328 if(fBN>1)
02329 {
02330 frag=true;
02331 break;
02332 }
02333 }
02334 #ifdef debug
02335 G4cout<<"G4QFrag::Frag:=>Before Gamma Suppression<=, nH="<<nHadr<<",frag="<<frag<<G4endl;
02336 #endif
02337 if(nHadr>2 && frag) for(G4int h=nHadr-1; h>=0; h--)
02338 {
02339 G4QHadron* curHadr = (*theResult)[h];
02340 G4int hF = curHadr->GetNFragments();
02341 G4int hPDG = curHadr->GetPDGCode();
02342 G4LorentzVector h4M=curHadr->Get4Momentum();
02343 if(hPDG==89999003||hPDG==90002999)G4cout<<"-Warning-G4QFr::Fr:nD-/pD++="<<hPDG<<G4endl;
02344 #ifdef debug
02345 G4cout<<"G4QFragmentation::Fragm: h#"<<h<<", hPDG="<<hPDG<<", hNFrag="<<hF<<G4endl;
02346 #endif
02347 G4int hChg = curHadr->GetCharge();
02348 G4bool UCB = false;
02349 if(hChg > 0 && hPDG!=321)
02350 {
02351 G4int hBN = curHadr->GetBaryonNumber();
02352 G4double hCB=ResNucEnv.CoulombBarrier(hChg,hBN);
02353 if(h4M.e()-h4M.m() < hCB) UCB = true;
02354 }
02355 if(hF || hPDG==22 || UCB)
02356 {
02357 G4int last=theResult->size()-1;
02358 G4QHadron* theLast = (*theResult)[last];
02359 if(hPDG==22 || UCB)
02360 {
02361 sum+=h4M;
02362 sumCount++;
02363 if(UCB) sumQC+=curHadr->GetQC();
02364 #ifdef debug
02365 G4cout<<"G4QFragmentation::Frag: gam4M="<<h4M<<" is added to s4M="<<sum<<G4endl;
02366 #endif
02367 }
02368 nHadr = static_cast<G4int>(theResult->size())-1;
02369 if(h < last)
02370 {
02371 curHadr->SetNFragments(0);
02372 curHadr->Set4Momentum(theLast->Get4Momentum());
02373 G4QPDGCode lQP=theLast->GetQPDG();
02374 if(lQP.GetPDGCode()!=10) curHadr->SetQPDG(lQP);
02375 else curHadr->SetQC(theLast->GetQC());
02376 #ifdef debug
02377 G4cout<<"G4QFragmentation::Fragment: Exchange with the last is done"<<G4endl;
02378 #endif
02379 }
02380 theResult->pop_back();
02381 delete theLast;
02382 #ifdef debug
02383 G4cout<<"G4QFragmentation::Fragment: The last is compessed"<<G4endl;
02384 #endif
02385 }
02386 }
02387 #ifdef debug
02388 G4cout<<"G4QFragment::Frag: nH="<<nHadr<<"="<<theResult->size()<<", sum="<<sum<<G4endl;
02389 #endif
02390 if(nHadr > 1) for(unsigned hdr=0; hdr<theResult->size()-1; hdr++)
02391 {
02392 G4QHadron* curHadr = (*theResult)[hdr];
02393 #ifdef debug
02394 G4cout<<"G4QFrag::Frag: h#"<<hdr<<"<"<<nHadr<<", hPDG="<<curHadr->GetPDGCode()<<G4endl;
02395 #endif
02396 G4QHadron* theLast = (*theResult)[theResult->size()-1];
02397 G4int hB = curHadr->GetBaryonNumber();
02398 G4int lB = theLast->GetBaryonNumber();
02399 #ifdef debug
02400 G4cout<<"G4QFra::Fra:hBN="<<hB<<"<lBN="<<lB<<",lstPDG="<<theLast->GetPDGCode()<<G4endl;
02401 #endif
02402 if(lB < hB)
02403 {
02404 G4QPDGCode hQPDG = curHadr->GetQPDG();
02405 G4LorentzVector h4m= curHadr->Get4Momentum();
02406 curHadr->Set4Momentum(theLast->Get4Momentum());
02407 G4QPDGCode lQP=theLast->GetQPDG();
02408 if(lQP.GetPDGCode()!=10) curHadr->SetQPDG(lQP);
02409 else curHadr->SetQC(theLast->GetQC());
02410 theLast->Set4Momentum(h4m);
02411 theLast->SetQPDG(hQPDG);
02412 }
02413 }
02414 nHadr=theResult->size();
02415 if(sumCount)
02416 {
02417 G4QHadron* theLast = (*theResult)[nHadr-1];
02418 G4int nucEnvBN=theLast->GetBaryonNumber();
02419 if ( nucEnvBN > 0 )
02420 {
02421 G4QHadron* theNew = new G4QHadron(theLast);
02422 #ifdef debug
02423 G4cout<<"G4QFra::Fra:BeforeLastSub,n="<<nHadr<<",PDG="<<theNew->GetPDGCode()<<G4endl;
02424 #endif
02425 theResult->pop_back();
02426 delete theLast;
02427 nHadr--;
02428 G4QContent newQC=theNew->GetQC();
02429 G4LorentzVector new4M=theNew->Get4Momentum();
02430 #ifdef debug
02431 G4cout<<"G4QFra::Fra:gSum4M="<<sum<<" is added to "<<new4M<<", QC="<<newQC<<G4endl;
02432 #endif
02433 G4LorentzVector exRes4M = new4M + sum;
02434 G4QContent exResQC = newQC + sumQC;
02435 theNew->Set4Momentum(exRes4M);
02436 theNew->SetQC(exResQC);
02437 #ifdef debug
02438 G4cout<<"G4QFra::Fra:BeforeEvap, n="<<nHadr<<", nPDG="<<theNew->GetPDGCode()<<G4endl;
02439 #endif
02440 EvaporateResidual(theNew);
02441 nHadr=theResult->size();
02442 }
02443 else G4cout<<"-Warning-G4QFragmentation::Fragment:RA="<<nucEnvBN<<",E/M cons?"<<G4endl;
02444 }
02445
02446 return theResult;
02447 }
02448
02449 void G4QFragmentation::Breeder()
02450 {
02451 static const G4double eps = 0.001;
02452
02453 static const G4LorentzVector vac4M(0.,0.,0.,0.);
02454
02455
02456
02457 #ifdef edebug
02458 G4LorentzVector totLS4M=theNucleus.Get4Momentum();
02459 G4int totChg=theNucleus.GetZ();
02460 G4int totBaN=theNucleus.GetA();
02461 G4int nStri=strings.size();
02462 G4cout<<"-EMCLS-G4QFr::Breed: CHECKRecovery #ofS="<<nStri<<",N4M(E=M)="<<totLS4M<<G4endl;
02463 for(G4int i=0; i<nStri; i++)
02464 {
02465 G4LorentzVector strI4M=strings[i]->Get4Momentum();
02466 totLS4M+=strI4M;
02467 G4int sChg=strings[i]->GetCharge();
02468 totChg+=sChg;
02469 G4int sBaN=strings[i]->GetBaryonNumber();
02470 totBaN+=sBaN;
02471 G4cout<<"-EMCLS-G4QFragm::Breeder: St#"<<i<<", 4M="<<strI4M<<", M="<<strI4M.m()
02472 <<", C="<<sChg<<", B="<<sBaN<<G4endl;
02473 }
02474 #endif
02475 G4int nOfStr=strings.size();
02476 #ifdef debug
02477 G4cout<<"G4QFragmentation::Breeder: BeforeFragmentation, #OfStr="<<nOfStr<<G4endl;
02478 #endif
02479 G4LorentzVector ft4M(0.,0.,0.,0.);
02480 G4QContent ftQC(0,0,0,0,0,0);
02481 G4bool ftBad=false;
02482 for(G4int i=0; i < nOfStr; ++i)
02483 {
02484 G4QString* crStr=strings[i];
02485 G4LorentzVector pS4M=crStr->Get4Momentum();
02486 ft4M+=pS4M;
02487 G4QContent pSQC=crStr->GetQC();
02488 ftQC+=pSQC;
02489 if(pS4M.m2() < 0.) ftBad=true;
02490 #ifdef debug
02491 G4cout<<">G4QFrag::Br:1stTest,S#"<<i<<",P="<<crStr<<",4M="<<pS4M<<",QC="<<pSQC<<G4endl;
02492 #endif
02493 }
02494 if(ftBad)
02495 {
02496 G4Quasmon* stringQuasmon = new G4Quasmon(ftQC, ft4M);
02497 #ifdef debug
02498 G4cout<<"->G4QFragmentation::Breeder:*TotQ*,QC="<<ftQC<<",4M="<<ft4M<<ft4M.m()<<G4endl;
02499 #endif
02500 theQuasmons.push_back(stringQuasmon);
02501 G4LorentzVector r4M=theNucleus.Get4Momentum();
02502 G4int rPDG=theNucleus.GetPDG();
02503 G4QHadron* resNuc = new G4QHadron(rPDG,r4M);
02504 theResult->push_back(resNuc);
02505 return;
02506 }
02507 for (G4int astring=0; astring < nOfStr; astring++)
02508 {
02509 #ifdef edebug
02510 G4LorentzVector sum=theNucleus.Get4Momentum();
02511 G4int rChg=totChg-theNucleus.GetZ();
02512 G4int rBaN=totBaN-theNucleus.GetA();
02513 G4int nOfHadr=theResult->size();
02514 G4cout<<"-EMCLS-G4QFragmentation::Breeder:#ofSt="<<nOfStr<<",#ofHad="<<nOfHadr<<G4endl;
02515 for(G4int i=astring; i<nOfStr; i++)
02516 {
02517 G4LorentzVector strI4M=strings[i]->Get4Momentum();
02518 sum+=strI4M;
02519 G4int sChg=strings[i]->GetCharge();
02520 rChg-=sChg;
02521 G4int sBaN=strings[i]->GetBaryonNumber();
02522 rBaN-=sBaN;
02523 G4cout<<"-EMCLS-G4QF::Breed:S#"<<i<<",4M="<<strI4M<<",C="<<sChg<<",B="<<sBaN<<G4endl;
02524 }
02525 for(G4int i=0; i<nOfHadr; i++)
02526 {
02527 G4LorentzVector hI4M=(*theResult)[i]->Get4Momentum();
02528 sum+=hI4M;
02529 G4int hChg=(*theResult)[i]->GetCharge();
02530 rChg-=hChg;
02531 G4int hBaN=(*theResult)[i]->GetBaryonNumber();
02532 rBaN-=hBaN;
02533 G4cout<<"-EMCLS-G4QFr::Breed: H#"<<i<<",4M="<<hI4M<<",C="<<hChg<<",B="<<hBaN<<G4endl;
02534 }
02535 G4cout<<"....-EMCLS-G4QFrag::Br:r4M="<<sum-totLS4M<<",rC="<<rChg<<",rB="<<rBaN<<G4endl;
02536 #endif
02537 G4QString* curString=strings[astring];
02538 if(!curString->GetDirection()) continue;
02539 G4int curStrChg = curString->GetCharge();
02540 G4int curStrBaN = curString->GetBaryonNumber();
02541 G4LorentzVector curString4M = curString->Get4Momentum();
02542 #ifdef debug
02543 G4cout<<"=--=>G4QFragmentation::Breeder: String#"<<astring<<",s4M/m="<<curString4M
02544 <<curString4M.m()<<", LPart="<<curString->GetLeftParton()->GetPDGCode()
02545 <<", RPart="<<curString->GetRightParton()->GetPDGCode()<<G4endl;
02546 #endif
02547 G4QHadronVector* theHadrons = 0;
02548 theHadrons=curString->FragmentString(true);
02549 if (!theHadrons)
02550 {
02551
02552 G4QParton* cLeft=curString->GetLeftParton();
02553 G4QParton* cRight=curString->GetRightParton();
02554 G4int sPDG=cLeft->GetPDGCode();
02555 G4int nPDG=cRight->GetPDGCode();
02556 G4int LT=cLeft->GetType();
02557 G4int RT=cRight->GetType();
02558 G4int LS=LT+RT;
02559 if(LT==2 && RT==2)
02560 {
02561 #ifdef debug
02562 G4cout<<"G4QFragmentation::Breeder:TryReduceString, L="<<sPDG<<",R="<<nPDG<<G4endl;
02563 #endif
02564 if( cLeft->ReduceDiQADiQ(cLeft, cRight) )
02565 {
02566 LT=1;
02567 RT=1;
02568 LS=2;
02569 sPDG=cLeft->GetPDGCode();
02570 nPDG=cRight->GetPDGCode();
02571 #ifdef debug
02572 G4cout<<"G4QFragmentation::Breeder:AfterReduction,L="<<sPDG<<",R="<<nPDG<<G4endl;
02573 #endif
02574 theHadrons=curString->FragmentString(true);
02575 cLeft=curString->GetLeftParton();
02576 cRight=curString->GetRightParton();
02577 #ifdef debug
02578 G4cout<<"G4QFrag::Breed:L="<<cLeft->Get4Momentum()<<",R="<<cRight->Get4Momentum()
02579 <<G4endl;
02580 #endif
02581 }
02582 #ifdef debug
02583 else G4cout<<"^G4QFragmentation::Breeder: DQ-aDQ reduction to Q-aQ Failed"<<G4endl;
02584 #endif
02585 }
02586 #ifdef debug
02587 G4cout<<"G4QFrag::Breed:AfterRedAttempt, theH="<<theHadrons<<", L4M="
02588 <<cLeft->Get4Momentum()<<", R4M="<<cRight->Get4Momentum()<<G4endl;
02589 #endif
02590 unsigned next=astring+1;
02591 if (!theHadrons)
02592 {
02593 G4int fusionDONE=0;
02594 if(next < strings.size())
02595 {
02596 G4int fustr=0;
02597 G4int swap=0;
02598 G4double Vmin=DBL_MAX;
02599 G4int dPDG=nPDG;
02600 G4int qPDG=sPDG;
02601 if(dPDG<-99 || (dPDG>0&&dPDG<7) || qPDG>99 || (qPDG<0 && qPDG>-7))
02602 {
02603 swap=qPDG;
02604 qPDG=dPDG;
02605 dPDG=swap;
02606 }
02607 if(dPDG>99) dPDG/=100;
02608 if(qPDG<-99) qPDG=-(-qPDG)/100;
02609 #ifdef debug
02610 G4cout<<"G4QFrag::Breed:TryFuseStringS, q="<<qPDG<<", a="<<dPDG<<", n="<<next
02611 <<G4endl;
02612 #endif
02613 G4ThreeVector curV=curString4M.vect()/curString4M.e();
02614 G4int reduce=0;
02615 G4int restr=0;
02616 G4int MPS=0;
02617 for (restr=next; restr < nOfStr; restr++)
02618 {
02619 reduce=0;
02620 G4QString* reString=strings[restr];
02621 G4QParton* Left=reString->GetLeftParton();
02622 G4QParton* Right=reString->GetRightParton();
02623 G4int uPDG=Left->GetPDGCode();
02624 G4int mPDG=Right->GetPDGCode();
02625 G4int PLT =Left->GetType();
02626 G4int PRT =Right->GetType();
02627 G4int aPDG=mPDG;
02628 G4int rPDG=uPDG;
02629 if(aPDG<-99 || (aPDG>0 && aPDG<7) || rPDG>99 || (rPDG<0 && rPDG>-7))
02630 {
02631 swap=rPDG;
02632 rPDG=aPDG;
02633 aPDG=swap;
02634 }
02635 if(aPDG > 99) aPDG/=100;
02636 if(rPDG <-99) rPDG=-(-rPDG)/100;
02637
02638 #ifdef debug
02639 G4cout<<"G4QFragm::Breed: TryReduce #"<<restr<<",q="<<rPDG<<",a="<<aPDG<<G4endl;
02640 #endif
02641 if(LT==2 && RT==2 && PLT==2 && PRT==2)
02642 {
02643 G4int cQ1=(-qPDG)/10;
02644 G4int cQ2=(-qPDG)%10;
02645 G4int cA1=dPDG/10;
02646 G4int cA2=dPDG%10;
02647 G4int pQ1=(-rPDG)/10;
02648 G4int pQ2=(-rPDG)%10;
02649 G4int pA1=aPDG/10;
02650 G4int pA2=aPDG%10;
02651 #ifdef debug
02652 G4cout<<"G4QFragment::Breeder: cQ="<<cQ1<<","<<cQ2<<", cA="<<cA1<<","<<cA2
02653 <<", pQ="<<pQ1<<","<<pQ2<<", pA="<<pA1<<","<<pA2<<G4endl;
02654 #endif
02655 G4bool iQA = (cA1==pQ1 || cA1==pQ2 || cA2==pQ1 || cA2==pQ2);
02656 G4bool iAQ = (cQ1==pA1 || cQ1==pA2 || cQ2==pA1 || cQ2==pA2);
02657 if(iQA) reduce++;
02658 if(iAQ) reduce++;
02659 if (reduce==2)
02660 {
02661 if(sPDG>0 && uPDG<0)
02662 {
02663 std::pair<G4int,G4int> resLL=ReducePair(sPDG/100, (-uPDG)/100);
02664 G4int newCL=resLL.first;
02665 G4int newPL=resLL.second;
02666 if(!newCL || !newPL)
02667 {
02668 G4cerr<<"*G4QFragmentation::Breeder:CL="<<newCL<<",PL="<<newPL<<G4endl;
02669 G4Exception("G4QFragmentation::Breeder:","72",FatalException,"2-LL-");
02670 }
02671 std::pair<G4int,G4int> resRR=ReducePair((-nPDG)/100, mPDG/100);
02672 G4int newCR=resRR.first;
02673 G4int newPR=resRR.second;
02674 if(!newCR || !newPR)
02675 {
02676 G4cerr<<"*G4QFragmentation::Breeder:CR="<<newCR<<",PR="<<newPR<<G4endl;
02677 G4Exception("G4QFragmentation::Breeder:","72",FatalException,"2-RR-");
02678 }
02679 cLeft->SetPDGCode(newCL);
02680 cRight->SetPDGCode(-newCR);
02681 Left->SetPDGCode(-newPL);
02682 Right->SetPDGCode(newPR);
02683 break;
02684 }
02685 else if(sPDG<0 && uPDG>0)
02686 {
02687 std::pair<G4int,G4int> resLL=ReducePair((-sPDG)/100, uPDG/100);
02688 G4int newCL=resLL.first;
02689 G4int newPL=resLL.second;
02690 if(!newCL || !newPL)
02691 {
02692 G4cerr<<"*G4QFragmentation::Breeder:CL="<<newCL<<",PL="<<newPL<<G4endl;
02693 G4Exception("G4QFragmentation::Breeder:","72",FatalException,"2-LL+");
02694 }
02695 std::pair<G4int,G4int> resRR=ReducePair(nPDG/100, (-mPDG)/100);
02696 G4int newCR=resRR.first;
02697 G4int newPR=resRR.second;
02698 if(!newCR || !newPR)
02699 {
02700 G4cerr<<"*G4QFragmentation::Breeder:CR="<<newCR<<",PR="<<newPR<<G4endl;
02701 G4Exception("G4QFragmentation::Breeder:","72",FatalException,"2-RR+");
02702 }
02703 cLeft->SetPDGCode(-newCL);
02704 cRight->SetPDGCode(newCR);
02705 Left->SetPDGCode(newPL);
02706 Right->SetPDGCode(-newPR);
02707 break;
02708 }
02709 else if(sPDG>0 && mPDG<0)
02710 {
02711 std::pair<G4int,G4int> resLL=ReducePair(sPDG/100, (-mPDG)/100);
02712 G4int newCL=resLL.first;
02713 G4int newPR=resLL.second;
02714 if(!newCL || !newPR)
02715 {
02716 G4cerr<<"*G4QFragmentation::Breeder:CL="<<newCL<<",PR="<<newPR<<G4endl;
02717 G4Exception("G4QFragmentation::Breeder:","72",FatalException,"2-LR");
02718 }
02719 std::pair<G4int,G4int> resRR=ReducePair((-nPDG)/100, uPDG/100);
02720 G4int newCR=resRR.first;
02721 G4int newPL=resRR.second;
02722 if(!newCR || !newPR)
02723 {
02724 G4cerr<<"*G4QFragmentation::Breeder:CR="<<newCR<<",PL="<<newPL<<G4endl;
02725 G4Exception("G4QFragmentation::Breeder:","72",FatalException,"2-LR");
02726 }
02727 cLeft->SetPDGCode(newCL);
02728 cRight->SetPDGCode(-newCR);
02729 Left->SetPDGCode(newPL);
02730 Right->SetPDGCode(-newPR);
02731 break;
02732 }
02733 else
02734 {
02735 std::pair<G4int,G4int> resLL=ReducePair((-sPDG)/100, mPDG/100);
02736 G4int newCL=resLL.first;
02737 G4int newPR=resLL.second;
02738 if(!newCL || !newPR)
02739 {
02740 G4cerr<<"*G4QFragmentation::Breeder:CL="<<newCL<<",PR="<<newPR<<G4endl;
02741 G4Exception("G4QFragmentation::Breeder:","72",FatalException,"2-RL");
02742 }
02743 std::pair<G4int,G4int> resRR=ReducePair(nPDG/100, (-uPDG)/100);
02744 G4int newCR=resRR.first;
02745 G4int newPL=resRR.second;
02746 if(!newCR || !newPR)
02747 {
02748 G4cerr<<"*G4QFragmentation::Breeder:CR="<<newCR<<",PL="<<newPL<<G4endl;
02749 G4Exception("G4QFragmentation::Breeder:","72",FatalException,"2-RL");
02750 }
02751 cLeft->SetPDGCode(-newCL);
02752 cRight->SetPDGCode(newCR);
02753 Left->SetPDGCode(-newPL);
02754 Right->SetPDGCode(newPR);
02755 break;
02756 }
02757 }
02758 }
02759
02760
02761 #ifdef debug
02762 G4cout<<"G4QFragm::Breed:TryFuse/w #"<<restr<<",q="<<rPDG<<",a="<<aPDG<<G4endl;
02763 #endif
02764 G4int PLS=PLT+PRT;
02765 if( (LS==2 && PLS==2) ||
02766 ( ( (LS==2 && PLS==3) || (LS==3 && PLS==2) ) &&
02767 ( (aPDG> 7 && (-dPDG==aPDG/10 || -dPDG==aPDG%10) ) ||
02768 (dPDG> 7 && (-aPDG==dPDG/10 || -aPDG==dPDG%10) ) ||
02769 (rPDG<-7 && (qPDG==(-rPDG)/10 || qPDG==(-rPDG)%10) ) ||
02770 (qPDG<-7 && (rPDG==(-qPDG)/10 || rPDG==(-qPDG)%10) )
02771
02772 )
02773 ) ||
02774 ( ( LS==2 && PLS==4 &&
02775 (aPDG> 7 && (-dPDG == aPDG/10 || -dPDG == aPDG%10) ) &&
02776 (rPDG<-7 && (qPDG==(-rPDG)/10 || qPDG==(-rPDG)%10) )
02777 ) ||
02778 ( LS==4 && PLS==2 &&
02779 (dPDG> 7 && (-aPDG == dPDG/10 || -aPDG == dPDG%10) ) &&
02780 (qPDG<-7 && (rPDG==(-qPDG)/10 || rPDG==(-qPDG)%10) )
02781 )
02782 ) ||
02783 ( LS==3 && PLS==3 &&
02784 ( (aPDG> 7 && (-dPDG == aPDG/10 || -dPDG == aPDG%10) &&
02785 qPDG<-7 && (rPDG==(-qPDG)/10 || rPDG==(-qPDG)%10)
02786 ) ||
02787 (dPDG> 7 && (-aPDG == dPDG/10 || -aPDG == dPDG%10) &&
02788 rPDG<-7 && (dPDG==(-rPDG)/10 || dPDG==(-rPDG)%10)
02789 )
02790 )
02791 )
02792 )
02793 {
02794 G4LorentzVector reString4M = reString->Get4Momentum();
02795 G4ThreeVector reV = reString4M.vect()/reString4M.e();
02796 G4double dV=(curV-reV).mag2();
02797 #ifdef debug
02798 G4cout<<"G4QFragmentation::Breeder: StringCand#"<<restr<<", q="<<rPDG<<", a="
02799 <<aPDG<<", L="<<uPDG<<", R="<<mPDG<<",dV="<<dV<<G4endl;
02800 #endif
02801 if(dV < Vmin)
02802 {
02803 Vmin=dV;
02804 fustr=restr;
02805 MPS=PLS;
02806 }
02807 }
02808 }
02809 if(reduce==2)
02810 {
02811 #ifdef debug
02812 G4cout<<"-G4QFragmentation::Breeder:Reduced #"<<astring<<" & #"<<restr<<G4endl;
02813 #endif
02814 astring--;
02815 continue;
02816 }
02817 if(fustr)
02818 {
02819 #ifdef debug
02820 G4cout<<"G4QFragm::Breeder: StPartner#"<<fustr<<",LT="<<LT<<",RT="<<RT<<G4endl;
02821 #endif
02822 G4QString* fuString=strings[fustr];
02823 G4QParton* Left=fuString->GetLeftParton();
02824 G4QParton* Right=fuString->GetRightParton();
02825 G4int uPDG=Left->GetPDGCode();
02826 G4int mPDG=Right->GetPDGCode();
02827 G4int rPDG=uPDG;
02828 G4int aPDG=mPDG;
02829 if(aPDG<-99 || (aPDG>0 && aPDG<7) || rPDG>99 || (rPDG<0 && rPDG>-7))
02830 {
02831 swap=rPDG;
02832 rPDG=aPDG;
02833 aPDG=swap;
02834 }
02835 if(aPDG > 99) aPDG/=100;
02836 if(rPDG <-99) rPDG=-(-rPDG)/100;
02837
02838 G4LorentzVector resL4M(0.,0.,0.,0.);
02839 G4LorentzVector resR4M(0.,0.,0.,0.);
02840 G4LorentzVector L4M=cLeft->Get4Momentum();
02841 G4LorentzVector R4M=cRight->Get4Momentum();
02842 #ifdef debug
02843 G4cout<<"G4QFragmentation::Breeder:BeforeFuDir,sL="<<sPDG<<",nR="<<nPDG<<",uL="
02844 <<uPDG<<",mR="<<mPDG<<",L4M="<<L4M<<",R4M="<<R4M<<G4endl;
02845 #endif
02846 G4LorentzVector PL4M=Left->Get4Momentum();
02847 G4LorentzVector PR4M=Right->Get4Momentum();
02848 fusionDONE=AnnihilationOrder(LS, MPS, uPDG, mPDG, sPDG, nPDG);
02849 if (fusionDONE>0)
02850 {
02851 if(fusionDONE>1)
02852 {
02853 if ( (uPDG<0 || nPDG<0) && -uPDG==nPDG ) Left->SetPDGCode(sPDG);
02854 else if( (mPDG<0 || sPDG<0) && -mPDG==sPDG ) Right->SetPDGCode(nPDG);
02855 else if( (uPDG<0 || sPDG<0) && -uPDG==sPDG ) Left->SetPDGCode(nPDG);
02856 else if( (mPDG<0 || nPDG<0) && -mPDG==nPDG ) Right->SetPDGCode(sPDG);
02857 }
02858 {
02859 Left->SetPDGCode(SumPartonPDG(uPDG, sPDG));
02860 Right->SetPDGCode(SumPartonPDG(mPDG, nPDG));
02861 }
02862 Left->Set4Momentum(L4M+PL4M);
02863 Right->Set4Momentum(R4M+PR4M);
02864 #ifdef debug
02865 G4cout<<"G4QFragmentation::Breeder:LL/RR s4M="<<fuString->Get4Momentum()
02866 <<",S="<<L4M+PL4M+R4M+PR4M<<", L="<<Left->Get4Momentum()<<", R="
02867 <<Right->Get4Momentum()<<G4endl;
02868 #endif
02869 }
02870 else if(fusionDONE<0)
02871 {
02872 Left->SetPDGCode(SumPartonPDG(uPDG, nPDG));
02873 Left->Set4Momentum(L4M+PR4M);
02874 Right->SetPDGCode(SumPartonPDG(mPDG, sPDG));
02875 Right->Set4Momentum(R4M+PL4M);
02876 #ifdef debug
02877 G4cout<<"G4QFragmentation::Breeder:LR/RL s4M="<<fuString->Get4Momentum()
02878 <<",S="<<L4M+PL4M+R4M+PR4M<<", L="<<Left->Get4Momentum()<<", R="
02879 <<Right->Get4Momentum()<<G4endl;
02880 #endif
02881 }
02882 #ifdef debug
02883 else G4cout<<"-Warning-G4QFragmentation::Breeder: WrongStringFusion"<<G4endl;
02884 #endif
02885 #ifdef edebug
02886 G4cout<<"#EMC#G4QFragmentation::Breeder:StringFused,F="<<fusionDONE<<",L="<<L4M
02887 <<",R="<<R4M<<",pL="<<PL4M<<",pR="<<PR4M<<",nL="<<Left->Get4Momentum()
02888 <<",nR="<<Right->Get4Momentum()<<",S="<<fuString->Get4Momentum()<<G4endl;
02889 #endif
02890 if(fusionDONE)
02891 {
02892 #ifdef debug
02893 G4cout<<"###G4QFragmentation::Breeder: Str#"<<astring<<" fused/w Str#"<<fustr
02894 <<"->"<<fuString->Get4Momentum()<<fuString->Get4Momentum().m()<<G4endl;
02895 #endif
02896 continue;
02897 }
02898 }
02899 #ifdef debug
02900 else
02901 {
02902
02903 G4cout<<"**G4QFragmentation::Breeder:*NoPart*M="<<curString->Get4Momentum().m()
02904 <<", F="<<fusionDONE<<", LPDG="<<curString->GetLeftParton()->GetPDGCode()
02905 <<", RPDG="<<curString->GetRightParton()->GetPDGCode()<<G4endl;
02906 }
02907 #endif
02908 }
02909 if(!fusionDONE)
02910 {
02911 G4int nHadr=theResult->size();
02912 #ifdef debug
02913 G4cout<<"+++4QFragmentation::Breeder:*TryHadr* M="<<curString->Get4Momentum().m()
02914 <<", nH="<<nHadr<<", LPDG="<<curString->GetLeftParton()->GetPDGCode()
02915 <<", RPDG="<<curString->GetRightParton()->GetPDGCode()<<G4endl;
02916 #endif
02917
02918 while( (nHadr=theResult->size()) && !theHadrons)
02919 {
02920 #ifdef edebug
02921 for(G4int i=0; i<nHadr; i++)
02922 {
02923 G4LorentzVector h4M=(*theResult)[i]->Get4Momentum();
02924 G4int hPDG=(*theResult)[i]->GetPDGCode();
02925 G4QContent hQC=(*theResult)[i]->GetQC();
02926 G4cout<<"-EMC-G4QFrag::Breed:H#"<<i<<",4M="<<h4M<<hQC<<",PDG="<<hPDG<<G4endl;
02927 }
02928 #endif
02929 G4int fusDONE=0;
02930 G4int fuhad=-1;
02931 G4int newPDG=0;
02932 G4int secPDG=0;
02933 G4double maM2=-DBL_MAX;
02934 G4LorentzVector selH4M(0.,0.,0.,0.);
02935 G4QHadron* selHP=0;
02936 G4QString* cString=strings[astring];
02937 G4LorentzVector cString4M = cString->Get4Momentum();
02938 cLeft=cString->GetLeftParton();
02939 cRight=cString->GetRightParton();
02940 G4int sumT=cLeft->GetType()+cRight->GetType();
02941 sPDG=cLeft->GetPDGCode();
02942 nPDG=cRight->GetPDGCode();
02943 G4int cMax=0;
02944 for (G4int reh=0; reh < nHadr; reh++)
02945 {
02946 G4QHadron* curHadr=(*theResult)[reh];
02947 G4int curPDG=curHadr->GetPDGCode();
02948 G4QContent curQC=curHadr->GetQC();
02949 if(curPDG==331 && sPDG!=3 && nPDG!=3 && sPDG!=-3 && nPDG!=-3)
02950 {
02951 if(sPDG==2 || sPDG==-2 || nPDG==2 || nPDG==-2)
02952 curQC=G4QContent(0,1,0,0,1,0);
02953 else curQC=G4QContent(1,0,0,1,0,0);
02954 }
02955 else if(curPDG==221 && sPDG!=2 && nPDG!=2 && sPDG!=-2 && nPDG!=-2)
02956 curQC=G4QContent(1,0,0,1,0,0);
02957 else if(curPDG==111 && sPDG!=1 && nPDG!=1 && sPDG!=-1 && nPDG!=-1)
02958 curQC=G4QContent(0,1,0,0,1,0);
02959 G4int partPDG1=curQC.AddParton(sPDG);
02960 G4int partPDG2=curQC.AddParton(nPDG);
02961 #ifdef debug
02962 G4cout<<"G4QFragmentation::Breeder: Hadron # "<<reh<<", QC="<<curQC
02963 <<", P1="<<partPDG1<<", P2="<<partPDG2<<G4endl;
02964 #endif
02965 if(partPDG1 || partPDG2)
02966 {
02967 G4int cCur=1;
02968 if(sumT>3 && partPDG1 && partPDG2) cCur=2;
02969 G4LorentzVector curHadr4M = curHadr->Get4Momentum();
02970 G4double M2=(cString4M+curHadr4M).m2();
02971 #ifdef debug
02972 G4cout<<"G4QFragmentation::Breeder:*IN*Hadron#"<<reh<<",M2="<<M2<<G4endl;
02973 #endif
02974 if( (sumT<4 || cCur>=cMax) && M2 > maM2)
02975 {
02976 maM2=M2;
02977 fuhad=reh;
02978 if(partPDG1)
02979 {
02980 fusDONE= 1;
02981 newPDG=partPDG1;
02982 if(partPDG2)
02983 {
02984 secPDG=partPDG2;
02985 cMax=cCur;
02986 }
02987 }
02988 else
02989 {
02990 fusDONE=-1;
02991 newPDG=partPDG2;
02992 }
02993 #ifdef debug
02994 G4cout<<"G4QFrag::Br:*Selected*,P1="<<partPDG1<<",P2="<<partPDG2<<G4endl;
02995 #endif
02996 selH4M=curHadr4M;
02997 selHP=curHadr;
02998 }
02999 }
03000 }
03001 #ifdef debug
03002 G4cout<<"G4QFragmentation::Breeder: fuh="<<fuhad<<",fus="<<fusDONE<<G4endl;
03003 #endif
03004 if(fuhad>-1)
03005 {
03006 if (fusDONE>0)
03007 {
03008 cLeft->SetPDGCode(newPDG);
03009 G4LorentzVector newLeft=cLeft->Get4Momentum()+selH4M;
03010 cLeft->Set4Momentum(newLeft);
03011 if(secPDG && cMax>1)
03012 {
03013 #ifdef debug
03014 G4cout<<"G4QFragm::Br:TryReduce,nPDG="<<newPDG<<",sPDG="<<secPDG<<G4endl;
03015 #endif
03016 cLeft->ReduceDiQADiQ(cLeft, cRight);
03017 }
03018 #ifdef debug
03019 G4cout<<"G4QFragmentation::Breeder: Left, s4M="<<curString->Get4Momentum()
03020 <<", L4M="<<newLeft<<", R4M="<<cRight->Get4Momentum()<<", h4M="
03021 <<selH4M<<", newL="<<newPDG<<", oldR="<<cRight->GetPDGCode()<<G4endl;
03022 #endif
03023 }
03024 else if(fusDONE<0)
03025 {
03026 cRight->SetPDGCode(newPDG);
03027 G4LorentzVector newRight=cRight->Get4Momentum()+selH4M;
03028 cRight->Set4Momentum(newRight);
03029 #ifdef debug
03030 G4cout<<"G4QFragmentation::Breeder: Right, s4M="<<curString->Get4Momentum()
03031 <<", L4M="<<cLeft->Get4Momentum()<<", R4M="<<newRight<<", h4M="
03032 <<selH4M<<", newR="<<newPDG<<", oldL="<<cLeft->GetPDGCode()<<G4endl;
03033 #endif
03034 }
03035 #ifdef debug
03036 else G4cout<<"-G4QFragmentation::Breeder: Wrong String+HadronFusion"<<G4endl;
03037 #endif
03038 #ifdef debug
03039 if(fusDONE) G4cout<<"####G4QFragmentation::Breeder: String #"<<astring
03040 <<" is fused with Hadron #"<<fuhad
03041 <<", new4Mom="<<curString->Get4Momentum()
03042 <<", M2="<<curString->Get4Momentum().m2()
03043 <<", QC="<<curString->GetQC()<<G4endl;
03044 #endif
03045 }
03046 else
03047 {
03048 #ifdef debug
03049 G4cout<<"**G4QFragmentation::Breeder:*NoH*,M="<<curString->Get4Momentum().m()
03050 <<", LPDG="<<curString->GetLeftParton()->GetPDGCode()
03051 <<", RPDG="<<curString->GetRightParton()->GetPDGCode()<<G4endl;
03052
03053
03054 #endif
03055 break;
03056 }
03057
03058 #ifdef debug
03059 G4cout<<"G4QFragmentation::Breeder: before HR, nH="<<theResult->size()<<G4endl;
03060 G4int icon=0;
03061 #endif
03062 G4QHadronVector::iterator ih;
03063 #ifdef debug
03064 G4bool found=false;
03065 #endif
03066 for(ih = theResult->begin(); ih != theResult->end(); ih++)
03067 {
03068 #ifdef debug
03069 G4cout<<"G4QFrag::Breeder:#"<<icon<<", i="<<(*ih)<<", sH="<<selHP<<G4endl;
03070 icon++;
03071 #endif
03072 if((*ih)==selHP)
03073 {
03074 #ifdef debug
03075 G4cout<<"G4QFragm::Breed: *HadronFound*,PDG="<<selHP->GetPDGCode()<<G4endl;
03076 #endif
03077 G4LorentzVector p4M=selHP->Get4Momentum();
03078
03079 curString4M+=p4M;
03080 G4int Chg=selHP->GetCharge();
03081 G4int BaN=selHP->GetBaryonNumber();
03082 curStrChg+=Chg;
03083 curStrBaN+=BaN;
03084 #ifdef edebug
03085 G4cout<<"-EMC->>G4QFragmentation::Breeder: S+=H, 4M="<<curString4M<<", M="
03086 <<curString4M.m()<<", Charge="<<curStrChg<<", B="<<curStrBaN<<G4endl;
03087 #endif
03088 delete selHP;
03089 theResult->erase(ih);
03090 #ifdef debug
03091 found=true;
03092 #endif
03093 break;
03094 }
03095 }
03096 #ifdef debug
03097 if(!found) G4cout<<"*G4QFragmentation::Breeder:nH="<<theResult->size()<<G4endl;
03098 #endif
03099
03100 theHadrons=curString->FragmentString(true);
03101 #ifdef debug
03102 G4cout<<"G4QFrag::Breeder: tH="<<theHadrons<<",nH="<<theResult->size()<<G4endl;
03103 #endif
03104 }
03105 #ifdef debug
03106 G4cout<<"*G4QFragmentation::Breeder: *CanTryToDecay?* nH="<<theHadrons<<", next="
03107 <<next<<" =? nS="<<strings.size()<<", nR="<<theResult->size()<<G4endl;
03108 #endif
03109 if(!theHadrons && next == strings.size() && !(theResult->size()))
03110 {
03111 G4QContent miQC=curString->GetQC();
03112 G4int miPDG=miQC.GetSPDGCode();
03113 if(miPDG == 10)
03114 {
03115 G4QChipolino QCh(miQC);
03116 G4QPDGCode h1QPDG=QCh.GetQPDG1();
03117 G4double h1M =h1QPDG.GetMass();
03118 G4QPDGCode h2QPDG=QCh.GetQPDG2();
03119 G4double h2M =h2QPDG.GetMass();
03120 G4double ttM =curString4M.m();
03121 if(h1M+h2M<ttM+eps)
03122 {
03123 G4LorentzVector h14M(0.,0.,0.,h1M);
03124 G4LorentzVector h24M(0.,0.,0.,h2M);
03125 if(std::fabs(ttM-h1M-h2M)<=eps)
03126 {
03127 G4double part1=h1M/(h1M+h2M);
03128 h14M=part1*curString4M;
03129 h24M=curString4M-h14M;
03130 }
03131 else
03132 {
03133 if(!G4QHadron(curString4M).DecayIn2(h14M,h24M))
03134 {
03135 G4cerr<<"***G4QFragmentation::Breeder: tM="<<ttM<<"->h1="<<h1QPDG<<"("
03136 <<h1M<<")+h2="<<h1QPDG<<"("<<h2M<<")="<<h1M+h2M<<G4endl;
03137 G4Exception("G4QFragmentation::Breeder:","72",FatalException,"ChiDec");
03138 }
03139 }
03140 G4QHadron* h1H = new G4QHadron(h1QPDG.GetPDGCode(),h14M);
03141 theResult->push_back(h1H);
03142 #ifdef debug
03143 G4LorentzVector f4M=h1H->Get4Momentum();
03144 G4int fPD=h1H->GetPDGCode();
03145 G4int fCg=h1H->GetCharge();
03146 G4int fBN=h1H->GetBaryonNumber();
03147 G4cout<<"-EMC->>G4QFragment::Breeder: String=Hadr ChiPro1 is filled, f4M="
03148 <<f4M<<", fPDG="<<fPD<<", fCg="<<fCg<<", fBN="<<fBN<<G4endl;
03149 #endif
03150 G4QHadron* h2H = new G4QHadron(h2QPDG.GetPDGCode(),h24M);
03151 theResult->push_back(h2H);
03152 #ifdef debug
03153 G4LorentzVector s4M=h2H->Get4Momentum();
03154 G4int sPD=h2H->GetPDGCode();
03155 G4int sCg=h2H->GetCharge();
03156 G4int sBN=h2H->GetBaryonNumber();
03157 G4cout<<"-EMC->>G4QFragmentation::Breeder: String=Hadr ChiPro2 is filled, s4M="
03158 <<s4M<<", sPDG="<<sPD<<", sCg="<<sCg<<", sBN="<<sBN<<G4endl;
03159 #endif
03160 #ifdef edebug
03161 G4cout<<"-EMC-..Chi..G4QFragmentation::Breeder: DecayCHECK, Ch4M="
03162 <<curString4M<<", d4M="<<curString4M-h14M-h24M<<G4endl;
03163 #endif
03164 break;
03165 }
03166 else
03167 {
03168 G4Quasmon* stringQuasmon = new G4Quasmon(miQC, curString4M);
03169 theQuasmons.push_back(stringQuasmon);
03170 break;
03171 }
03172 }
03173 else if(miPDG)
03174 {
03175 if (miPDG>0 && miPDG%10 < 3) miPDG+=2;
03176 else if(miPDG<0 && (-miPDG)%10< 3) miPDG-=2;
03177 G4Quasmon Quasm;
03178 G4QHadron* sHad = new G4QHadron(miPDG,curString4M);
03179 G4QHadronVector* tmpQHadVec=Quasm.DecayQHadron(sHad);
03180 G4int tmpN=tmpQHadVec->size();
03181 #ifdef debug
03182 G4cout<<"G4QFragmentation::Breeder: Decay the Last, Res#H="<<tmpN<<G4endl;
03183 #endif
03184 if(tmpN>1)
03185 {
03186 for(G4int aH=0; aH < tmpN; aH++)
03187 {
03188 theResult->push_back((*tmpQHadVec)[aH]);
03189 #ifdef debug
03190 G4QHadron* prodH =(*tmpQHadVec)[aH];
03191 G4LorentzVector p4M=prodH->Get4Momentum();
03192 G4int PDG=prodH->GetPDGCode();
03193 G4int Chg=prodH->GetCharge();
03194 G4int BaN=prodH->GetBaryonNumber();
03195 G4cout<<"-EMC->>G4QFragment::Breeder:String=Hadr,H#"<<aH<<" filled, 4M="
03196 <<p4M<<", PDG="<<PDG<<", Chg="<<Chg<<", BaN="<<BaN<<G4endl;
03197 #endif
03198 }
03199 }
03200 else
03201 {
03202 G4Quasmon* stringQuasmon = new G4Quasmon(miQC, curString4M);
03203 #ifdef debug
03204 G4cout<<"G4QFragmentat::Breeder:==> to Quasm="<<miQC<<curString4M<<", Nuc="
03205 <<theNucleus<<theNucleus.Get4Momentum()<<", NString="<<strings.size()
03206 <<", nR="<<theResult->size()<<", nQ="<<theQuasmons.size()<<G4endl;
03207 #endif
03208 theQuasmons.push_back(stringQuasmon);
03209 delete sHad;
03210 tmpQHadVec->clear();
03211 delete tmpQHadVec;
03212 break;
03213 }
03214 tmpQHadVec->clear();
03215 delete tmpQHadVec;
03216 break;
03217 }
03218 }
03219 }
03220 }
03221
03222 #ifdef debug
03223 G4cout<<"G4QFragmentation::Breeder: theH="<<theHadrons<<"?=0, next="<<next<<G4endl;
03224 #endif
03225 if(!theHadrons && next < strings.size())
03226 {
03227
03228 G4QContent miQC=curString->GetQC();
03229 G4int miPDG=miQC.GetSPDGCode();
03230 #ifdef debug
03231 G4cout<<"---->>G4QFragmentation::Breeder: SQC="<<miQC<<", miSPDG="<<miPDG<<G4endl;
03232 #endif
03233 G4double miM=0.;
03234 if(miPDG!=10) miM=G4QPDGCode(miPDG).GetMass();
03235 else
03236 {
03237 G4QChipolino QCh(miQC);
03238 miM=QCh.GetQPDG1().GetMass()+QCh.GetQPDG2().GetMass();
03239 }
03240 G4double cM2=curString4M.m2();
03241 #ifdef debug
03242 G4cout<<"---->>G4QFragmentation::Breeder: minMass="<<miM<<", realM2="<<cM2<<G4endl;
03243 #endif
03244 G4double cM=0.;
03245 if(cM2>0.)
03246 {
03247 cM=std::sqrt(cM2);
03248 if(std::fabs(cM-miM) < eps)
03249 {
03250 if(miPDG!=10)
03251 {
03252 G4QHadron* sHad = new G4QHadron(miPDG,curString4M);
03253 theResult->push_back(sHad);
03254 #ifdef debug
03255 G4cout<<"----->>G4QFragmentation::Breeder:S->H="<<miPDG<<curString4M<<G4endl;
03256 #endif
03257 }
03258 else
03259 {
03260 G4QChipolino QCh(miQC);
03261 G4QPDGCode h1QPDG=QCh.GetQPDG1();
03262 G4double h1M =h1QPDG.GetMass();
03263 G4QPDGCode h2QPDG=QCh.GetQPDG2();
03264 G4double h2M =h2QPDG.GetMass();
03265 G4double pt1 =h1M/(h1M+h2M);
03266 G4LorentzVector h14M=pt1*curString4M;
03267 G4LorentzVector h24M=curString4M-h14M;
03268 G4QHadron* h1H = new G4QHadron(h1QPDG.GetPDGCode(),h14M);
03269 theResult->push_back(h1H);
03270 #ifdef debug
03271 G4LorentzVector f4M=h1H->Get4Momentum();
03272 G4int fPD=h1H->GetPDGCode();
03273 G4int fCg=h1H->GetCharge();
03274 G4int fBN=h1H->GetBaryonNumber();
03275 G4cout<<"-EMC->>G4QFragmentation::Breeder:Str=2HadrAR Prod-F is filled, f4M="
03276 <<f4M<<", fPDG="<<fPD<<", fCg="<<fCg<<", fBN="<<fBN<<G4endl;
03277 #endif
03278 G4QHadron* h2H = new G4QHadron(h2QPDG.GetPDGCode(),h24M);
03279 theResult->push_back(h2H);
03280 #ifdef debug
03281 G4LorentzVector s4M=h2H->Get4Momentum();
03282 G4int sPD=h2H->GetPDGCode();
03283 G4int sCg=h2H->GetCharge();
03284 G4int sBN=h2H->GetBaryonNumber();
03285 G4cout<<"-EMC->>G4QFragmentation::Breeder:Str=2HadrAR Prod-S is filled, s4M="
03286 <<s4M<<", sPDG="<<sPD<<", sCg="<<sCg<<", sBN="<<sBN<<G4endl;
03287 #endif
03288 }
03289 continue;
03290 }
03291 else
03292 {
03293 G4ThreeVector cP=curString4M.vect();
03294 G4double cE=curString4M.e();
03295 G4ThreeVector curV=cP/cE;
03296 G4double miM2=miM*miM;
03297 G4int restr=0;
03298 G4int fustr=0;
03299 G4double selX=0.;
03300 G4double maD=-DBL_MAX;
03301 G4double Vmin=DBL_MAX;
03302 G4LorentzVector s4M(0.,0.,0.,0.);
03303 #ifdef debug
03304 G4cout<<"G4QFr::Breed:TryRecover,V="<<curV<<",cM2="<<cM2<<",miM="<<miM<<G4endl;
03305 #endif
03306 nOfStr=strings.size();
03307 for(restr=next; restr < nOfStr; ++restr) if(restr != astring)
03308 {
03309 #ifdef debug
03310 G4cout<<"G4QFragmentation::Breeder: rS="<<restr<<", nS="<<nOfStr<<G4endl;
03311 #endif
03312 G4QString* pString=strings[restr];
03313 #ifdef debug
03314 G4cout<<"G4QFragmentation::Breeder: pString="<<pString<<G4endl;
03315 #endif
03316 G4LorentzVector p4M=pString->Get4Momentum();
03317 #ifdef debug
03318 G4cout<<"G4QFragmentation::Breeder: p4M="<<p4M<<G4endl;
03319 #endif
03320 G4ThreeVector pP=p4M.vect();
03321 G4double pE=p4M.e();
03322 G4double D2=cE*pE-cP.dot(pP);
03323 G4double pM2=p4M.m2();
03324 #ifdef debug
03325 G4cout<<"G4QFrag::Breeder: pM2="<<pM2<<",miM2="<<miM2<<",cM2="<<cM2<<G4endl;
03326 #endif
03327 G4double dM4=pM2*(miM2-cM2);
03328 G4double D=D2*D2+dM4;
03329 #ifdef debug
03330 G4cout<<"G4QFragmentation::Breeder: D="<<D<<",dM4="<<dM4<<",D2="<<D2<<G4endl;
03331 #endif
03332 G4double x=-1.;
03333 if(D > 0. && pM2>.01) x=(std::sqrt(D)-D2)/pM2;
03334 #ifdef debug
03335 else G4cout<<"G4QFragment::Breeder: D="<<D<<",D2="<<D2<<",pM4="<<dM4<<G4endl;
03336 G4cout<<"G4QFragmentation::Breeder: pM2="<<pM2<<",D2="<<D2<<",x="<<x<<G4endl;
03337 #endif
03338 if(x > 0. && x < 1.)
03339 {
03340 G4QContent pQC=pString->GetQC();
03341 G4int pPDG=pQC.GetSPDGCode();
03342 G4double pM=0.;
03343 if(pPDG==10)
03344 {
03345 G4QChipolino QCh(pQC);
03346 pM=QCh.GetQPDG1().GetMass()+QCh.GetQPDG2().GetMass();
03347 }
03348 else pM=G4QPDGCode(pPDG).GetMass();
03349 G4double rM=std::sqrt(pM2);
03350 G4double delta=(1.-x)*rM-pM;
03351 if(delta > 0. && delta > maD)
03352 {
03353 maD=delta;
03354 #ifdef debug
03355 G4cout<<"G4QFragmentation::Breeder: Subtr,S#"<<restr<<",d="<<maD<<G4endl;
03356 #endif
03357 fustr=restr;
03358 selX=x;
03359 s4M=p4M;
03360 }
03361 }
03362 else if(x <= 0.)
03363 {
03364 G4ThreeVector pV=pP/pE;
03365 G4double dV=(curV-pV).mag2();
03366 if(dV < Vmin)
03367 {
03368 #ifdef debug
03369 G4cout<<"G4QFragmentation::Breeder: FreeAdd,S#"<<restr<<",x="<<x<<G4endl;
03370 #endif
03371 Vmin=dV;
03372 fustr=restr;
03373 selX=x;
03374 s4M=p4M;
03375 }
03376 }
03377 #ifdef debug
03378 G4cout<<"G4QFragmentation::Breeder:EndOfLOOP r="<<restr<<"<"<<nOfStr<<G4endl;
03379 #endif
03380 }
03381 #ifdef debug
03382 G4cout<<"G4QFragmentation::Breeder: AfterLOOP fustr="<<fustr<<G4endl;
03383 #endif
03384 if(fustr)
03385 {
03386 #ifdef edebug
03387 G4LorentzVector sum4M=s4M+curString4M;
03388 G4cout<<"G4QFragmentation::Breeder: Found Sum4M="<<sum4M<<G4endl;
03389 #endif
03390 G4QString* pString=strings[fustr];
03391 curString4M+=selX*s4M;
03392 if(std::abs(miPDG)%10 > 2)
03393 {
03394 G4Quasmon Quasm;
03395 G4QHadron* sHad = new G4QHadron(miPDG,curString4M);
03396 G4QHadronVector* tmpQHadVec=Quasm.DecayQHadron(sHad);
03397 #ifdef debug
03398 G4cout<<"G4QFragmentation::Breeder:DecStH,nH="<<tmpQHadVec->size()<<G4endl;
03399 #endif
03400 for(unsigned aH=0; aH < tmpQHadVec->size(); aH++)
03401 {
03402 theResult->push_back((*tmpQHadVec)[aH]);
03403 #ifdef debug
03404 G4QHadron* prodH =(*tmpQHadVec)[aH];
03405 G4LorentzVector p4M=prodH->Get4Momentum();
03406 G4int PDG=prodH->GetPDGCode();
03407 G4int Chg=prodH->GetCharge();
03408 G4int BaN=prodH->GetBaryonNumber();
03409 G4cout<<"-EMC->>G4QFragmentation::Breeder:St=Had,pH#"<<aH<<" filled, 4M="
03410 <<p4M<<", PDG="<<PDG<<", Chg="<<Chg<<", BaN="<<BaN<<G4endl;
03411 #endif
03412 }
03413 tmpQHadVec->clear();
03414 delete tmpQHadVec;
03415 }
03416 else if(miPDG == 10)
03417 {
03418 G4QChipolino QCh(miQC);
03419 G4QPDGCode h1QPDG=QCh.GetQPDG1();
03420 G4double h1M =h1QPDG.GetMass();
03421 G4QPDGCode h2QPDG=QCh.GetQPDG2();
03422 G4double h2M =h2QPDG.GetMass();
03423 G4double ttM =curString4M.m();
03424 if(h1M+h2M<ttM+eps)
03425 {
03426 G4LorentzVector h14M(0.,0.,0.,h1M);
03427 G4LorentzVector h24M(0.,0.,0.,h2M);
03428 if(std::fabs(ttM-h1M-h2M)<=eps)
03429 {
03430 G4double part1=h1M/(h1M+h2M);
03431 h14M=part1*curString4M;
03432 h24M=curString4M-h14M;
03433 }
03434 else
03435 {
03436 if(!G4QHadron(curString4M).DecayIn2(h14M,h24M))
03437 {
03438 G4cerr<<"***G4QFragmentation::Breeder: tM="<<ttM<<"->h1="<<h1QPDG
03439 <<"("<<h1M<<")+h2="<<h1QPDG<<"("<<h2M<<")="<<h1M+h2M<<G4endl;
03440 G4Exception("G4QFragmentation::Breeder:","72",FatalException,"ChDe");
03441 }
03442 }
03443 G4QHadron* h1H = new G4QHadron(h1QPDG.GetPDGCode(),h14M);
03444 theResult->push_back(h1H);
03445 #ifdef debug
03446 G4LorentzVector f4M=h1H->Get4Momentum();
03447 G4int fPD=h1H->GetPDGCode();
03448 G4int fCg=h1H->GetCharge();
03449 G4int fBN=h1H->GetBaryonNumber();
03450 G4cout<<"-EMC->>G4QFragmentation::Breeder:Str=Hadr Prod-F's filled, f4M="
03451 <<f4M<<", fPDG="<<fPD<<", fCg="<<fCg<<", fBN="<<fBN<<G4endl;
03452 #endif
03453 G4QHadron* h2H = new G4QHadron(h2QPDG.GetPDGCode(),h24M);
03454 theResult->push_back(h2H);
03455 #ifdef debug
03456 G4LorentzVector s4M=h2H->Get4Momentum();
03457 G4int sPD=h2H->GetPDGCode();
03458 G4int sCg=h2H->GetCharge();
03459 G4int sBN=h2H->GetBaryonNumber();
03460 G4cout<<"-EMC->>G4QFragmentation::Breeder:Str=Hadr Prod-S's filled, s4M="
03461 <<s4M<<", sPDG="<<sPD<<", sCg="<<sCg<<", sBN="<<sBN<<G4endl;
03462 #endif
03463 #ifdef edebug
03464 G4cout<<"-EMC-Chipo.G4QFragmentation::Breeder:DecCHECK,c4M="<<curString4M
03465 <<", ChQC="<<miQC<<", d4M="<<curString4M-h14M-h24M<<G4endl;
03466 #endif
03467 }
03468 else
03469 {
03470 G4cerr<<"***G4QFragm::Breeder:tM="<<ttM<<miQC<<"->h1="<<h1QPDG<<"(" <<h1M
03471 <<")+h2="<<h1QPDG<<"("<<h2M<<") = "<<h1M+h2M<<G4endl;
03472 G4Exception("G4QFragmentation::Breeder:","72",FatalException,"ChiDecMa");
03473 }
03474 }
03475 else
03476 {
03477 G4QHadron* sHad = new G4QHadron(miPDG,curString4M);
03478 theResult->push_back(sHad);
03479 #ifdef debug
03480 G4cout<<"-EMC->>G4QFragmentation::Breeder:Str=Hadr Filled, 4M="
03481 <<curString4M<<", PDG="<<miPDG<<G4endl;
03482 #endif
03483 }
03484 G4double corF=1-selX;
03485 G4QParton* Left=pString->GetLeftParton();
03486 G4QParton* Right=pString->GetRightParton();
03487 Left->Set4Momentum(corF*Left->Get4Momentum());
03488 Right->Set4Momentum(corF*Right->Get4Momentum());
03489 #ifdef edebug
03490 G4cout<<"-EMC-...Cor...G4QFragmentation::Breeder:CorCHECK Sum="<<sum4M
03491 <<" =? "<<curString4M+pString->Get4Momentum()<<", M="<<miM<<" ?= "
03492 <<curString4M.m()<<G4endl;
03493 #endif
03494 #ifdef debug
03495 G4cout<<"---->>G4QFragmentation::Breeder:*Corrected* String->Hadr="<<miPDG
03496 <<curString4M<<" by String #"<<fustr<<G4endl;
03497 #endif
03498 continue;
03499 }
03500 }
03501 }
03502 }
03503 #ifdef debug
03504 else G4cout<<"***G4QFragmentation::Breeder: **No SSCorrection**,next="<<next<<G4endl;
03505 #endif
03506
03507 G4QParton* lpcS=curString->GetLeftParton();
03508 G4QParton* rpcS=curString->GetRightParton();
03509 G4int lPDGcS=lpcS->GetPDGCode();
03510 G4int rPDGcS=rpcS->GetPDGCode();
03511 if (lPDGcS==3 && rPDGcS==-3)
03512 {
03513 lpcS->SetPDGCode( 1);
03514 rpcS->SetPDGCode(-1);
03515 }
03516 else if(lPDGcS==-3 && rPDGcS==3)
03517 {
03518 lpcS->SetPDGCode(-1);
03519 rpcS->SetPDGCode( 1);
03520 }
03521
03522 G4int nofRH=theResult->size();
03523 #ifdef debug
03524 G4cout<<"G4QFragmentation::Breeder: theH="<<theHadrons<<", #OfH="<<nofRH<<G4endl;
03525 #endif
03526 if(!theHadrons && nofRH)
03527 {
03528 #ifdef debug
03529 G4cout<<"!G4QFragmentation::Breeder:Can try SHCor, nH="<<theResult->size()<<G4endl;
03530 #endif
03531
03532 G4QContent miQC=curString->GetQC();
03533 G4int miPDG=miQC.GetSPDGCode();
03534 G4double miM=0.;
03535 if(miPDG==10)
03536 {
03537 G4QChipolino QCh(miQC);
03538 miM=QCh.GetQPDG1().GetMass()+QCh.GetQPDG2().GetMass();
03539 }
03540 else miM=G4QPDGCode(miPDG).GetMass();
03541 G4double spM=0.;
03542 G4ThreeVector cP=curString4M.vect();
03543 G4double cE=curString4M.e();
03544 G4ThreeVector curV=cP/cE;
03545 G4int reha=0;
03546 G4int fuha=-1;
03547 G4double dMmin=DBL_MAX;
03548 G4LorentzVector s4M(0.,0.,0.,0.);
03549 G4double sM=0.;
03550 for (reha=0; reha < nofRH; reha++)
03551 {
03552 G4QHadron* pHadron=(*theResult)[reha];
03553 G4LorentzVector p4M=pHadron->Get4Momentum();
03554 G4double pM=p4M.m();
03555 G4LorentzVector t4M=p4M+curString4M;
03556 G4double tM2=t4M.m2();
03557 if(tM2 >= sqr(pM+miM+eps))
03558 {
03559 G4double tM=std::sqrt(tM2);
03560 G4double dM=tM-pM-miM;
03561 if(dM < dMmin)
03562 {
03563 dMmin=dM;
03564 fuha=reha;
03565 spM=pM;
03566 s4M=t4M;
03567 sM=tM;
03568 }
03569 }
03570 #ifdef debug
03571 else G4cout<<"G4QFragmentation::Breeder:H# "<<reha<<",tM="<<std::sqrt(tM2)<<" < "
03572 <<" mS="<<miM<<" + mH="<<pM<<" = "<<pM+miM<<G4endl;
03573 #endif
03574 }
03575 #ifdef debug
03576 G4cout<<"G4QFragment::Breeder: fuha="<<fuha<<", dMmin="<<dMmin<<G4endl;
03577 #endif
03578 if(fuha>-1)
03579 {
03580 G4QHadron* pHadron=(*theResult)[fuha];
03581 G4LorentzVector mi4M(0.,0.,0.,miM);
03582 if(miM+spM<sM+eps)
03583 {
03584 G4LorentzVector sp4M(0.,0.,0.,spM);
03585 if(std::fabs(sM-miM-spM)<=eps)
03586 {
03587 G4double part1=miM/(miM+spM);
03588 mi4M=part1*s4M;
03589 sp4M=s4M-mi4M;
03590 }
03591 else
03592 {
03593 if(!G4QHadron(s4M).DecayIn2(mi4M,sp4M))
03594 {
03595 G4cerr<<"***G4QFragmentation::Breeder: *SH*, tM="<<sM<<"->h1=("<<miPDG<<")"
03596 <<miM<<" + h2="<<spM<<" = "<<miM+spM<<G4endl;
03597 G4Exception("G4QFragmentation::Breeder:","72",FatalException,"SHChipoDec");
03598 }
03599 }
03600 pHadron->Set4Momentum(sp4M);
03601 #ifdef debug
03602 G4cout<<"-EMC->...>G4QFragmentation::Breeder: H# "<<fuha<<" is updated, new4M="
03603 <<sp4M<<G4endl;
03604 #endif
03605 }
03606 else
03607 {
03608 G4cerr<<"***G4QFragm::Breeder: HS Failed, tM="<<sM<<"->h1M=("<<miPDG<<")"<<miM
03609 <<"+h2M="<<spM<<" = "<<miM+spM<<G4endl;
03610 G4Exception("G4QFragmentation::Breeder:","72",FatalException,"HSChipoDecMass");
03611 }
03612 if(std::abs(miPDG)%10 > 2)
03613 {
03614 G4Quasmon Quasm;
03615 G4QHadron* sHad = new G4QHadron(miPDG,mi4M);
03616 G4QHadronVector* tmpQHadVec=Quasm.DecayQHadron(sHad);
03617 #ifdef debug
03618 G4cout<<"G4QFragment::Breeder:*HS* DecStrHad, nH="<<tmpQHadVec->size()<<G4endl;
03619 #endif
03620 for(unsigned aH=0; aH < tmpQHadVec->size(); aH++)
03621 {
03622 theResult->push_back((*tmpQHadVec)[aH]);
03623 #ifdef debug
03624 G4QHadron* prodH =(*tmpQHadVec)[aH];
03625 G4LorentzVector p4M=prodH->Get4Momentum();
03626 G4int PDG=prodH->GetPDGCode();
03627 G4int Chg=prodH->GetCharge();
03628 G4int BaN=prodH->GetBaryonNumber();
03629 G4cout<<"-EMC->>G4QFragmentation::Breeder: Str+Hadr PrH#"<<aH<<" filled, 4M="
03630 <<p4M<<", PDG="<<PDG<<", Chg="<<Chg<<", BaN="<<BaN<<G4endl;
03631 #endif
03632 }
03633 tmpQHadVec->clear();
03634 delete tmpQHadVec;
03635 }
03636 else if(miPDG == 10)
03637 {
03638 G4QChipolino QCh(miQC);
03639 G4QPDGCode h1QPDG=QCh.GetQPDG1();
03640 G4double h1M =h1QPDG.GetMass();
03641 G4QPDGCode h2QPDG=QCh.GetQPDG2();
03642 G4double h2M =h2QPDG.GetMass();
03643 G4double ttM =curString4M.m();
03644 if(h1M+h2M<miM+eps)
03645 {
03646 G4LorentzVector h14M(0.,0.,0.,h1M);
03647 G4LorentzVector h24M(0.,0.,0.,h2M);
03648 if(std::fabs(ttM-h1M-h2M)<=eps)
03649 {
03650 G4double part1=h1M/(h1M+h2M);
03651 h14M=part1*mi4M;
03652 h24M=mi4M-h14M;
03653 }
03654 else
03655 {
03656 if(!G4QHadron(mi4M).DecayIn2(h14M,h24M))
03657 {
03658 G4cerr<<"***G4QFragmentation::Breeder: HS tM="<<ttM<<"->h1="<<h1QPDG<<"("
03659 <<h1M<<")+h2="<<h1QPDG<<"("<<h2M<<")="<<h1M+h2M<<G4endl;
03660 G4Exception("G4QFragmentation::Breeder:","72",FatalException,"ChipoDec");
03661 }
03662 }
03663 G4QHadron* h1H = new G4QHadron(h1QPDG.GetPDGCode(),h14M);
03664 theResult->push_back(h1H);
03665 #ifdef debug
03666 G4LorentzVector f4M=h1H->Get4Momentum();
03667 G4int fPD=h1H->GetPDGCode();
03668 G4int fCg=h1H->GetCharge();
03669 G4int fBN=h1H->GetBaryonNumber();
03670 G4cout<<"-EMC->>G4QFragmentation::Breeder: CorStrHadr Prod-1 is filled, f4M="
03671 <<f4M<<", fPDG="<<fPD<<", fCg="<<fCg<<", fBN="<<fBN<<G4endl;
03672 #endif
03673 G4QHadron* h2H = new G4QHadron(h2QPDG.GetPDGCode(),h24M);
03674 theResult->push_back(h2H);
03675 #ifdef debug
03676 G4LorentzVector n4M=h2H->Get4Momentum();
03677 G4int nPD=h2H->GetPDGCode();
03678 G4int nCg=h2H->GetCharge();
03679 G4int nBN=h2H->GetBaryonNumber();
03680 G4cout<<"-EMC->>G4QFragmentation::Breeder: CorStrHadr Prod-2 is filled, n4M="
03681 <<n4M<<", nPDG="<<nPD<<", nCg="<<nCg<<", nBN="<<nBN<<G4endl;
03682 #endif
03683 #ifdef edebug
03684 G4cout<<"-EMC-...HS-Chipo...G4QFragmentation::Breeder:DecCHECK, Ch4M="<<mi4M
03685 <<", ChQC="<<miQC<<", d4M="<<mi4M-h14M-h24M<<G4endl;
03686 #endif
03687 }
03688 }
03689 else
03690 {
03691 G4QHadron* sHad = new G4QHadron(miPDG, mi4M);
03692 theResult->push_back(sHad);
03693 #ifdef debug
03694 G4cout<<"----->>G4QFragmentation::Breeder: CorStr=Hadr is Filled, 4M="
03695 <<curString4M<<", StPDG="<<miPDG<<G4endl;
03696 #endif
03697 }
03698 #ifdef edebug
03699 G4cout<<"-EMC-...Cor...G4QFragmentation::Breeder:StHadCor CHECK Sum="<<s4M
03700 <<" =? "<<mi4M+pHadron->Get4Momentum()<<G4endl;
03701 #endif
03702 #ifdef debug
03703 G4cout<<"------->>G4QFragmentation::Breeder: *Corrected* String+Hadr="<<miPDG
03704 <<mi4M<<" by Hadron #"<<reha<<G4endl;
03705 #endif
03706 continue;
03707 }
03708 else
03709 {
03710 #ifdef debug
03711 G4cout<<"G4QFragmentation::Breeder: Str+Hadr Failed, 4M="<<curString4M
03712 <<", PDG="<<miPDG<<" -> Now try to recover the string as a hadron"<<G4endl;
03713 #endif
03714
03715
03716
03717
03718
03719
03720
03721
03722
03723
03724
03725
03726
03727
03728
03729
03730
03731
03732
03733
03734
03735 }
03736
03737 G4QContent curStringQC=curString->GetQC();
03738 G4Quasmon* stringQuasmon = new G4Quasmon(curStringQC, curString4M);
03739 theQuasmons.push_back(stringQuasmon);
03740 continue;
03741 }
03742 }
03743 G4Quasmon tmpQ;
03744 G4int nHfin=0;
03745 if(theHadrons) nHfin=theHadrons->size();
03746 else
03747 {
03748 G4LorentzVector ss4M(0.,0.,0.,0.);
03749 G4QContent ssQC(0,0,0,0,0,0);
03750 G4int tnSt=strings.size();
03751 for(G4int i=astring; i < tnSt; ++i)
03752 {
03753 G4LorentzVector pS4M=strings[i]->Get4Momentum();
03754 ss4M+=pS4M;
03755 G4QContent pSQC=strings[i]->GetQC();
03756 ssQC+=pSQC;
03757 #ifdef debug
03758 G4cout<<"=--=>G4QFragmentation::Breeder:S#"<<i<<",4M="<<pS4M<<",QC="<<pSQC<<G4endl;
03759 #endif
03760 }
03761 #ifdef debug
03762 G4cout<<"==>G4QFragmentation::Breeder:AllStrings are summed up in a Quasmon"<<G4endl;
03763 #endif
03764 G4Quasmon* stringQuasmon = new G4Quasmon(ssQC, ss4M);
03765 theQuasmons.push_back(stringQuasmon);
03766 break;
03767 }
03768 #ifdef debug
03769 G4cout<<"G4QFragmentation::Breeder: Trying to decay hadrons #ofHRes="<<nHfin<<G4endl;
03770 #endif
03771 for(G4int aTrack=0; aTrack<nHfin; aTrack++)
03772 {
03773 G4QHadron* curHadron=(*theHadrons)[aTrack];
03774 G4int hPDG=curHadron->GetPDGCode();
03775 G4LorentzVector curH4M=curHadron->Get4Momentum();
03776 G4int curHCh=curHadron->GetCharge();
03777 G4int curHBN=curHadron->GetBaryonNumber();
03778 #ifdef debug
03779 G4cout<<"----->>G4QFragmentation::Breeder:S#"<<astring<<",H#"<<aTrack<<",PDG="<<hPDG
03780 <<",4M="<<curHadron->Get4Momentum()<<G4endl;
03781 #endif
03782 if(std::abs(hPDG)%10 > 2)
03783 {
03784 G4QHadronVector* tmpQHadVec=tmpQ.DecayQHadron(curHadron);
03785 #ifdef debug
03786 G4cout<<"G4QFragmentation::Breeder:-DECAY'S DONE-,nH="<<tmpQHadVec->size()<<G4endl;
03787 #endif
03788 for(unsigned aH=0; aH < tmpQHadVec->size(); aH++)
03789 {
03790 theResult->push_back((*tmpQHadVec)[aH]);
03791
03792 G4QHadron* prodH =(*tmpQHadVec)[aH];
03793 G4LorentzVector p4M=prodH->Get4Momentum();
03794 G4int Chg=prodH->GetCharge();
03795 G4int BaN=prodH->GetBaryonNumber();
03796 curString4M-=p4M;
03797 curStrChg-=Chg;
03798 curStrBaN-=BaN;
03799 curH4M-=p4M;
03800 curHCh-=Chg;
03801 curHBN-=BaN;
03802 #ifdef edebug
03803 G4int PDG=prodH->GetPDGCode();
03804 G4cout<<"-EMC->>G4QFragmentation::Breeder:String*Filled, 4M="<<p4M<<", PDG="<<PDG
03805 <<", Chg="<<Chg<<", BaN="<<BaN<<G4endl;
03806 #endif
03807 }
03808 #ifdef edebug
03809 G4cout<<"-EMC-.G4QFr::Br:Dec,r4M="<<curH4M<<",rC="<<curHCh<<",rB="<<curHBN<<G4endl;
03810 #endif
03811 tmpQHadVec->clear();
03812 delete tmpQHadVec;
03813 }
03814 else
03815 {
03816 theResult->push_back(curHadron);
03817
03818 curString4M-=curH4M;
03819 G4int curCh=curHadron->GetCharge();
03820 G4int curBN=curHadron->GetBaryonNumber();
03821 curStrChg-=curCh;
03822 curStrBaN-=curBN;
03823 #ifdef edebug
03824 G4cout<<"-EMC->>-->>G4QFragmentation::Breeder: curH filled 4M="<<curH4M<<",PDG="
03825 <<curHadron->GetPDGCode()<<", Chg="<<curCh<<", BaN="<<curBN<<G4endl;
03826 #endif
03827 }
03828 }
03829
03830 if(theHadrons) delete theHadrons;
03831 #ifdef edebug
03832 G4cout<<"-EMC-.........G4QFragmentation::Breeder: StringDecay CHECK, r4M="<<curString4M
03833 <<", rChg="<<curStrChg<<", rBaN="<<curStrBaN<<G4endl;
03834 #endif
03835
03836 if(curStrChg || curStrBaN || curString4M.t() > eps || std::fabs(curString4M.x()) > eps
03837 || std::fabs(curString4M.y()) > eps || std::fabs(curString4M.z()) > eps )
03838 {
03839 G4double dEn=curString4M.t();
03840 G4double dPx=curString4M.x();
03841 G4double dPy=curString4M.y();
03842 G4double dPz=curString4M.z();
03843 G4int nHadr=theResult->size();
03844 G4double hEn=0.;
03845 G4double hPx=0.;
03846 G4double hPy=0.;
03847 G4double hPz=0.;
03848 G4int hCh=0;
03849 G4int hBN=0;
03850 G4double mEn=0.;
03851 G4double mPx=0.;
03852 G4double mPy=0.;
03853 G4double mPz=0.;
03854 G4int mCh=0;
03855 G4int mBN=0;
03856 for(G4int i=0; i<nHadr; i++)
03857 {
03858 mEn=hEn;
03859 mPx=hPx;
03860 mPy=hPy;
03861 mPz=hPz;
03862 mCh=hCh;
03863 mBN=hBN;
03864 G4QHadron* curHadr = (*theResult)[i];
03865 G4LorentzVector hI4M = curHadr->Get4Momentum();
03866 hEn=hI4M.t();
03867 hPx=hI4M.x();
03868 hPy=hI4M.y();
03869 hPz=hI4M.z();
03870 hCh=curHadr->GetCharge();
03871 hBN=curHadr->GetBaryonNumber();
03872 G4cout<<"G4QFragmentation::Breeder: H#"<<i<<", d4M="<<curString4M+hI4M
03873 <<",dCh="<<hCh+curStrChg<<",dBN="<<hBN+curStrBaN<<G4endl;
03874 if( !(hCh+curStrChg) && !(hBN+curStrBaN) && std::fabs(dEn+hEn)<eps &&
03875 std::fabs(dPx+hPx)<eps && std::fabs(dPy+hPy)<eps && std::fabs(dPz+hPz)<eps )
03876 {
03877 G4cout<<"G4QFragmentation::Breeder: ***Cured*** Redundent Hadron # "<<i<<G4endl;
03878 G4QHadron* theLast = (*theResult)[nHadr-1];
03879 curHadr->Set4Momentum(theLast->Get4Momentum());
03880 G4QPDGCode lQP=theLast->GetQPDG();
03881 if(lQP.GetPDGCode()!=10) curHadr->SetQPDG(lQP);
03882 else curHadr->SetQC(theLast->GetQC());
03883 theResult->pop_back();
03884 delete theLast;
03885 break;
03886 }
03887 if( !(hCh+mCh+curStrChg) && !(hBN+mBN+curStrBaN) && std::fabs(dEn+hEn+mEn)<eps &&
03888 std::fabs(dPx+hPx+mPx)<eps && std::fabs(dPy+hPy+mPy)<eps &&
03889 std::fabs(dPz+hPz+mPz)<eps && i>0)
03890 {
03891 G4cout<<"G4QFragmentation::Breeder:***Cured*** Redundent 2Hadrons i="<<i<<G4endl;
03892 G4QHadron* preHadr = (*theResult)[i-1];
03893 G4QHadron* theLast = (*theResult)[nHadr-1];
03894 if(i < nHadr-1)
03895 {
03896 preHadr->Set4Momentum(theLast->Get4Momentum());
03897 G4QPDGCode lQP=theLast->GetQPDG();
03898 if(lQP.GetPDGCode()!=10) preHadr->SetQPDG(lQP);
03899 else preHadr->SetQC(theLast->GetQC());
03900 }
03901 theResult->pop_back();
03902 delete theLast;
03903 theLast = (*theResult)[nHadr-2];
03904 if(i < nHadr-2)
03905 {
03906 curHadr->Set4Momentum(theLast->Get4Momentum());
03907 G4QPDGCode lQP=theLast->GetQPDG();
03908 if(lQP.GetPDGCode()!=10) curHadr->SetQPDG(lQP);
03909 else curHadr->SetQC(theLast->GetQC());
03910 }
03911 theResult->pop_back();
03912 delete theLast;
03913 nHadr=theResult->size();
03914 break;
03915 }
03916
03917 G4cout<<"*Warning*G4QFragmentation::Breeder: Nonconservation isn't cured!"<<G4endl;
03918 }
03919 }
03920
03921 }
03922 G4LorentzVector r4M=theNucleus.Get4Momentum();
03923 G4int rPDG=theNucleus.GetPDG();
03924 G4QHadron* resNuc = new G4QHadron(rPDG,r4M);
03925 theResult->push_back(resNuc);
03926 #ifdef edebug
03927 G4LorentzVector s4M(0.,0.,0.,0.);
03928 G4int rCh=totChg;
03929 G4int rBN=totBaN;
03930 G4int nHadr=theResult->size();
03931 G4int nQuasm=theQuasmons.size();
03932 G4cout<<"-EMCLS-G4QFragmentation::Breeder:#ofHadr="<<nHadr<<", #OfQuasm="<<nQuasm<<",rN="
03933 <<r4M.m()<<"="<<G4QNucleus(rPDG).GetGSMass()<<G4endl;
03934 for(G4int i=0; i<nHadr; i++)
03935 {
03936 G4LorentzVector hI4M=(*theResult)[i]->Get4Momentum();
03937 s4M+=hI4M;
03938 G4int hChg=(*theResult)[i]->GetCharge();
03939 rCh-=hChg;
03940 G4int hBaN=(*theResult)[i]->GetBaryonNumber();
03941 rBN-=hBaN;
03942 G4cout<<"-EMCLS-G4QFragmentation::Breeder:(1) Hadron#"<<i<<", 4M="<<hI4M<<", PDG="
03943 <<(*theResult)[i]->GetPDGCode()<<", C="<<hChg<<", B="<<hBaN<<G4endl;
03944 }
03945 for(G4int i=0; i<nQuasm; i++)
03946 {
03947 G4LorentzVector hI4M=theQuasmons[i]->Get4Momentum();
03948 s4M+=hI4M;
03949 G4int hChg=theQuasmons[i]->GetCharge();
03950 rCh-=hChg;
03951 G4int hBaN=theQuasmons[i]->GetBaryonNumber();
03952 rBN-=hBaN;
03953 G4cout<<"-EMCLS-G4QFragmentation::Breeder:(1) Quasmon#"<<i<<", 4M="<<hI4M<<", C="<<hChg
03954 <<", B="<<hBaN<<G4endl;
03955 }
03956 G4cout<<"-EMCLS-G4QFragm::Breed: LS r4M="<<s4M-totLS4M<<",rC="<<rCh<<",rB="<<rBN<<G4endl;
03957 #endif
03958
03959 G4int nRes=theResult->size();
03960 #ifdef ppdebug
03961 G4cout<<"G4QFragmentation::Breeder: Strings4M="<<ft4M<<", nRes="<<nRes<<G4endl;
03962 #endif
03963 G4ThreeVector LS3Mom=ft4M.v();
03964 G4ThreeVector LSDir=LS3Mom.unit();
03965 if(nRes > 2 && maxEn > 0.)
03966 {
03967 std::list<std::pair<G4double, G4QHadron*> > theSorted;
03968 std::list<std::pair<G4double, G4QHadron*> >::iterator current;
03969 for(G4int secondary = 0; secondary<nRes-1; ++secondary)
03970 {
03971 G4QHadron* ih =theResult->operator[](secondary);
03972 G4LorentzVector h4M =ih->Get4Momentum();
03973 G4double hM2 =ih->GetMass2();
03974 G4ThreeVector h3M =h4M.v();
03975 G4double toSort=DBL_MAX;
03976 if(hM2>0.00001) toSort=h4M.e()+h3M.dot(LSDir)/std::sqrt(hM2);
03977 #ifdef ppdebug
03978 G4cout<<"G4QFragmentation::Breeder:#"<<secondary<<",M2="<<hM2<<",s="<<toSort<<G4endl;
03979 #endif
03980 std::pair<G4double, G4QHadron*> it;
03981 it.first = toSort;
03982 it.second = ih;
03983 G4bool inserted = false;
03984 for(current = theSorted.begin(); current!=theSorted.end(); ++current)
03985 {
03986 if((*current).first > toSort)
03987 {
03988 theSorted.insert(current, it);
03989 inserted = true;
03990 break;
03991 }
03992 }
03993 if(!inserted) theSorted.push_back(it);
03994 }
03995 theResult->clear();
03996 G4LorentzVector q4M(0.,0.,0.,0.);
03997 G4QContent qQC(0,0,0,0,0,0);
03998 for(current = theSorted.begin(); current!=theSorted.end(); ++current)
03999 {
04000 G4QHadron* ih= (*current).second;
04001 G4LorentzVector h4M= ih->Get4Momentum();
04002 G4int hPDG= ih->GetPDGCode();
04003 G4double dE= 0.;
04004 G4bool tested=true;
04005 if (hPDG> 1111 && hPDG< 3335) dE=h4M.e()-ih->GetMass();
04006 else if(hPDG>-1111 && hPDG<1111 && hPDG!=22) dE=h4M.e();
04007
04008 else tested=false;
04009 #ifdef ppdebug
04010 G4cout<<"G4QFragmentation::Breeder:dE="<<dE<<",mE="<<maxEn<<",t="<<tested<<G4endl;
04011 #endif
04012 if(tested && dE < maxEn)
04013 {
04014 maxEn-=dE;
04015 q4M+=h4M;
04016 qQC+=ih->GetQC();
04017 #ifdef debug
04018 G4cout<<"%->G4QFragmentation::Breeder:Exclude,4M="<<h4M<<",dE="<<maxEn<<G4endl;
04019 #endif
04020 delete ih;
04021 }
04022 else theResult->push_back(ih);
04023 }
04024 G4Quasmon* softQuasmon = new G4Quasmon(qQC, q4M);
04025 #ifdef debug
04026 G4cout<<"%->G4QFragmentation::Breeder:QuasmonIsFilled,4M="<<q4M<<",QC="<<qQC<<G4endl;
04027 #endif
04028 if(q4M != vac4M) theQuasmons.push_back(softQuasmon);
04029 else delete softQuasmon;
04030 theResult->push_back(resNuc);
04031 #ifdef edebug
04032 G4LorentzVector f4M(0.,0.,0.,0.);
04033 G4int fCh=totChg;
04034 G4int fBN=totBaN;
04035 G4int nHd=theResult->size();
04036 G4int nQm=theQuasmons.size();
04037 G4cout<<"-EMCLS-G4QFragmentation::Breeder:#ofHadr="<<nHd<<", #OfQuasm="<<nQm<<",rN="
04038 <<r4M.m()<<"="<<G4QNucleus(rPDG).GetGSMass()<<G4endl;
04039 for(G4int i=0; i<nHd; i++)
04040 {
04041 G4LorentzVector hI4M=(*theResult)[i]->Get4Momentum();
04042 f4M+=hI4M;
04043 G4int hChg=(*theResult)[i]->GetCharge();
04044 fCh-=hChg;
04045 G4int hBaN=(*theResult)[i]->GetBaryonNumber();
04046 fBN-=hBaN;
04047 G4cout<<"-EMCLS-G4QFragmentation::Breeder:(2) Hadron#"<<i<<", 4M="<<hI4M<<", PDG="
04048 <<(*theResult)[i]->GetPDGCode()<<", C="<<hChg<<", B="<<hBaN<<G4endl;
04049 }
04050 for(G4int i=0; i<nQm; i++)
04051 {
04052 G4LorentzVector hI4M=theQuasmons[i]->Get4Momentum();
04053 f4M+=hI4M;
04054 G4int hChg=theQuasmons[i]->GetCharge();
04055 fCh-=hChg;
04056 G4int hBaN=theQuasmons[i]->GetBaryonNumber();
04057 fBN-=hBaN;
04058 G4cout<<"-EMCLS-G4QFragmentation::Breeder:(2) Quasmon#"<<i<<", 4M="<<hI4M<<", C="
04059 <<hChg<<", B="<<hBaN<<G4endl;
04060 }
04061 G4cout<<"-EMCLS-G4QFragm::Breed:, r4M="<<f4M-totLS4M<<",rC="<<fCh<<",rB="<<fBN<<G4endl;
04062 #endif
04063 }
04064 return;
04065 }
04066
04067
04068 G4bool G4QFragmentation::ExciteDiffParticipants(G4QHadron* projectile,
04069 G4QHadron* target) const
04070 {
04071 G4LorentzVector Pprojectile=projectile->Get4Momentum();
04072 G4double Mprojectile=projectile->GetMass();
04073 G4double Mprojectile2=Mprojectile*Mprojectile;
04074 G4LorentzVector Ptarget=target->Get4Momentum();
04075 G4double Mtarget=target->GetMass();
04076 G4double Mtarget2=Mtarget*Mtarget;
04077 #ifdef debug
04078 G4cout<<"G4QFragm::ExciteDiffPartici:Ep="<<Pprojectile.e()<<",Et="<<Ptarget.e()<<G4endl;
04079 #endif
04080
04081 G4LorentzVector Psum=Pprojectile+Ptarget;
04082 G4LorentzRotation toCms(-Psum.boostVector());
04083 G4LorentzVector Ptmp=toCms*Pprojectile;
04084 if(Ptmp.pz()<=0.)
04085 {
04086 #ifdef debug
04087 G4cout<<"G4QFragmentation::ExciteDiffParticipants: *1* abort Collision!! *1*"<<G4endl;
04088 #endif
04089 return false;
04090 }
04091 toCms.rotateZ(-Ptmp.phi());
04092 toCms.rotateY(-Ptmp.theta());
04093 #ifdef debug
04094 G4cout<<"G4QFragment::ExciteDiffParticipantts:BeforBoost Pproj="<<Pprojectile<<", Ptarg="
04095 <<Ptarget<<G4endl;
04096 #endif
04097 G4LorentzRotation toLab(toCms.inverse());
04098 Pprojectile.transform(toCms);
04099 Ptarget.transform(toCms);
04100 #ifdef debug
04101 G4cout<< "G4QFragment::ExciteDiffParticipantts: AfterBoost Pproj="<<Pprojectile<<"Ptarg="
04102 <<Ptarget<<", cms4M="<<Pprojectile+Ptarget<<G4endl;
04103 G4cout<<"G4QFragment::ExciteDiffParticipants: ProjX+="<<Pprojectile.plus()<<", ProjX-="
04104 <<Pprojectile.minus()<<", TargX+="<< Ptarget.plus()<<", TargX-="<<Ptarget.minus()
04105 <<G4endl;
04106 #endif
04107 G4LorentzVector Qmomentum(0.,0.,0.,0.);
04108 G4int whilecount=0;
04109 #ifdef debug
04110 G4cout<<"G4QFragmentation::ExciteDiffParticipants: Before DO"<<G4endl;
04111 #endif
04112 do
04113 {
04114
04115 G4double maxPtSquare=sqr(Ptarget.pz());
04116 #ifdef debug
04117 G4cout<<"G4QFragmentation::ExciteDiffParticipants: maxPtSq="<<maxPtSquare<<G4endl;
04118 if(whilecount++>=500 && whilecount%100==0)
04119 G4cout<<"G4QFragmentation::ExciteDiffParticipantts: can loop, loopCount="<<whilecount
04120 <<", maxPtSquare="<<maxPtSquare<<G4endl;
04121 #endif
04122 if(whilecount>1000)
04123 {
04124 #ifdef debug
04125 G4cout<<"G4QFragmentation::ExciteDiffParticipants: *2* abort Loop!! *2*"<<G4endl;
04126 #endif
04127 return false;
04128 }
04129 Qmomentum=G4LorentzVector(GaussianPt(widthOfPtSquare,maxPtSquare),0);
04130 #ifdef debug
04131 G4cout<<"G4QFragment::ExciteDiffParticipants: generated Pt="<<Qmomentum<<", ProjPt="
04132 <<Pprojectile+Qmomentum<<", TargPt="<<Ptarget-Qmomentum<<G4endl;
04133 #endif
04134
04135 G4double Xmin=0.;
04136 G4double Xmax=1.;
04137 G4double Xplus =ChooseX(Xmin,Xmax);
04138 G4double Xminus=ChooseX(Xmin,Xmax);
04139 #ifdef debug
04140 G4cout<<"G4QFragment::ExciteDiffParticip: X-plus="<<Xplus<<",X-minus="<<Xminus<<G4endl;
04141 #endif
04142 G4double pt2=Qmomentum.vect().mag2();
04143 G4double Qplus =-pt2/Xminus/Ptarget.minus();
04144 G4double Qminus= pt2/Xplus /Pprojectile.plus();
04145 Qmomentum.setPz((Qplus-Qminus)/2);
04146 Qmomentum.setE( (Qplus+Qminus)/2);
04147 #ifdef debug
04148 G4cout<<"G4QFragment::ExciteDiffParticip: Qplus="<<Qplus<<", Qminus="<<Qminus<<", pt2="
04149 <<pt2<<", Qmomentum="<<Qmomentum<<", ProjM="<<(Pprojectile+Qmomentum).mag()
04150 <<", TargM="<<(Ptarget-Qmomentum).mag()<<G4endl;
04151 #endif
04152 } while((Pprojectile+Qmomentum).mag2()<=Mprojectile2 ||
04153 (Ptarget-Qmomentum).mag2()<=Mtarget2);
04154 Pprojectile += Qmomentum;
04155 Ptarget -= Qmomentum;
04156 #ifdef debug
04157 G4cout<<"G4QFragment::ExciteDiffParticipan: Proj(Q)="<<Pprojectile<<", Targ(Q)="<<Ptarget
04158 <<", Proj(back)="<<toLab*Pprojectile<<", Targ(bac)="<< toLab*Ptarget << G4endl;
04159 #endif
04160
04161 Pprojectile.transform(toLab);
04162 Ptarget.transform(toLab);
04163 #ifdef debug
04164 G4cout<< "G4QFragmentation::ExciteDiffParticipants: TargetMass="<<Ptarget.mag()<<G4endl;
04165 #endif
04166 target->Set4Momentum(Ptarget);
04167 #ifdef debug
04168 G4cout<<"G4QFragment::ExciteDiffParticipants:ProjectileMass="<<Pprojectile.mag()<<G4endl;
04169 #endif
04170 projectile->Set4Momentum(Pprojectile);
04171 return true;
04172 }
04173
04174
04175
04176 G4bool G4QFragmentation::ExciteSingDiffParticipants(G4QHadron* projectile,
04177 G4QHadron* target) const
04178 {
04179 G4LorentzVector Pprojectile=projectile->Get4Momentum();
04180 G4double Mprojectile=projectile->GetMass();
04181 G4double Mprojectile2=Mprojectile*Mprojectile;
04182 G4LorentzVector Ptarget=target->Get4Momentum();
04183 G4double Mtarget=target->GetMass();
04184 G4double Mtarget2=Mtarget*Mtarget;
04185 #ifdef debug
04186 G4cout<<"G4QFragm::ExSingDiffPartici:Ep="<<Pprojectile.e()<<",Et="<<Ptarget.e()<<G4endl;
04187 #endif
04188 G4bool KeepProjectile= G4UniformRand() > 0.5;
04189
04190 if(KeepProjectile )
04191 {
04192 #ifdef debug
04193 G4cout<<"--1/2--G4QFragmentation::ExSingDiffParticipants: Projectile is fixed"<<G4endl;
04194 #endif
04195 Mprojectile2 = projectile->GetMass2()*(1.-perCent);
04196 }
04197 else
04198 {
04199 #ifdef debug
04200 G4cout<<"---1/2---G4QFragmentation::ExSingDiffParticipants: Target is fixed"<<G4endl;
04201 #endif
04202 Mtarget2 = target->GetMass2()*(1.-perCent);
04203 }
04204
04205
04206 G4LorentzVector Psum=Pprojectile+Ptarget;
04207 G4LorentzRotation toCms(-Psum.boostVector());
04208 G4LorentzVector Ptmp=toCms*Pprojectile;
04209 if(Ptmp.pz()<=0.)
04210 {
04211 #ifdef debug
04212 G4cout<<"G4QFragment::ExciteSingDiffParticipants: *1* abort Collision!! *1*"<<G4endl;
04213 #endif
04214 return false;
04215 }
04216 toCms.rotateZ(-Ptmp.phi());
04217 toCms.rotateY(-Ptmp.theta());
04218 #ifdef debug
04219 G4cout<<"G4QFragm::ExciteSingDiffParticipantts: Be4Boost Pproj="<<Pprojectile<<",Ptarg="
04220 <<Ptarget<<G4endl;
04221 #endif
04222 G4LorentzRotation toLab(toCms.inverse());
04223 Pprojectile.transform(toCms);
04224 Ptarget.transform(toCms);
04225 #ifdef debug
04226 G4cout<< "G4QFragment::ExciteDiffParticipantts: AfterBoost Pproj="<<Pprojectile<<"Ptarg="
04227 <<Ptarget<<", cms4M="<<Pprojectile+Ptarget<<G4endl;
04228
04229 G4cout<<"G4QFragment::ExciteDiffParticipantts: ProjX+="<<Pprojectile.plus()<<", ProjX-="
04230 <<Pprojectile.minus()<<", TargX+="<< Ptarget.plus()<<", TargX-="<<Ptarget.minus()
04231 <<G4endl;
04232 #endif
04233 G4LorentzVector Qmomentum(0.,0.,0.,0.);
04234 G4int whilecount=0;
04235 do
04236 {
04237
04238 G4double maxPtSquare=sqr(Ptarget.pz());
04239 if(whilecount++>=500 && whilecount%100==0)
04240 #ifdef debug
04241 G4cout<<"G4QFragment::ExciteSingDiffParticipantts: can loop, loopCount="<<whilecount
04242 <<", maxPtSquare="<<maxPtSquare<<G4endl;
04243 #endif
04244 if(whilecount>1000)
04245 {
04246 #ifdef debug
04247 G4cout<<"G4QFragmentation::ExciteSingDiffParticipants: *2* abort Loop!! *2*"<<G4endl;
04248 #endif
04249 return false;
04250 }
04251 Qmomentum=G4LorentzVector(GaussianPt(widthOfPtSquare,maxPtSquare),0);
04252 #ifdef debug
04253 G4cout<<"G4QFragm::ExciteSingDiffParticipants: generated Pt="<<Qmomentum<<", ProjPt="
04254 <<Pprojectile+Qmomentum<<", TargPt="<<Ptarget-Qmomentum<<G4endl;
04255 #endif
04256
04257 G4double Xmin=0.;
04258 G4double Xmax=1.;
04259 G4double Xplus =ChooseX(Xmin,Xmax);
04260 G4double Xminus=ChooseX(Xmin,Xmax);
04261 #ifdef debug
04262 G4cout<<"G4QFragm::ExciteSingDiffPartici:X-plus="<<Xplus<<",X-minus="<<Xminus<<G4endl;
04263 #endif
04264 G4double pt2=G4ThreeVector(Qmomentum.vect()).mag2();
04265 G4double Qplus =-pt2/Xminus/Ptarget.minus();
04266 G4double Qminus= pt2/Xplus /Pprojectile.plus();
04267 if (KeepProjectile)
04268 Qminus=(projectile->GetMass2()+pt2)/(Pprojectile.plus()+Qplus) - Pprojectile.minus();
04269 else Qplus=Ptarget.plus() - (target->GetMass2()+pt2)/(Ptarget.minus()-Qminus);
04270 Qmomentum.setPz((Qplus-Qminus)/2);
04271 Qmomentum.setE( (Qplus+Qminus)/2);
04272 #ifdef debug
04273 G4cout<<"G4QFragm::ExciteDiffParticip: Qplus="<<Qplus<<", Qminus="<<Qminus<<", pt2="
04274 <<pt2<<", Qmomentum="<<Qmomentum<<", ProjM="<<(Pprojectile+Qmomentum).mag()
04275 <<", TargM="<<(Ptarget-Qmomentum).mag()<<G4endl;
04276 #endif
04277
04278
04279
04280 } while((Ptarget-Qmomentum).mag2()<=Mtarget2 ||
04281 (Pprojectile+Qmomentum).mag2()<=Mprojectile2 ||
04282 (Ptarget-Qmomentum).e() < 0. || (Pprojectile+Qmomentum).e() < 0.);
04283 Pprojectile += Qmomentum;
04284 Ptarget -= Qmomentum;
04285 #ifdef debug
04286 G4cout<<"G4QFragmentation::ExciteSingDiffParticipan: Proj(Q)="<<Pprojectile<<"(E="
04287 <<Pprojectile.e()<<"), Targ(Q)="<<Ptarget<<"(E="<<Ptarget.e()
04288 <<"), Proj(back)="<<toLab*Pprojectile<<", Targ(bac)="<< toLab*Ptarget << G4endl;
04289 #endif
04290
04291 Pprojectile.transform(toLab);
04292 Ptarget.transform(toLab);
04293 #ifdef debug
04294 G4cout<< "G4QFragm::ExciteSingleDiffParticipants: TargetMass="<<Ptarget.mag()<<G4endl;
04295 #endif
04296 target->Set4Momentum(Ptarget);
04297 #ifdef debug
04298 G4cout<<"G4QFragm::ExciteSingleParticipants:ProjectileMass="<<Pprojectile.mag()<<G4endl;
04299 #endif
04300 projectile->Set4Momentum(Pprojectile);
04301 return true;
04302 }
04303
04304 void G4QFragmentation::SetParameters(G4int nC, G4double stT, G4double tbD, G4double SigPt)
04305 {
04306 nCutMax = nC;
04307 stringTension = stT;
04308 tubeDensity = tbD;
04309 widthOfPtSquare = -2*SigPt*SigPt;
04310 }
04311
04312 G4double G4QFragmentation::ChooseX(G4double Xmin, G4double Xmax) const
04313 {
04314
04315
04316 if(Xmax == Xmin) return Xmin;
04317 if( Xmin < 0. || Xmax < Xmin)
04318 {
04319 G4cerr<<"***G4QFragmentation::ChooseX: Xmin="<<Xmin<<", Xmax="<<Xmax<< G4endl;
04320 G4Exception("G4QFragmentation::ChooseX:","72",FatalException,"Bad X or X-Range");
04321 }
04322
04323
04324 G4double sxi=std::sqrt(Xmin);
04325 G4double x=sqr(sxi+G4UniformRand()*(std::sqrt(Xmax)-sxi));
04326 #ifdef debug
04327 G4cout<<"G4QFragmentation::ChooseX: DiffractiveX="<<x<<G4endl;
04328 #endif
04329 return x;
04330 }
04331
04332
04333 G4ThreeVector G4QFragmentation::GaussianPt(G4double widthSq, G4double maxPtSquare) const
04334 {
04335 #ifdef debug
04336 G4cout<<"G4QFragmentation::GaussianPt:width2="<<widthSq<<",maxPt2="<<maxPtSquare<<G4endl;
04337 #endif
04338 G4double pt2=0.;
04339 G4double rm=maxPtSquare/widthSq;
04340 if(rm>-.01) pt2=widthSq*(std::sqrt(1.-G4UniformRand()*(1.-sqr(1.+rm)))-1.);
04341 else pt2=widthSq*std::log(1.-G4UniformRand()*(1.-std::exp(rm)));
04342 pt2=std::sqrt(pt2);
04343 G4double phi=G4UniformRand()*twopi;
04344 return G4ThreeVector(pt2*std::cos(phi),pt2*std::sin(phi),0.);
04345 }
04346
04347 G4int G4QFragmentation::SumPartonPDG(G4int PDG1, G4int PDG2) const
04348 {
04349 if (PDG1 < 7 && PDG1 > 0 && PDG2 < 7 && PDG2 > 0)
04350 {
04351 if(PDG1 > PDG2) return PDG1*1000+PDG2*100+1;
04352 else return PDG2*1000+PDG1*100+1;
04353 }
04354 else if (PDG1 >-7 && PDG1 < 0 && PDG2 >-7 && PDG2 < 0)
04355 {
04356 if(-PDG1 > -PDG2) return PDG1*1000+PDG2*100-1;
04357 else return PDG2*1000+PDG1*100-1;
04358 }
04359 else if (PDG1 <-99 && PDG2 < 7 && PDG2 > 0)
04360 {
04361 G4int PDG=-PDG1/100;
04362 if(PDG2==PDG/10) return -PDG%10;
04363 if(PDG2==PDG%10) return -PDG/10;
04364 else
04365 {
04366 G4cerr<<"***4QFragmentation::SumPartonPDG: PDG1="<<PDG1<<", PDG2="<<PDG2<<G4endl;
04367 G4Exception("G4QFragmentation::SumPartonPDG:","72",FatalException,"Q&ADiQ notMatch");
04368 }
04369 }
04370 else if (PDG2 <-99 && PDG1 < 7 && PDG1 > 0)
04371 {
04372 G4int PDG=-PDG2/100;
04373 if(PDG1==PDG/10) return -PDG%10;
04374 if(PDG1==PDG%10) return -PDG/10;
04375 else
04376 {
04377 G4cerr<<"***4QFragmentation::SumPartonPDG: PDG1="<<PDG1<<", PDG2="<<PDG2<<G4endl;
04378 G4Exception("G4QFragmentation::SumPartonPDG:","72",FatalException,"ADiQ&Q notMatch");
04379 }
04380 }
04381 else if (PDG1 > 99 && PDG2 >-7 && PDG2 < 0)
04382 {
04383 G4int PDG=PDG1/100;
04384 if(PDG2==-PDG/10) return PDG%10;
04385 if(PDG2==-PDG%10) return PDG/10;
04386 else
04387 {
04388 G4cerr<<"***4QFragmentation::SumPartonPDG: PDG1="<<PDG1<<", PDG2="<<PDG2<<G4endl;
04389 G4Exception("G4QFragmentation::SumPartonPDG:","72",FatalException,"DiQ&AQ notMatch");
04390 }
04391 }
04392 else if (PDG2 > 99 && PDG1 >-7 && PDG1 < 0)
04393 {
04394 G4int PDG=PDG2/100;
04395 if(PDG1==-PDG/10) return PDG%10;
04396 if(PDG1==-PDG%10) return PDG/10;
04397 else
04398 {
04399 G4cerr<<"***G4QFragmentation::SumPartonPDG: PDG1="<<PDG1<<", PDG2="<<PDG2<<G4endl;
04400 G4Exception("G4QFragmentation::SumPartonPDG:","72",FatalException,"AQ&DiQ notMatch");
04401 }
04402 }
04403 else
04404 {
04405 G4cerr<<"***G4QFragmentation::SumPartonPDG: PDG1="<<PDG1<<", PDG2="<<PDG2<<G4endl;
04406 G4Exception("G4QFragmentation::SumPartonPDG:","72",FatalException,"Can'tSumUpPartons");
04407 }
04408 return 0;
04409 }
04410
04411
04412 std::pair<G4int,G4int> G4QFragmentation::ReducePair(G4int P1, G4int P2) const
04413 {
04414 #ifdef debug
04415 G4cout<<"G4QFragmentation::ReducePair: **Called** P1="<<P1<<", P2="<<P2<<G4endl;
04416 #endif
04417 G4int P11=P1/10;
04418 G4int P12=P1%10;
04419 G4int P21=P2/10;
04420 G4int P22=P2%10;
04421 if (P11==P21) return std::make_pair(P12,P22);
04422 else if(P11==P22) return std::make_pair(P12,P21);
04423 else if(P12==P21) return std::make_pair(P11,P22);
04424 else if(P12==P22) return std::make_pair(P11,P21);
04425
04426 G4cout<<"-Warning-G4QFragmentation::ReducePair:**Failed**, P1="<<P1<<", P2="<<P2<<G4endl;
04427
04428 return std::make_pair(0,0);
04429 }
04430
04431
04432 G4int G4QFragmentation::AnnihilationOrder(G4int LS, G4int MPS, G4int uPDG, G4int mPDG,
04433 G4int sPDG, G4int nPDG)
04434 {
04435 G4int Ord=0;
04436
04437 if (LS==2 && MPS==2 )
04438 {
04439 #ifdef debug
04440 G4cout<<"G4QFragmentation::AnnihilationOrder: QaQ/QaQ->DiQ-aDiQ, uPDG="<<uPDG<<G4endl;
04441 #endif
04442 if ( (uPDG>0 && sPDG>0 && mPDG<0 && nPDG<0) ||
04443 (uPDG<0 && sPDG<0 && mPDG>0 && nPDG>0) ) Ord= 1;
04444 else if( (uPDG>0 && nPDG>0 && mPDG<0 && sPDG<0) ||
04445 (uPDG<0 && nPDG<0 && mPDG>0 && sPDG>0) ) Ord=-1;
04446 else G4cerr<<"-Warning-G4QFragmentation::AnnihilationOrder: Wrong 22 fusion, L="<<uPDG
04447 <<",R="<<mPDG<<",cL="<<sPDG<<",cR="<<nPDG<<G4endl;
04448 }
04449 else if ( LS==2 && MPS==3 )
04450 {
04451 if (uPDG > 7)
04452 {
04453 #ifdef debug
04454 G4cout<<"G4QFragmentation::AnnihOrder: pLDiQ, sPDG="<<sPDG<<", nPDG="<<nPDG<<G4endl;
04455 #endif
04456 if ( sPDG<0 && (-sPDG==uPDG/1000 || -sPDG==(uPDG/100)%10) ) Ord= 1;
04457 else if( nPDG<0 && (-nPDG==uPDG/1000 || -nPDG==(uPDG/100)%10) ) Ord=-1;
04458 else G4cerr<<"-Warning-G4QFragmentation::AnnihilationOrder: Wrong pLDiQ, L="<<uPDG
04459 <<",R="<<mPDG<<",cL="<<sPDG<<",cR="<<nPDG<<G4endl;
04460 }
04461 else if (mPDG > 7)
04462 {
04463 #ifdef debug
04464 G4cout<<"G4QFragmentation::AnnihOrder: pRDiQ, sPDG="<<sPDG<<", nPDG="<<nPDG<<G4endl;
04465 #endif
04466 if ( sPDG<0 && (-sPDG==mPDG/1000 || -sPDG==(mPDG/100)%10) ) Ord=-1;
04467 else if( nPDG<0 && (-nPDG==mPDG/1000 || -nPDG==(mPDG/100)%10) ) Ord= 1;
04468 else G4cerr<<"-Warning-G4QFragmentation::AnnihilationOrder: Wrong pRDiQ, L="<<uPDG
04469 <<",R="<<mPDG<<",cL="<<sPDG<<",cR="<<nPDG<<G4endl;
04470 }
04471 else if (uPDG <-7)
04472 {
04473 #ifdef debug
04474 G4cout<<"G4QFragmentation::AnnihOrder: pLaDiQ, sPDG="<<sPDG<<", nPDG="<<nPDG<<G4endl;
04475 #endif
04476 if ( sPDG>0 && (sPDG==(-uPDG)/1000 || sPDG==((-uPDG)/100)%10) ) Ord= 1;
04477 else if( nPDG>0 && (nPDG==(-uPDG)/1000 || nPDG==((-uPDG)/100)%10) ) Ord=-1;
04478 else G4cerr<<"-Warning-G4QFragmentation::AnnihilationOrder: Wrong pLaDiQ, L="<<uPDG
04479 <<", R="<<mPDG<<", cL="<<sPDG<<", cR="<<nPDG<<G4endl;
04480 }
04481 else if (mPDG <-7)
04482 {
04483 #ifdef debug
04484 G4cout<<"G4QFragmentation::AnnihOrder: pRaDiQ, sPDG="<<sPDG<<", nPDG="<<nPDG<<G4endl;
04485 #endif
04486 if ( sPDG>0 && (sPDG==(-mPDG)/1000 || sPDG==((-mPDG)/100)%10) ) Ord=-1;
04487 else if( nPDG>0 && (nPDG==(-mPDG)/1000 || nPDG==((-mPDG)/100)%10) ) Ord= 1;
04488 else G4cerr<<"-Warning-G4QFragmentation::AnnihilationOrder: Wrong pRaDiQ, L="<<uPDG
04489 <<", R="<<mPDG<<", cL="<<sPDG<<", cR="<<nPDG<<G4endl;
04490 }
04491 else if( (sPDG<0 && (-sPDG==mPDG || -sPDG==uPDG) ) ||
04492 (nPDG<0 && (-nPDG==mPDG || -nPDG==uPDG) ) ) Ord= 2;
04493 #ifdef debug
04494 else G4cout<<"-Warning-G4QFragmentation::AnnihilatOrder: Wrong 23StringFusion"<<G4endl;
04495 G4cout<<"G4QFragmentation::AnnihilationOrder: Ord="<<Ord<<",sPDG="<<sPDG<<",nPDG="
04496 <<nPDG<<", uPDG="<<uPDG<<",mPDG="<<mPDG<<G4endl;
04497 #endif
04498 }
04499 else if ( LS==3 && MPS==2 )
04500 {
04501 if (sPDG > 7)
04502 {
04503 #ifdef debug
04504 G4cout<<"G4QFragmentation::AnnihOrder: cLDiQ, uPDG="<<uPDG<<", mPDG="<<mPDG<<G4endl;
04505 #endif
04506 if ( uPDG<0 && (-uPDG==sPDG/1000 || -uPDG==(sPDG/100)%10) ) Ord= 1;
04507 else if( mPDG<0 && (-mPDG==sPDG/1000 || -mPDG==(sPDG/100)%10) ) Ord=-1;
04508 else G4cerr<<"-Warning-G4QFragmentation::AnnihilationOrder: Wrong cLDiQ, L="<<uPDG
04509 <<",R="<<mPDG<<",cL="<<sPDG<<",cR="<<nPDG<<G4endl;
04510 }
04511 else if (nPDG > 7)
04512 {
04513 #ifdef debug
04514 G4cout<<"G4QFragmentation::AnnihOrder: cRDiQ, uPDG="<<uPDG<<", mPDG="<<mPDG<<G4endl;
04515 #endif
04516 if ( uPDG<0 && (-uPDG==nPDG/1000 || -uPDG==(nPDG/100)%10) ) Ord=-1;
04517 else if( mPDG<0 && (-mPDG==nPDG/1000 || -mPDG==(nPDG/100)%10) ) Ord= 1;
04518 else G4cerr<<"-Warning-G4QFragmentation::AnnihilationOrder: Wrong cRDiQ, L="<<uPDG
04519 <<",R="<<mPDG<<",cL="<<sPDG<<",cR="<<nPDG<<G4endl;
04520 }
04521 else if (sPDG <-7)
04522 {
04523 #ifdef debug
04524 G4cout<<"G4QFragmentation::AnnihOrder: cLaDiQ, uPDG="<<uPDG<<", mPDG="<<mPDG<<G4endl;
04525 #endif
04526 if ( uPDG>0 && (uPDG==(-sPDG)/1000 || uPDG==((-sPDG)/100)%10) ) Ord= 1;
04527 else if( mPDG>0 && (mPDG==(-sPDG)/1000 || mPDG==((-sPDG)/100)%10) ) Ord=-1;
04528 else G4cerr<<"-Warning-G4QFragmentation::AnnihilationOrder: Wrong cLaDiQ, L="<<uPDG
04529 <<", R="<<mPDG<<", cL="<<sPDG<<", cR="<<nPDG<<G4endl;
04530 }
04531 else if (nPDG <-7)
04532 {
04533 #ifdef debug
04534 G4cout<<"G4QFragmentation::AnnihOrder: cRaDiQ, uPDG="<<uPDG<<", mPDG="<<mPDG<<G4endl;
04535 #endif
04536 if ( uPDG>0 && (uPDG==(-nPDG)/1000 || uPDG==((-nPDG)/100)%10) ) Ord=-1;
04537 else if( mPDG>0 && (mPDG==(-nPDG)/1000 || mPDG==((-nPDG)/100)%10) ) Ord= 1;
04538 else G4cerr<<"-Warning-G4QFragmentation::AnnihilationOrder: Wrong cRaDiQ, L="<<uPDG
04539 <<", R="<<mPDG<<", cL="<<sPDG<<", cR="<<nPDG<<G4endl;
04540 }
04541 else if( (uPDG<0 && (-uPDG==sPDG || -uPDG==nPDG) ) ||
04542 (mPDG<0 && (-mPDG==sPDG || -mPDG==nPDG) ) ) Ord=2;
04543 #ifdef debug
04544 else G4cout<<"-Warning-G4QFragmentation::AnnihilatOrder: Wrong 32StringFusion"<<G4endl;
04545 G4cout<<"G4QFragmentation::AnnihilationOrder: Ord="<<Ord<<",sPDG="<<sPDG<<",nPDG="
04546 <<nPDG<<", uPDG="<<uPDG<<",mPDG="<<mPDG<<G4endl;
04547 #endif
04548 }
04549 else if ( (LS==2 && MPS==4) || (LS==4 && MPS==2) )
04550 {
04551 if (uPDG > 7)
04552 {
04553 #ifdef debug
04554 G4cout<<"G4QFragmentation::AnnihilOrder:pLDiQ, sPDG="<<sPDG<<",nPDG="<<nPDG<<G4endl;
04555 #endif
04556 if ( sPDG<0 && (-sPDG==uPDG/1000 || -sPDG==(uPDG/100)%10) &&
04557 (nPDG==(-mPDG)/1000 || nPDG==((-mPDG)/100)%10) ) Ord= 1;
04558 else if( nPDG<0 && (-nPDG==uPDG/1000 || -nPDG==(uPDG/100)%10) &&
04559 (sPDG==(-mPDG)/1000 || sPDG==((-mPDG)/100)%10) ) Ord=-1;
04560 else G4cerr<<"-Warning-G4QFragmentation::AnnihilationOrder: Wrong 24 pLDiQ, L="<<uPDG
04561 <<",R="<<mPDG<<",cL="<<sPDG<<",cR="<<nPDG<<G4endl;
04562 }
04563 else if (mPDG >7)
04564 {
04565 #ifdef debug
04566 G4cout<<"G4QFragmentation::AnnihilOrder:PRDiQ, sPDG="<<sPDG<<",nPDG="<<nPDG<<G4endl;
04567 #endif
04568 if ( sPDG<0 && (-sPDG==mPDG/1000 || -sPDG==(mPDG/100)%10) &&
04569 (nPDG==(-uPDG)/1000 || nPDG==((-uPDG)/100)%10) ) Ord=-1;
04570 else if( nPDG<0 && (-nPDG==mPDG/1000 || -nPDG==(mPDG/100)%10) &&
04571 (sPDG==(-uPDG)/1000 || sPDG==((-uPDG)/100)%10) ) Ord= 1;
04572 else G4cerr<<"-Warning-G4QFragmentation::AnnihilationOrder: Wrong 24 pLDiQ, L="<<uPDG
04573 <<",R="<<mPDG<<",cL="<<sPDG<<",cR="<<nPDG<<G4endl;
04574 }
04575 else if (sPDG > 7)
04576 {
04577 #ifdef debug
04578 G4cout<<"G4QFragmentation::AnnihilOrder:cLDiQ, uPDG="<<uPDG<<",mPDG="<<mPDG<<G4endl;
04579 #endif
04580 if ( uPDG<0 && (-uPDG==sPDG/1000 || -uPDG==(sPDG/100)%10) &&
04581 (mPDG==(-nPDG)/1000 || mPDG==((-nPDG)/100)%10) ) Ord= 1;
04582 else if( mPDG<0 && (-mPDG==sPDG/1000 || -mPDG==(sPDG/100)%10) &&
04583 (uPDG==(-nPDG)/1000 || uPDG==((-nPDG)/100)%10) ) Ord=-1;
04584 else G4cerr<<"-Warning-G4QFragmentation::AnnihilationOrder: Wrong 24 cLDiQ, L="<<uPDG
04585 <<",R="<<mPDG<<",cL="<<sPDG<<",cR="<<nPDG<<G4endl;
04586 }
04587 else if (nPDG > 7)
04588 {
04589 #ifdef debug
04590 G4cout<<"G4QFragmentation::AnnihilOrder:cRDiQ, uPDG="<<uPDG<<",mPDG="<<mPDG<<G4endl;
04591 #endif
04592 if ( uPDG<0 && (-uPDG==nPDG/1000 || -uPDG==(nPDG/100)%10) &&
04593 (mPDG==(-sPDG)/1000 || mPDG==((-sPDG)/100)%10) ) Ord=-1;
04594 else if( mPDG<0 && (-mPDG==nPDG/1000 || -mPDG==(nPDG/100)%10) &&
04595 (uPDG==(-sPDG)/1000 || uPDG==((-sPDG)/100)%10) ) Ord= 1;
04596 else G4cerr<<"-Warning-G4QFragmentation::AnnihilationOrder: Wrong 24 cRDiQ, L="<<uPDG
04597 <<",R="<<mPDG<<",cL="<<sPDG<<",cR="<<nPDG<<G4endl;
04598 }
04599 #ifdef debug
04600 else G4cout<<"-Warning-G4QFragmentation::AnnihilOrder: Wrong 24 StringFusion"<<G4endl;
04601 G4cout<<"G4QFragmentation::AnnihilationOrder: Ord="<<Ord<<",sPDG="<<sPDG<<",nPDG="
04602 <<nPDG<<", uPDG="<<uPDG<<",mPDG="<<mPDG<<G4endl;
04603 #endif
04604 }
04605 else if ( LS==3 && MPS==3 )
04606 {
04607 if (uPDG > 7)
04608 {
04609 #ifdef debug
04610 G4cout<<"G4QFragmentation::AnnihilOrder: pLDiQ, sPDG="<<sPDG<<",nPDG="<<nPDG<<G4endl;
04611 #endif
04612 if ( sPDG<-7 && (-nPDG==uPDG/1000 || -nPDG==(uPDG/100)%10) &&
04613 (mPDG==(-sPDG)/1000 || mPDG==((-sPDG)/100)%10) ) Ord=-1;
04614 else if( nPDG<-7 && (-sPDG==uPDG/1000 || -sPDG==(uPDG/100)%10) &&
04615 (mPDG==(-nPDG)/1000 || mPDG==((-nPDG)/100)%10) ) Ord= 1;
04616 else G4cerr<<"-Warning-G4QFragmentation::AnnihilationOrder: Wrong 33 pLDiQ, L="<<uPDG
04617 <<",R="<<mPDG<<",cL="<<sPDG<<",cR="<<nPDG<<G4endl;
04618 }
04619 else if(mPDG > 7)
04620 {
04621 #ifdef debug
04622 G4cout<<"G4QFragmentation::AnnihilOrder: pRDiQ, sPDG="<<sPDG<<",nPDG="<<nPDG<<G4endl;
04623 #endif
04624 if ( sPDG<-7 && (-nPDG==mPDG/1000 || -nPDG==(mPDG/100)%10) &&
04625 (uPDG==(-sPDG)/1000 || uPDG==((-sPDG)/100)%10) ) Ord= 1;
04626 else if( nPDG<-7 && (-sPDG==mPDG/1000 || -sPDG==(mPDG/100)%10) &&
04627 (uPDG==(-nPDG)/1000 || uPDG==((-nPDG)/100)%10) ) Ord=-1;
04628 else G4cerr<<"-Warning-G4QFragmentation::AnnihilationOrder: Wrong 33 pRDiQ, L="<<uPDG
04629 <<",R="<<mPDG<<",cL="<<sPDG<<",cR="<<nPDG<<G4endl;
04630 }
04631 else if(sPDG > 7)
04632 {
04633 #ifdef debug
04634 G4cout<<"G4QFragmentation::AnnihilOrder: cLDiQ, uPDG="<<uPDG<<",mPDG="<<mPDG<<G4endl;
04635 #endif
04636 if ( uPDG<-7 && (-mPDG==sPDG/1000 || -mPDG==(sPDG/100)%10) &&
04637 (nPDG==(-uPDG)/1000 || nPDG==((-uPDG)/100)%10) ) Ord=-1;
04638 else if( mPDG<-7 && (-uPDG==sPDG/1000 || -uPDG==(sPDG/100)%10) &&
04639 (nPDG==(-mPDG)/1000 || nPDG==((-mPDG)/100)%10) ) Ord= 1;
04640 else G4cerr<<"-Warning-G4QFragmentation::AnnihilationOrder: Wrong 33 cLDiQ, L="<<uPDG
04641 <<",R="<<mPDG<<",cL="<<sPDG<<",cR="<<nPDG<<G4endl;
04642 }
04643 else if(nPDG > 7)
04644 {
04645 #ifdef debug
04646 G4cout<<"G4QFragmentation::AnnihilOrder: cRDiQ, uPDG="<<uPDG<<",mPDG="<<mPDG<<G4endl;
04647 #endif
04648 if ( uPDG<-7 && (-mPDG==nPDG/1000 || -mPDG==(nPDG/100)%10) &&
04649 (nPDG==(-uPDG)/1000 || nPDG==((-uPDG)/100)%10) ) Ord= 1;
04650 else if( mPDG<-7 && (-uPDG==nPDG/1000 || -sPDG==(nPDG/100)%10) &&
04651 (sPDG==(-mPDG)/1000 || sPDG==((-mPDG)/100)%10) ) Ord=-1;
04652 else G4cerr<<"-Warning-G4QFragmentation::AnnihilationOrder: Wrong 33 cRDiQ, L="<<uPDG
04653 <<",R="<<mPDG<<",cL="<<sPDG<<",cR="<<nPDG<<G4endl;
04654 }
04655 #ifdef debug
04656 else G4cout<<"-Warning-G4QFragmentation::AnnihilOrder: Wrong 33 StringFusion"<<G4endl;
04657 G4cout<<"G4QFragmentation::AnnihilationOrder: Ord="<<Ord<<",sPDG="<<sPDG<<",nPDG="
04658 <<nPDG<<", uPDG="<<uPDG<<",mPDG="<<mPDG<<G4endl;
04659 #endif
04660 }
04661 return Ord;
04662 }
04663
04664 void G4QFragmentation::SwapPartons()
04665 {
04666 static const G4double baryM=800.;
04667 G4QStringVector::iterator ist;
04668 for(ist = strings.begin(); ist < strings.end(); ist++)
04669 {
04670 G4QParton* cLeft=(*ist)->GetLeftParton();
04671 G4QParton* cRight=(*ist)->GetRightParton();
04672 G4LorentzVector cL4M=cLeft->Get4Momentum();
04673 G4LorentzVector cR4M=cRight->Get4Momentum();
04674 G4LorentzVector cS4M=cL4M+cR4M;
04675 G4double cSM2=cS4M.m2();
04676 if(std::fabs(cSM2)<.01)
04677 {
04678 G4double dM2=.001-cSM2;
04679 G4double E=cS4M.e();
04680 G4double dE=std::sqrt(E*E+dM2)-E;
04681 G4double LE=cL4M.e();
04682 G4double RE=cR4M.e();
04683 if(LE<RE) cLeft->Set4Momentum( G4LorentzVector(LE+dE) );
04684 else cRight->Set4Momentum( G4LorentzVector(RE+dE) );
04685 cSM2=.001;
04686 }
04687 if(cSM2<0.011)
04688 {
04689 G4int cLPDG=cLeft->GetPDGCode();
04690 G4int cRPDG=cRight->GetPDGCode();
04691 G4int cLT=cLeft->GetType();
04692 G4int cRT=cRight->GetType();
04693 G4QStringVector::iterator sst;
04694 G4QStringVector::iterator pst;
04695 G4double maxM=-DBL_MAX;
04696 G4int sDir=0;
04697 #ifdef debug
04698 G4cout<<"G4QFragmentation::SwapPartons: M2=="<<cSM2<<", 4M="<<cS4M<<",LPDG="<<cLPDG
04699 <<",RPDG="<<cRPDG<<G4endl;
04700 #endif
04701 for(pst = strings.begin(); pst < strings.end(); pst++) if(pst != ist)
04702 {
04703 G4QParton* pLeft=(*pst)->GetLeftParton();
04704 G4QParton* pRight=(*pst)->GetRightParton();
04705 G4int pLPDG=pLeft->GetPDGCode();
04706 G4int pRPDG=pRight->GetPDGCode();
04707 G4LorentzVector pL4M=pLeft->Get4Momentum();
04708 G4LorentzVector pR4M=pRight->Get4Momentum();
04709 G4int pLT=pLeft->GetType();
04710 G4int pRT=pRight->GetType();
04711 #ifdef debug
04712 G4cout<<"G4QFragmentation::SwapPartons: p4M="<<cS4M<<",pM2="<<cS4M.m2()<<",LPDG="
04713 <<pLPDG<<",RPDG="<<pRPDG<<G4endl;
04714 #endif
04715 G4double LM=0.;
04716 G4double RM=0.;
04717 if(((cLPDG<-7 || (cLPDG>0 && cLPDG< 7) ) && (pLPDG<-7 || (pLPDG>0 && pLPDG< 7) ))||
04718 ((cLPDG> 7 || (cLPDG<0 && cLPDG>-7) ) && (pLPDG> 7 || (pLPDG<0 && cLPDG>-7) )))
04719 {
04720 G4double pLM2=(cL4M+pR4M).m2();
04721 G4double cLM2=(cR4M+pL4M).m2();
04722 if(pLM2>0. && cLM2>0.)
04723 {
04724 G4double pLM=std::sqrt(pLM2);
04725 if(cLT+pRT==3) pLM-=baryM;
04726 G4double cLM=std::sqrt(cLM2);
04727 if(cRT+pLT==3) cLM-=baryM;
04728 LM=std::min(pLM2,cLM2);
04729 }
04730 }
04731 if(((cRPDG<-7 || (cRPDG>0 && cRPDG< 7) ) && (pRPDG<-7 || (pRPDG>0 && pRPDG< 7) ))||
04732 ((cRPDG> 7 || (cRPDG<0 && cRPDG>-7) ) && (pRPDG> 7 || (pRPDG<0 && cRPDG>-7) )) )
04733 {
04734 G4double pRM2=(cR4M+pL4M).m2();
04735 G4double cRM2=(cL4M+pR4M).m2();
04736 if(pRM2>0. && cRM2>0.)
04737 {
04738 G4double pRM=std::sqrt(pRM2);
04739 if(cRT+pLT==3) pRM-=baryM;
04740 G4double cRM=std::sqrt(cRM2);
04741 if(cLT+pRT==3) cRM-=baryM;
04742 RM=std::min(pRM,cRM);
04743 }
04744 }
04745 G4int dir=0;
04746 G4double sM=0.;
04747 if( LM && LM > RM )
04748 {
04749 dir= 1;
04750 sM=LM;
04751 }
04752 else if(RM)
04753 {
04754 dir=-1;
04755 sM=RM;
04756 }
04757 if(sM > maxM)
04758 {
04759 sst=pst;
04760 maxM=sM;
04761 sDir=dir;
04762 }
04763 }
04764 if(sDir)
04765 {
04766 G4QParton* pLeft=(*sst)->GetLeftParton();
04767 G4QParton* pRight=(*sst)->GetRightParton();
04768 G4QParton* swap=pLeft;
04769 if(sDir>0)
04770 {
04771 (*sst)->SetLeftParton(cLeft);
04772 (*ist)->SetLeftParton(swap);
04773 }
04774 else
04775 {
04776 swap=pRight;
04777 (*sst)->SetRightParton(cRight);
04778 (*ist)->SetRightParton(swap);
04779 }
04780 }
04781 #ifdef debug
04782 else G4cout<<"***G4QFragmentation::SwapPartons:**Failed**,cLPDG="<<cLPDG<<",cRPDG="
04783 <<cRPDG<<",-->cM2="<<cSM2<<G4endl;
04784 #endif
04785 }
04786 }
04787 }
04788
04789
04790 void G4QFragmentation::EvaporateResidual(G4QHadron* qH)
04791 {
04792 static const G4double mAlph = G4QPDGCode(2112).GetNuclMass(2,2,0);
04793 static const G4double mDeut = G4QPDGCode(2112).GetNuclMass(1,1,0);
04794 static const G4double mNeut = G4QPDGCode(2112).GetMass();
04795 static const G4double mProt = G4QPDGCode(2212).GetMass();
04796 static const G4double mAlPr = mAlph+mProt;
04797 static const G4double mAlNt = mAlph+mNeut;
04798 static const G4double dProt = mProt+mProt;
04799 static const G4double dNeut = mNeut+mNeut;
04800 static const G4double dAlph = mAlph+mAlph;
04801 static const G4double eps=.003;
04802 G4QEnvironment envir(theNucleus);
04803 G4int thePDG = qH->GetPDGCode();
04804 G4int theBN = qH->GetBaryonNumber();
04805 G4QContent theQC = qH->GetQC();
04806 G4int theS=theQC.GetStrangeness();
04807 #ifdef debug
04808 G4cout<<"G4QFragment::EvaRes:-Called- PDG="<<thePDG<<",4M="<<qH->Get4Momentum()
04809 <<",QC="<<theQC<<", BN="<<theBN<<G4endl;
04810 #endif
04811 if(thePDG==10)
04812 {
04813 #ifdef debug
04814 G4cout<<"G4QFragment::EvaRes: Cgipolino QC="<<theQC<<qH->Get4Momentum()<<G4endl;
04815 #endif
04816 G4QContent chQC=qH->GetQC();
04817 G4QChipolino QCh(chQC);
04818 G4LorentzVector ch4M=qH->Get4Momentum();
04819 G4QPDGCode h1QPDG=QCh.GetQPDG1();
04820 G4double h1M =h1QPDG.GetMass();
04821 G4QPDGCode h2QPDG=QCh.GetQPDG2();
04822 G4double h2M =h2QPDG.GetMass();
04823 G4double chM2 =ch4M.m2();
04824 if( sqr(h1M+h2M) < chM2 )
04825 {
04826 G4LorentzVector h14M(0.,0.,0.,h1M);
04827 G4LorentzVector h24M(0.,0.,0.,h2M);
04828 if(!G4QHadron(ch4M).DecayIn2(h14M,h24M))
04829 {
04830 G4cerr<<"***G4QFrag::EvaporateResid: CM="<<std::sqrt(chM2)<<" -> h1="<<h1QPDG<<"("
04831 <<h1M<<") + h2="<<h1QPDG<<"("<<h2M<<") = "<<h1M+h2M<<" **Failed**"<<G4endl;
04832
04833 G4Exception("G4QFragmentation::EvaporateResidual()", "HAD_CHPS_0000",
04834 FatalException, "QChipolino DecIn2 error");
04835 }
04836 delete qH;
04837 G4QHadron* h1H = new G4QHadron(h1QPDG.GetPDGCode(),h14M);
04838 theResult->push_back(h1H);
04839 #ifdef debug
04840 G4cout<<"G4QFragm::EvaporateResidual: Chipolino -> H1="<<h1QPDG<<h14M<<G4endl;
04841 #endif
04842 qH = new G4QHadron(h2QPDG.GetPDGCode(),h24M);
04843 theResult->push_back(qH);
04844 #ifdef debug
04845 G4cout<<"G4QE::EvaporateResidual: Chipolino -> H2="<<h2QPDG<<h24M<<G4endl;
04846 #endif
04847 }
04848 else
04849 {
04850 G4cerr<<"***G4QFragment::EvaporateResid: Chipolino="<<qH->GetQC()<<qH->Get4Momentum()
04851 <<", chipoM="<<std::sqrt(chM2)<<" < m1="<<h1M<<"("<<h1QPDG<<") + m2="<<h2M
04852 <<"("<<h2QPDG<<") = "<<h1M+h2M<<G4endl;
04853
04854 G4Exception("G4QFragmentation::EvaporateResidual()", "HAD_CHPS_0001",
04855 FatalException, "LowMassChipolino in Input");
04856 }
04857 return;
04858 }
04859 else if(theS<0)
04860 {
04861 #ifdef debug
04862 G4cout<<"G4QFragment::EvaRes: AntistrangeNucleus="<<thePDG<<qH->Get4Momentum()<<G4endl;
04863 #endif
04864 envir.DecayAntistrange(qH, theResult);
04865 return;
04866 }
04867 else if(theBN==1)
04868 {
04869 #ifdef debug
04870 G4cout<<"G4QFragmentation::EvaporateResid:Baryon="<<thePDG<<qH->Get4Momentum()<<G4endl;
04871 #endif
04872 envir.DecayBaryon(qH, theResult);
04873 return;
04874 }
04875 else if(!theBN)
04876 {
04877 #ifdef debug
04878 G4LorentzVector mesLV=qH->Get4Momentum();
04879 G4cout<<"G4QFragmentation::EvaRes:(!)Meson(!) PDG="<<thePDG<<",4M="<<mesLV<<mesLV.m()
04880 <<",QC="<<qH->GetQC()<<",MPDG="<<G4QPDGCode(thePDG).GetMass()<<G4endl;
04881 #endif
04882 envir.DecayMeson(qH, theResult);
04883 return;
04884 }
04885 G4int theC=theQC.GetCharge();
04886 #ifdef debug
04887 G4cout<<"G4QFragment::EvaRes: qH.Charge = "<<theC<<G4endl;
04888 #endif
04889 if(!thePDG) thePDG = theQC.GetSPDGCode();
04890 if( thePDG == 10 && theBN > 0 ) thePDG=theQC.GetZNSPDGCode();
04891 if(theS>0) thePDG-=theS*999999;
04892 #ifdef debug
04893 G4cout<<"G4QFragment::EvaRes: S="<<theS<<", newPDG="<<thePDG<<G4endl;
04894 #endif
04895 G4double totGSM = G4QNucleus(thePDG).GetGSMass();
04896 #ifdef debug
04897 G4cout<<"G4QFragment::EvaRes: totGSM="<<totGSM<<G4endl;
04898 #endif
04899 if(theBN==2)
04900 {
04901 if(!theC) totGSM=dNeut;
04902 else if(theC==2) totGSM=dProt;
04903 else totGSM=mDeut;
04904 }
04905 else if(theBN==5)
04906 {
04907 if (theC==3) totGSM=mAlPr;
04908 else if(theC==2) totGSM=mAlNt;
04909 }
04910 else if(theBN==8) totGSM=dAlph;
04911
04912 G4LorentzVector q4M = qH->Get4Momentum();
04913 G4double totMass = q4M.m();
04914 #ifdef debug
04915 G4cout<<"G4QFragment::EvaRes: Excitation = "<<totMass-totGSM<<G4endl;
04916 #endif
04917 if(std::fabs(totMass-totGSM) < eps)
04918 {
04919 theResult->push_back(qH);
04920 }
04921 else if(totMass > totGSM)
04922 {
04923 #ifdef debug
04924 G4cout<<"G4QFragment::EvaRes: try Evaporate Nucleus PDG="<<thePDG<<G4endl;
04925 #endif
04926 theNucleus.EvaporateNucleus(qH, theResult);
04927 #ifdef debug
04928 G4cout<<"G4QFragment::EvaRes: ** Evaporation is done **"<<G4endl;
04929 #endif
04930
04931 qH=0;
04932 }
04933 else
04934 {
04935 #ifdef debug
04936 G4cout<<"-War-G4QFr::EvaRes:*Must correct* "<<theQC<<q4M<<totMass<<"<"<<totGSM<<G4endl;
04937 #endif
04938 }
04939 #ifdef qdebug
04940 if (qH)
04941 {
04942 G4cout<<"-W-G4QFragmentation::EvaporateResid:*Deleted*,PDG="<<qH->GetPDGCode()<<G4endl;
04943 delete qH;
04944 }
04945 #endif
04946 return;
04947 }
