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 "G4QContent.hh"
00050 #include <cmath>
00051 using namespace std;
00052
00053
00054
00055
00056 G4QContent::G4QContent(G4int d, G4int u, G4int s_value, G4int ad, G4int au, G4int as):
00057 nD(d),nU(u),nS(s_value),nAD(ad),nAU(au),nAS(as)
00058 {
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076 }
00077
00078
00079 G4QContent::G4QContent(std::pair<G4int,G4int> PP): nD(0),nU(0),nS(0),nAD(0),nAU(0),nAS(0)
00080 {
00081 G4int P1=PP.first;
00082 G4int P2=PP.second;
00083 if(!P1 || !P2)
00084 {
00085
00086
00087 G4ExceptionDescription ed;
00088 ed << "Wrong parton pair: zero parton P1=" << P1 << ", P2=" << P2 << G4endl;
00089 G4Exception("G4QContent::G4QContent()", "HAD_CHPS_0072", FatalException, ed);
00090 }
00091 #ifdef debug
00092 G4cout<<"G4QContent::PairConstr: P1="<<P1<<", P2="<<P2<<G4endl;
00093 #endif
00094 G4bool suc=true;
00095 G4int A1=P1;
00096 if (P1 > 7) A1= P1/100;
00097 else if(P1 <-7) A1=(-P1)/100;
00098 else if(P1 < 0) A1= -P1;
00099 G4int P11=0;
00100 G4int P12=0;
00101 if(A1>7)
00102 {
00103 P11=A1/10;
00104 P12=A1%10;
00105 }
00106 if(P1>0)
00107 {
00108 if(!P11)
00109 {
00110 if (A1==1) ++nD;
00111 else if(A1==2) ++nU;
00112 else if(A1==3) ++nS;
00113 else suc=false;
00114 }
00115 else
00116 {
00117 if (P11==1) ++nD;
00118 else if(P11==2) ++nU;
00119 else if(P11==3) ++nS;
00120 else suc=false;
00121 if (P12==1) ++nD;
00122 else if(P12==2) ++nU;
00123 else if(P12==3) ++nS;
00124 else suc=false;
00125 }
00126 }
00127 else
00128 {
00129 if(!P11)
00130 {
00131 if (A1==1) ++nAD;
00132 else if(A1==2) ++nAU;
00133 else if(A1==3) ++nAS;
00134 else suc=false;
00135 }
00136 else
00137 {
00138 if (P11==1) ++nAD;
00139 else if(P11==2) ++nAU;
00140 else if(P11==3) ++nAS;
00141 else suc=false;
00142 if (P12==1) ++nAD;
00143 else if(P12==2) ++nAU;
00144 else if(P12==3) ++nAS;
00145 else suc=false;
00146 }
00147 }
00148 #ifdef debug
00149 G4cout<<"G4QContent::PCo:1:"<<nD<<","<<nU<<","<<nS<<","<<nAD<<","<<nAU<<","<<nAS<<G4endl;
00150 #endif
00151 G4int A2=P2;
00152 if (P2 > 7) A2= P2/100;
00153 else if(P2 <-7) A2=(-P2)/100;
00154 else if(P2 < 0) A2= -P2;
00155 G4int P21=0;
00156 G4int P22=0;
00157 if(A2>7)
00158 {
00159 P21=A2/10;
00160 P22=A2%10;
00161 }
00162 if(P2>0)
00163 {
00164 if(!P21)
00165 {
00166 if (A2==1) ++nD;
00167 else if(A2==2) ++nU;
00168 else if(A2==3) ++nS;
00169 else suc=false;
00170 }
00171 else
00172 {
00173 if (P21==1) ++nD;
00174 else if(P21==2) ++nU;
00175 else if(P21==3) ++nS;
00176 else suc=false;
00177 if (P22==1) ++nD;
00178 else if(P22==2) ++nU;
00179 else if(P22==3) ++nS;
00180 else suc=false;
00181 }
00182 }
00183 else
00184 {
00185 if(!P21)
00186 {
00187 if (A2==1) ++nAD;
00188 else if(A2==2) ++nAU;
00189 else if(A2==3) ++nAS;
00190 else suc=false;
00191 }
00192 else
00193 {
00194 if (P21==1) ++nAD;
00195 else if(P21==2) ++nAU;
00196 else if(P21==3) ++nAS;
00197 else suc=false;
00198 if (P22==1) ++nAD;
00199 else if(P22==2) ++nAU;
00200 else if(P22==3) ++nAS;
00201 else suc=false;
00202 }
00203 }
00204 if(!suc)
00205 {
00206
00207
00208 G4ExceptionDescription ed;
00209 ed << "Impossible parton pair, P1=" << P1 << ",P2=" << P2 << G4endl;
00210 G4Exception("G4QContent::G4QContent()", "HAD_CHPS_0073", FatalException, ed);
00211 }
00212 #ifdef debug
00213 G4cout<<"G4QContent::PCo:2:"<<nD<<","<<nU<<","<<nS<<","<<nAD<<","<<nAU<<","<<nAS<<G4endl;
00214 #endif
00215 }
00216
00217 G4QContent::G4QContent(const G4QContent& right)
00218 {
00219 nU = right.nU;
00220 nD = right.nD;
00221 nS = right.nS;
00222 nAU = right.nAU;
00223 nAD = right.nAD;
00224 nAS = right.nAS;
00225 }
00226
00227 G4QContent::G4QContent(G4QContent* right)
00228 {
00229 nU = right->nU;
00230 nD = right->nD;
00231 nS = right->nS;
00232 nAU = right->nAU;
00233 nAD = right->nAD;
00234 nAS = right->nAS;
00235 }
00236
00237
00238 const G4QContent& G4QContent::operator=(const G4QContent &right)
00239 {
00240 if(this != &right)
00241 {
00242 nU = right.nU;
00243 nD = right.nD;
00244 nS = right.nS;
00245 nAU = right.nAU;
00246 nAD = right.nAD;
00247 nAS = right.nAS;
00248 }
00249 return *this;
00250 }
00251
00252
00253 ostream& operator<<(ostream& lhs, G4QContent& rhs)
00254 {
00255 lhs << "{" << rhs.GetD() << "," << rhs.GetU() << "," << rhs.GetS() << ","
00256 << rhs.GetAD() << "," << rhs.GetAU() << "," << rhs.GetAS() << "}";
00257 return lhs;
00258 }
00259
00260
00261 ostream& operator<<(ostream& lhs, const G4QContent& rhs)
00262 {
00263 lhs << "{" << rhs.GetD() << "," << rhs.GetU() << "," << rhs.GetS() << ","
00264 << rhs.GetAD() << "," << rhs.GetAU() << "," << rhs.GetAS() << "}";
00265 return lhs;
00266 }
00267
00268
00269 G4QContent G4QContent::operator-=(const G4QContent& rhs)
00270 {
00271 #ifdef debug
00272 G4cout<<"G4QC::-=(const): is called:"<<G4endl;
00273 #endif
00274 G4int rD=rhs.nD;
00275 G4int rU=rhs.nU;
00276 G4int rS=rhs.nS;
00277 G4int rAD=rhs.nAD;
00278 G4int rAU=rhs.nAU;
00279 G4int rAS=rhs.nAS;
00284 if(nU<rU||nAU<rAU||nD<rD||nAD<rAD)
00285 {
00286 G4int dU=rU-nU;
00287 G4int dAU=rAU-nAU;
00288 if(dU>0||dAU>0)
00289 {
00290 G4int kU=dU;
00291 if(kU<dAU) kU=dAU;
00292 G4int mU=rU;
00293 if(rAU<mU) mU=rAU;
00294 if(kU<=mU)
00295 {
00296 rU-=kU;
00297 rAU-=kU;
00298 rD+=kU;
00299 rAD+=kU;
00300 }
00301 else
00302 {
00303 rU-=mU;
00304 rAU-=mU;
00305 rD+=mU;
00306 rAD+=mU;
00307 }
00308 }
00309 G4int dD=rD-nD;
00310 G4int dAD=rAD-nAD;
00311 if(dD>0||dAD>0)
00312 {
00313 G4int kD=dD;
00314 if(kD<dAD) kD=dAD;
00315 G4int mD=rD;
00316 if(rAD<mD) mD=rAD;
00317 if(kD<=mD)
00318 {
00319 rD-=kD;
00320 rAD-=kD;
00321 rU+=kD;
00322 rAU+=kD;
00323 }
00324 else
00325 {
00326 rD-=mD;
00327 rAD-=mD;
00328 rU+=mD;
00329 rAU+=mD;
00330 }
00331 }
00332 }
00333 #ifdef debug
00334 G4cout<<"G4QC::-=:comp: "<<rD<<","<<rU<<","<<rS<<","<<rAD<<","<<rAU<<","<<rAS<<G4endl;
00335 #endif
00336 if(rS==1 && rAS==1 && (nS<1 || nAS<1))
00337 {
00338 rS =0;
00339 rAS=0;
00340 if(nU>rU&&nAU>rAU)
00341 {
00342 rU +=1;
00343 rAU+=1;
00344 }
00345 else
00346 {
00347 rD +=1;
00348 rAD+=1;
00349 }
00350 }
00351 nD -= rD;
00352 if (nD<0)
00353 {
00354 nAD -= nD;
00355 nD = 0;
00356 }
00357 nU -= rU;
00358 if (nU<0)
00359 {
00360 nAU -= nU;
00361 nU = 0;
00362 }
00363 nS -= rS;
00364 if (nS<0)
00365 {
00366 nAS -= nS;
00367 nS = 0;
00368 }
00369 nAD -= rAD;
00370 if (nAD<0)
00371 {
00372 nD -= nAD;
00373 nAD = 0;
00374 }
00375 nAU -= rAU;
00376 if (nAU<0)
00377 {
00378 nU -= nAU;
00379 nAU = 0;
00380 }
00381 nAS -= rAS;
00382 if (nAS<0)
00383 {
00384 nS -= nAS;
00385 nAS = 0;
00386 }
00387 return *this;
00388 }
00389
00390
00391 G4QContent G4QContent::operator-=(G4QContent& rhs)
00392 {
00393 #ifdef debug
00394 G4cout<<"G4QC::-=: is called:"<<G4endl;
00395 #endif
00396 G4int rD=rhs.nD;
00397 G4int rU=rhs.nU;
00398 G4int rS=rhs.nS;
00399 G4int rAD=rhs.nAD;
00400 G4int rAU=rhs.nAU;
00401 G4int rAS=rhs.nAS;
00402 G4int rQ =rD+rU+rS;
00403 G4int rAQ=rAD+rAU+rAS;
00404 G4int nQ =nD+nU+nS;
00405 G4int nAQ=nAD+nAU+nAS;
00406 if(nQ<rQ||nAQ<rAQ)
00407 {
00408 G4int dU=rU-nU;
00409 G4int dAU=rAU-nAU;
00410 if(dU>0||dAU>0)
00411 {
00412 G4int kU=dU;
00413 if(kU<dAU) kU=dAU;
00414 G4int mU=rU;
00415 if(rAU<mU) mU=rAU;
00416 if(kU<=mU)
00417 {
00418 rU-=kU;
00419 rAU-=kU;
00420 rD+=kU;
00421 rAD+=kU;
00422 }
00423 else
00424 {
00425 rU-=mU;
00426 rAU-=mU;
00427 rD+=mU;
00428 rAD+=mU;
00429 }
00430 }
00431 G4int dD=rD-nD;
00432 G4int dAD=rAD-nAD;
00433 if(dD>0||dAD>0)
00434 {
00435 G4int kD=dD;
00436 if(kD<dAD) kD=dAD;
00437 G4int mD=rD;
00438 if(rAD<mD) mD=rAD;
00439 if(kD<=mD)
00440 {
00441 rD-=kD;
00442 rAD-=kD;
00443 rU+=kD;
00444 rAU+=kD;
00445 }
00446 else
00447 {
00448 rD-=mD;
00449 rAD-=mD;
00450 rU+=mD;
00451 rAU+=mD;
00452 }
00453 }
00454 }
00455 if(rS==1 && rAS==1 && (nS<1 || nAS<1))
00456 {
00457 rS =0;
00458 rAS=0;
00459 if(nU>rU&&nAU>rAU)
00460 {
00461 rU +=1;
00462 rAU+=1;
00463 }
00464 else
00465 {
00466 rD +=1;
00467 rAD+=1;
00468 }
00469 }
00470 nD -= rD;
00471 if (nD<0)
00472 {
00473 nAD -= nD;
00474 nD = 0;
00475 }
00476 nU -= rU;
00477 if (nU<0)
00478 {
00479 nAU -= nU;
00480 nU = 0;
00481 }
00482 nS -= rS;
00483 if (nS<0)
00484 {
00485 nAS -= nS;
00486 nS = 0;
00487 }
00488 nAD -= rAD;
00489 if (nAD<0)
00490 {
00491 nD -= nAD;
00492 nAD = 0;
00493 }
00494 nAU -= rAU;
00495 if (nAU<0)
00496 {
00497 nU -= nAU;
00498 nAU = 0;
00499 }
00500 nAS -= rAS;
00501 if (nAS<0)
00502 {
00503 nS -= nAS;
00504 nAS = 0;
00505 }
00506 return *this;
00507 }
00508
00509
00510 G4QContent operator+(const G4QContent& lhs, const G4QContent& rhs)
00511 {
00512 G4QContent s_value = lhs;
00513 return s_value += rhs;
00514 }
00515
00516
00517 G4QContent operator-(const G4QContent& lhs, const G4QContent& rhs)
00518 {
00519 G4QContent s_value = lhs;
00520 return s_value -= rhs;
00521 }
00522
00523
00524 G4QContent operator*(const G4QContent& lhs, const G4int& rhs)
00525 {
00526 G4QContent s_value = lhs;
00527 return s_value *= rhs;
00528 }
00529
00530
00531 G4QContent operator*(const G4int& lhs, const G4QContent& rhs)
00532 {
00533 G4QContent s_value = rhs;
00534 return s_value *= lhs;
00535 }
00536
00537
00538 G4QContent::~G4QContent() {}
00539
00540
00541 G4bool G4QContent::SubtractPi0()
00542 {
00543 #ifdef debug
00544 G4cout<<"G4QC::SubtractPi0: U="<<nU<<", AU="<<nAU<<", D="<<nD<<", AD="<<nAD<<G4endl;
00545 #endif
00546 G4int tot=GetTot();
00547 G4int ab =GetBaryonNumber();
00548 if(ab){if(tot<3*ab+2) return false;}
00549 else if(tot<4) return false;
00550
00551 if(nU>0 && nAU>0)
00552 {
00553 nU--;
00554 nAU--;
00555 return true;
00556 }
00557 else if(nD>0 && nAD>0)
00558 {
00559 nD--;
00560 nAD--;
00561 return true;
00562 }
00563 return false;
00564 }
00565
00566
00567 G4bool G4QContent::SubtractPion()
00568 {
00569 #ifdef debug
00570 G4cout<<"G4QC::SubtractPion: U="<<nU<<", AU="<<nAU<<", D="<<nD<<", AD="<<nAD<<G4endl;
00571 #endif
00572 G4int tot=GetTot();
00573 G4int ab =GetBaryonNumber();
00574 if(ab){if(tot<3*ab+2) return false;}
00575 else if(tot<4) return false;
00576
00577 if((nU>nAU || (ab && nU>0))&& nAD>0)
00578 {
00579 nU--;
00580 nAD--;
00581 return true;
00582 }
00583 else if((nAU>nU || (ab && nAU>0)) && nD>0)
00584 {
00585 nAU--;
00586 nD--;
00587 return true;
00588 }
00589 return false;
00590 }
00591
00592
00593 G4bool G4QContent::SubtractHadron(G4QContent h)
00594 {
00595 #ifdef debug
00596 G4cout<<"G4QC::SubtractHadron "<<h<<" is called for QC="<<GetThis()<<G4endl;
00597 #endif
00598 if(h.GetU()<=nU && h.GetD()<=nD && h.GetS()<=nS&&
00599 h.GetAU()<=nAU&& h.GetAD()<=nAD&& h.GetAS()<=nAS) return true;
00600 return false;
00601 }
00602
00603
00604 G4bool G4QContent::SubtractKaon(G4double mQ)
00605 {
00606 #ifdef debug
00607 G4cout<<"G4QC::SubtractKaon is called: QC="<<GetThis()<<G4endl;
00608 #endif
00609 if(mQ<640.) return false;
00610 G4int tot=GetTot();
00611 G4int ab =GetBaryonNumber();
00612 if(ab){if(tot<3*ab+2) return false;}
00613 else if(tot<4) return false;
00614
00615 if((nS>nAS || (ab && nS>0)) && (nAD>0 || nAU>0))
00616 {
00617 nS--;
00618 if (nAU>0) nAU--;
00619 else nAD--;
00620 return true;
00621 }
00622 else if((nAS>nS || (ab && nAS>0)) && (nD>0 || nU>0))
00623 {
00624 nAS--;
00625 if (nU>0) nU--;
00626 else nD--;
00627 return true;
00628 }
00629 #ifdef debug
00630 G4cout<<"G4QCont::SubtractKaon Can't SubtractKaon: QC="<<GetThis()<<G4endl;
00631 #endif
00632 return false;
00633 }
00634
00635
00636 G4QContent G4QContent::SplitChipo (G4double mQ)
00637 {
00638 G4QContent Pi(0,1,0,1,0,0);
00639 if (nU>0&&nAU>0) Pi=G4QContent(0,1,0,0,1,0);
00640 else if (nD>0&&nAD>0) Pi=G4QContent(1,0,0,1,0,0);
00641 else if (nD>=nU&&nAU>=nAD) Pi=G4QContent(1,0,0,0,1,0);
00642 G4int bn=GetBaryonNumber();
00643 G4int b =abs(bn);
00644 if(!b && mQ<545.&&nS>0&&nAS>0)
00645 {
00646 G4int ss= nS;
00647 if(nAS<nS) ss=nAS;
00648 nS -=ss;
00649 nAS-=ss;
00650 }
00651 if (!b)DecQAQ(-4);
00652 else if(b==1)DecQAQ(-5);
00653 else DecQAQ(0);
00654 G4int tot=GetTot();
00655 G4int q=GetQ();
00656 G4int aq=GetAQ();
00657 G4QContent r=Pi;
00658 if((tot!=4||q!=2) && (tot!=5||(q!=1&&aq!=1)) && (tot!=6||abs(b)!=2))
00659 {
00660
00661 G4cerr<<"***G4QCont::SplitChipo: QC="<<GetThis()<<G4endl;
00662
00663 }
00664 else if(tot==4)
00665 {
00666 r=GetThis();
00667 if (r.SubtractPi0()) ;
00668 else if(r.SubtractPion()) ;
00669 else if(r.SubtractKaon(mQ)) ;
00670 else
00671 {
00672
00673 G4cerr<<"***G4QCont::SplitChipo:MesTot="<<tot<<",b="<<b<<",q="<<q<<",a="<<aq<<G4endl;
00674
00675 }
00676 }
00677 else if(b==1&&tot==5)
00678 {
00679 if(nU==3)
00680 {
00681 r.SetU(1);
00682 r+=IndAQ();
00683 }
00684 else if(nD==3)
00685 {
00686 r.SetD(1);
00687 r+=IndAQ();
00688 }
00689 else if(nS==3)
00690 {
00691 r.SetS(1);
00692 r+=IndAQ();
00693 }
00694 else if(nAU==3)
00695 {
00696 r.SetAU(1);
00697 r+=IndQ();
00698 }
00699 else if(nAD==3)
00700 {
00701 r.SetAD(1);
00702 r+=IndQ();
00703 }
00704 else if(nAS==3)
00705 {
00706 r.SetAS(1);
00707 r+=IndQ();
00708 }
00709 else if(q==1&&nU)
00710 {
00711 r.SetU(1);
00712 if(nAU) r.SetAU(1);
00713 else r.SetAD(1);
00714 }
00715 else if(q==1&&nD)
00716 {
00717 r.SetD(1);
00718 if(nAD) r.SetAD(1);
00719 else r.SetAU(1);
00720 }
00721 else if(q==1&&nS)
00722 {
00723 r.SetS(1);
00724 if(nAS) r.SetAS(1);
00725 else r.SetAU(1);
00726 }
00727 else if(aq==1&&nAU)
00728 {
00729 r.SetAU(1);
00730 if(nU) r.SetU(1);
00731 else r.SetD(1);
00732 }
00733 else if(aq==1&&nAD)
00734 {
00735 r.SetAD(1);
00736 if(nD) r.SetD(1);
00737 else r.SetU(1);
00738 }
00739 else if(aq==1&&nAS)
00740 {
00741 r.SetAS(1);
00742 if(nS) r.SetS(1);
00743 else r.SetU(1);
00744 }
00745 else
00746 {
00747
00748 G4cerr<<"***G4QCont::SplitChipo: Baryonic tot=5,b=1,qCont="<<GetThis()<<G4endl;
00749
00750 }
00751 }
00752 else if(tot==b*3)
00753 {
00754 r=GetThis();
00755 if (bn>0)
00756 {
00757 G4QContent la(1,1,1,0,0,0);
00758 G4QContent nt(2,1,0,0,0,0);
00759 G4QContent pr(1,2,0,0,0,0);
00760 G4QContent ks(0,1,2,0,0,0);
00761 if (nD>nU) ks=G4QContent(1,0,2,0,0,0);
00762 G4QContent dm(3,0,0,0,0,0);
00763 G4QContent dp(0,3,0,0,0,0);
00764 G4QContent om(0,0,3,0,0,0);
00765 if (nU>=nD&&nU>=nS)
00766 {
00767 if (r.SubtractHadron(pr)) r-=pr;
00768 else if(r.SubtractHadron(dp)) r-=dp;
00769 else if(r.SubtractHadron(nt)) r-=nt;
00770 else if(r.SubtractHadron(la)) r-=la;
00771 else if(r.SubtractHadron(dm)) r-=dm;
00772 else
00773 {
00774
00775 G4cerr<<"***G4QCont::SplitChipo:Dibar (1) tot=6, b=2, qCont="<<GetThis()<<G4endl;
00776
00777 }
00778 }
00779 else if(nD>=nU&&nD>=nS)
00780 {
00781 if (r.SubtractHadron(nt)) r-=nt;
00782 else if(r.SubtractHadron(dm)) r-=dm;
00783 else if(r.SubtractHadron(pr)) r-=pr;
00784 else if(r.SubtractHadron(dp)) r-=dp;
00785 else if(r.SubtractHadron(la)) r-=la;
00786 else
00787 {
00788
00789 G4cerr<<"***G4QContent::SplitChipo:Dib(2) tot=6, b=2, qCont="<<GetThis()<<G4endl;
00790
00791 }
00792 }
00793 else
00794 {
00795 if (r.SubtractHadron(la)) r-=la;
00796 else if(r.SubtractHadron(ks)) r-=ks;
00797 else if(r.SubtractHadron(om)) r-=om;
00798 else if(r.SubtractHadron(pr)) r-=pr;
00799 else if(r.SubtractHadron(nt)) r-=nt;
00800 else
00801 {
00802
00803 G4cerr<<"***G4QContent::SplitChipo:Dib(3) tot=6, b=2, qCont="<<GetThis()<<G4endl;
00804
00805 }
00806 }
00807 }
00808 else
00809 {
00810 G4QContent la(0,0,0,1,1,1);
00811 G4QContent pr(0,0,0,1,2,0);
00812 G4QContent nt(0,0,0,2,1,0);
00813 G4QContent ks(0,1,2,0,0,0);
00814 if(nAD>nAU)ks=G4QContent(0,0,0,1,0,2);
00815 G4QContent dm(0,0,0,3,0,0);
00816 G4QContent dp(0,0,0,0,3,0);
00817 G4QContent om(0,0,0,0,0,3);
00818 if (nAU>=nAD&&nAU>=nAS)
00819 {
00820 if (r.SubtractHadron(pr)) r-=pr;
00821 else if(r.SubtractHadron(dp)) r-=dp;
00822 else if(r.SubtractHadron(nt)) r-=nt;
00823 else if(r.SubtractHadron(la)) r-=la;
00824 else if(r.SubtractHadron(dm)) r-=dm;
00825 else
00826 {
00827
00828 G4cerr<<"***G4QContent::SplitChipo:ADib(1) tot=6,b=2, qCont="<<GetThis()<<G4endl;
00829
00830 }
00831 }
00832 else if(nAD>=nAU&&nAD>=nAS)
00833 {
00834 if (r.SubtractHadron(nt)) r-=nt;
00835 else if(r.SubtractHadron(dm)) r-=dm;
00836 else if(r.SubtractHadron(pr)) r-=pr;
00837 else if(r.SubtractHadron(dp)) r-=dp;
00838 else if(r.SubtractHadron(la)) r-=la;
00839 else
00840 {
00841
00842 G4cerr<<"***G4QContent::SplitChipo:ADib(2) tot=6,b=2, qCont="<<GetThis()<<G4endl;
00843
00844 }
00845 }
00846 else
00847 {
00848 if (r.SubtractHadron(la)) r-=la;
00849 else if(r.SubtractHadron(ks)) r-=ks;
00850 else if(r.SubtractHadron(om)) r-=om;
00851 else if(r.SubtractHadron(pr)) r-=pr;
00852 else if(r.SubtractHadron(nt)) r-=nt;
00853 else
00854 {
00855
00856 G4cerr<<"***G4QContent::SplitChipo:ADib(3) tot=6,b=2, qCont="<<GetThis()<<G4endl;
00857
00858 }
00859 }
00860 }
00861 }
00862 else
00863 {
00864
00865 G4cerr<<"*G4QContent::SplitChipolino:UnknownHadron with QuarkCont="<<GetThis()<<G4endl;
00866
00867 }
00868 return r;
00869 }
00870
00871
00872 G4QContent G4QContent::IndQ (G4int index)
00873 {
00874 #ifdef debug
00875 G4cout << "G4QC::IndQ is called"<<G4endl;
00876 #endif
00877 if(index<nD) return G4QContent(1,0,0,0,0,0);
00878 else if(index<nD+nU) return G4QContent(0,1,0,0,0,0);
00879 else if(index<nD+nU+nS) return G4QContent(0,0,1,0,0,0);
00880
00881 else G4cerr<<"***G4QC::IndQ:index="<<index<<" for the QuarkContent="<<GetThis()<<G4endl;
00882
00883
00884 return G4QContent(0,0,0,0,0,0);
00885 }
00886
00887
00888 G4QContent G4QContent::IndAQ (G4int index)
00889 {
00890 #ifdef debug
00891 G4cout << "G4QC::IndAQ is called"<<G4endl;
00892 #endif
00893 if(index<nAD) return G4QContent(0,0,0,1,0,0);
00894 else if(index<nAD+nAU) return G4QContent(0,0,0,0,1,0);
00895 else if(index<nAD+nAU+nAS) return G4QContent(0,0,0,0,0,1);
00896
00897 else G4cerr<<"***G4QC::IndAQ:index="<<index<<" for the QuarkContent="<<GetThis()<<G4endl;
00898
00899
00900 return G4QContent(0,0,0,0,0,0);
00901 }
00902
00903
00904 G4int G4QContent::DecQAQ(const G4int& nQAQ)
00905 {
00906 #ifdef debug
00907 G4cout<<"G4QCont::DecQC: n="<<nQAQ<<","<<GetThis()<<G4endl;
00908 #endif
00909 G4int ban = GetBaryonNumber();
00910 G4int tot = GetTot();
00911 if (tot==ban*3) return 0;
00912 G4int nUP=0;
00913 if (nU>=nAU) nUP=nAU;
00914 else nUP= nU;
00915
00916 G4int nDP=0;
00917 if (nD>=nAD) nDP=nAD;
00918 else nDP= nD;
00919
00920 G4int nSP=0;
00921 if (nS>=nAS) nSP=nAS;
00922 else nSP= nS;
00923
00924 G4int nLP =nUP+nDP;
00925 G4int nTotP=nLP+nSP;
00926 G4int nReal=nQAQ;
00927 G4int nRet =nTotP-nQAQ;
00928 if (nQAQ<0)
00929 {
00930 G4int res=tot+nQAQ;
00931 #ifdef debug
00932 G4cout<<"G4QC::DecQC: tot="<<tot<<", nTP="<<nTotP<<", res="<<res<<G4endl;
00933 #endif
00934 if(res<0)
00935 {
00936 IncQAQ(1,0.);
00937 return 1;
00938 }
00939 res -=nTotP+nTotP;
00940 if(res<0) nReal=nTotP+res/2;
00941 else nReal=nTotP;
00942 nRet = tot/2-nReal;
00943 }
00944 else if(!nQAQ)
00945 {
00946 nReal=nTotP;
00947 nRet =0;
00948 }
00949 else if(nRet<0) nReal=nTotP;
00950
00951 if (!nReal) return nRet;
00952
00953 #ifdef debug
00954 G4cout<<"G4QC::DecQC: demanded "<<nQAQ<<" pairs, executed "<<nReal<<" pairs"<<G4endl;
00955 #endif
00957 for (G4int i=0; i<nReal; i++)
00958 {
00959 G4double base = nTotP;
00960
00961 G4int j = static_cast<int>(base*G4UniformRand());
00962 if (nUP && j<nUP && (nRet>2 || nUP>1 || (nD<2 && nS<2)))
00963 {
00964 #ifdef debug
00965 G4cout<<"G4QC::DecQC: decrementing UAU pair UP="<<nUP<<", QC="<<GetThis()<<G4endl;
00966 #endif
00967 nU--;
00968 nAU--;
00969 nUP--;
00970 nLP--;
00971 nTotP--;
00972 }
00973 else if (nDP && j<nLP && (nRet>2 || nDP>1 || (nU<2 && nS<2)))
00974 {
00975 #ifdef debug
00976 G4cout<<"G4QC::DecQC: decrementing DAD pair DP="<<nDP<<", QC="<<GetThis()<<G4endl;
00977 #endif
00978 nD--;
00979 nAD--;
00980 nDP--;
00981 nLP--;
00982 nTotP--;
00983 }
00984 else if (nSP&& (nRet>2 || nSP>1 || (nU<2 && nD<2)))
00985 {
00986 #ifdef debug
00987 G4cout<<"G4QC::DecQC: decrementing SAS pair SP="<<nSP<<", QC="<<GetThis()<<G4endl;
00988 #endif
00989 nS--;
00990 nAS--;
00991 nSP--;
00992 nTotP--;
00993 }
00994 else if (nUP)
00995 {
00996 #ifdef debug
00997 G4cout<<"G4QC::DecQC:Decrement UAU pair (final) UP="<<nUP<<",QC="<<GetThis()<<G4endl;
00998 #endif
00999 nU--;
01000 nAU--;
01001 nUP--;
01002 nLP--;
01003 nTotP--;
01004 }
01005 else if (nDP)
01006 {
01007 #ifdef debug
01008 G4cout<<"G4QC::DecQC:Decrement DAD pair (final) DP="<<nDP<<",QC="<<GetThis()<<G4endl;
01009 #endif
01010 nD--;
01011 nAD--;
01012 nDP--;
01013 nLP--;
01014 nTotP--;
01015 }
01016 else if (nSP)
01017 {
01018 #ifdef debug
01019 G4cout<<"G4QC::DecQC: decrementing SAS pair SP="<<nSP<<", QC="<<GetThis()<<G4endl;
01020 #endif
01021 nS--;
01022 nAS--;
01023 nSP--;
01024 nTotP--;
01025 }
01026 else G4cout<<"***G4QC::DecQC:i="<<i<<",j="<<j<<",D="<<nDP<<",U="<<nUP<<",S="<<nSP
01027 <<",T="<<nTotP<<",nRet="<<nRet<<", QC="<<GetThis()<<G4endl;
01028 }
01029 #ifdef debug
01030 G4cout<<"G4QC::DecQC: >->-> OUT <-<-< nRet="<<nRet<<", QC="<<GetThis()<<G4endl;
01031 #endif
01032 return nRet;
01033 }
01034
01035
01036 void G4QContent::IncQAQ(const G4int& nQAQ, const G4double& sProb)
01037 {
01038 G4int tot = GetTot();
01039 G4QContent mQC = GetThis();
01040 for (int i=0; i<nQAQ; i++)
01041 {
01042 G4int j = static_cast<int>((2.+sProb)*G4UniformRand());
01043 #ifdef debug
01044 G4cout<<"IncQC:out QC="<<GetThis()<<",j="<<j<<" for i="<<i<<G4endl;
01045 #endif
01046
01047 if ( !j && (nU<=nD || nU<=nS))
01048 {
01049 nU++;
01050 nAU++;
01051 tot+=2;
01052 }
01053
01054 else if (j==1 && (nD<=nU || nD<=nS))
01055 {
01056 nD++;
01057 nAD++;
01058 tot+=2;
01059 }
01060
01061 else if (j>1&& (nS<=nU || nS<=nD))
01062 {
01063 nS++;
01064 nAS++;
01065 tot+=2;
01066 }
01067 else if (!j)
01068 {
01069 nD++;
01070 nAD++;
01071 tot+=2;
01072 }
01073 else if (j==1)
01074 {
01075 nU++;
01076 nAU++;
01077 tot+=2;
01078 }
01079 else
01080 {
01081 nS++;
01082 nAS++;
01083 tot+=2;
01084 }
01085
01086
01087
01088
01089
01090
01091
01092
01093
01094
01095
01096
01097 }
01098 }
01099
01100
01101 G4int G4QContent::GetP() const
01102 {
01103 G4int rD=nD-nAD;
01104 G4int rU=nU-nAU;
01105 G4int rS=nS-nAS;
01106 G4int dQ=rD-rU;
01107 G4int b3=rD+rU+rS;
01108 return (b3-3*(rS+dQ))/6;
01109 }
01110
01111
01112 G4int G4QContent::GetN() const
01113 {
01114 G4int rD=nD-nAD;
01115 G4int rU=nU-nAU;
01116 G4int rS=nS-nAS;
01117 G4int dQ=rD-rU;
01118 G4int b3=rD+rU+rS;
01119 return (b3-3*(rS-dQ))/6;
01120 }
01121
01122
01123 G4int G4QContent::GetL() const
01124 {
01125 G4int rS=nS-nAS;
01126 return rS;
01127 }
01128
01129
01130 G4int G4QContent::GetAP() const
01131 {
01132 G4int rD=nAD-nD;
01133 G4int rU=nAU-nU;
01134 G4int rS=nAS-nS;
01135 G4int dQ=rD-rU;
01136 G4int b3=rD+rU+rS;
01137 return (b3-3*(rS+dQ))/6;
01138 }
01139
01140
01141 G4int G4QContent::GetAN() const
01142 {
01143 G4int rD=nAD-nD;
01144 G4int rU=nAU-nU;
01145 G4int rS=nAS-nS;
01146 G4int dQ=rD-rU;
01147 G4int b3=rD+rU+rS;
01148 return (b3-3*(rS-dQ))/6;
01149 }
01150
01151
01152 G4int G4QContent::GetAL() const
01153 {
01154 G4int rS=nAS-nS;
01155 return rS;
01156 }
01157
01158
01159 G4int G4QContent::GetCharge() const
01160 {
01161 static const G4int cU = 2;
01162 static const G4int cD =-1;
01163 static const G4int cS =-1;
01164 static const G4int cAU =-2;
01165 static const G4int cAD = 1;
01166 static const G4int cAS = 1;
01167
01168 G4int c=0;
01169 if(nU) c+=nU*cU;
01170 if(nD) c+=nD*cD;
01171 if(nS) c+=nS*cS;
01172 if(nAU)c+=nAU*cAU;
01173 if(nAD)c+=nAD*cAD;
01174 if(nAS)c+=nAS*cAS;
01175
01176 if(c%3) G4cerr<<"***G4QCont:GetCharge:c="<<c<<"/3 isn't integer, QC="<<GetThis()<<G4endl;
01177
01178 return c/3;
01179 }
01180
01181
01182 G4int G4QContent::GetBaryonNumber() const
01183 {
01184 #ifdef pdebug
01185 G4cout<<"G4QContent::GetBarNum: U="<<nU<<", D="<<nD<<", S="<<nS<<", AU="<<nAU<<", AD="
01186 <<nAD<<", AS="<<nAS<<G4endl;
01187 #endif
01188 G4int b=nU+nD+nS-nAU-nAD-nAS;
01189
01190 if(b%3)
01191 {
01192
01193
01194 G4ExceptionDescription ed;
01195 ed << "Wrong Baryon Number: warning " << b << "/3 isn't an integer"
01196 << G4endl;
01197 G4Exception("G4QContent::GetBaryonNumber()", "HAD_CHPS_0072", FatalException, ed);
01198 }
01199
01200 return b/3;
01201 }
01202
01203
01204 G4int G4QContent::GetSPDGCode() const
01205 {
01206 G4int p = 0;
01207 G4int n = GetTot();
01208 if(!n) return 22;
01209 G4int mD=nD;
01210 if (nD<=0) mD=nAD;
01211 G4int mU=nU;
01212 if (nU<=0) mU=nAU;
01213 G4int mS=nS;
01214 if (nS<=0) mS= nAS;
01215
01216 if ( nU>nAU && nAU>0)
01217 {
01218 mU=nU-nAU;
01219 n-=nAU+nAU;
01220 }
01221 if (nAU>nU && nU>0)
01222 {
01223 mU=nAU-nU;
01224 n-=nU+nU;
01225 }
01226 if ( nD>nAD && nAD>0)
01227 {
01228 mD=nD-nAD;
01229 n-=nAD+nAD;
01230 }
01231 if (nAD>nD && nD>0)
01232 {
01233 mD=nAD-nD;
01234 n-=nD+nD;
01235 }
01236 if ( nS>nAS && nAS>0)
01237 {
01238 mS=nS-nAS;
01239 n-=nAS+nAS;
01240 }
01241 if (nAS>nS && nS>0)
01242 {
01243 mS= nAS-nS;
01244 n-=nS+nS;
01245 }
01246
01247 if (nAD==nD && nD>0)
01248 {
01249 G4int dD=nD+nD;
01250 if(n>dD)
01251 {
01252 mD=0;
01253 n-=dD;
01254 }
01255 else if (n==dD)
01256 {
01257 mD=2;
01258 n=2;
01259 }
01260 else
01261 {
01262 #ifdef debug
01263 G4cout<<"***G4QC::SPDG:CanD U="<<mU<<",D="<<mD<<",S="<<mS<<",QC="<<GetThis()<<G4endl;
01264 #endif
01265 return 0;
01266 }
01267 }
01268 if (nAU==nU && nU>0)
01269 {
01270 G4int dU=nU+nU;
01271 if(n>dU)
01272 {
01273 mU=0;
01274 n-=dU;
01275 }
01276 else if (n==dU)
01277 {
01278 mU=2;
01279 n=2;
01280 }
01281 else
01282 {
01283 #ifdef debug
01284 G4cout<<"***G4QC::SPDG:CanU U="<<mU<<",D="<<mD<<",S="<<mS<<",QC="<<GetThis()<<G4endl;
01285 #endif
01286 return 0;
01287 }
01288 }
01289 if (nAS==nS && nS>0)
01290 {
01291 G4int dS=nS+nS;
01292 if(n>dS)
01293 {
01294 mS=0;
01295 n-=dS;
01296 }
01297 else if (n==dS)
01298 {
01299 mS=2;
01300 n=2;
01301 }
01302 else
01303 {
01304 #ifdef debug
01305 G4cout<<"***G4QC::SPDG:CanS U="<<mU<<",D="<<mD<<",S="<<mS<<",QC="<<GetThis()<<G4endl;
01306 #endif
01307 return 0;
01308 }
01309 }
01310
01311 G4int b=GetBaryonNumber();
01312 G4int c=GetCharge();
01313 G4int s_value=GetStrangeness();
01314 #ifdef debug
01315 G4cout<<"G4QC::SPDGC:bef. b="<<b<<",n="<<n<<",c="<<c<<",s="<<s_value<<",Q="<<GetThis()<<G4endl;
01316 #endif
01317 if (b)
01318 {
01319
01320 G4int ab=abs(b);
01321 if(ab>=2 && n>=6)
01322 {
01323 G4int mI=nU-nAU-nD+nAD;
01324
01325
01326 if ( (b > 0 && s_value < -1) || (b < 0 && s_value > 1) ) return 10;
01327 else if (abs(mI) > 2 || mS > 2
01328 || (b > 0 && s_value < 0)
01329 || (b < 0 && s_value > 0)) return GetZNSPDGCode();
01330 else if(mU>=mS&&mD>=mS&&mU+mD+mS==3*b)
01331 {
01332 G4int mZ=(mU+mD-mS-mS+3*mI)/6;
01333 p = 90000000+1000*(1000*mS+mZ)+mZ-mI;
01334 if(b>0) return p;
01335 else return -p;
01336 }
01337 else return 10;
01338 }
01339
01340 if(n>5) return GetZNSPDGCode();
01341 if(n==5) return 10;
01342 if(mS>0)
01343 {
01344 p=3002;
01345 if (mS==3) p+=332;
01346 else if (mS==2)
01347 {
01348 if (mU==1 && mD==0) p+=320;
01349 else if (mU==0 && mD==1) p+=310;
01350 else
01351 {
01352 #ifdef debug
01353 G4cout<<"**G4QC::SPDG:ExoticBSS,U="<<mU<<",D="<<mD<<",S="<<mS<<GetThis()<<G4endl;
01354 #endif
01355 return GetZNSPDGCode();
01356 }
01357 }
01358 else if (mS==1)
01359 {
01360 if (mU==2 && mD==0) p+=220;
01361 else if (mU==1 && mD==1) p+=120;
01362 else if (mU==0 && mD==2) p+=110;
01363 else
01364 {
01365 #ifdef debug
01366 G4cout<<"***G4QC::SPDG:ExoticBS,U="<<mU<<",D="<<mD<<",S="<<mS<<GetThis()<<G4endl;
01367 #endif
01368 return GetZNSPDGCode();
01369 }
01370 }
01371 else
01372 {
01373 #ifdef debug
01374 G4cout<<"***G4QC::GetSPDG:ExoBarS,U="<<mU<<",D="<<mD<<",S="<<mS<<GetThis()<<G4endl;
01375 #endif
01376 return GetZNSPDGCode();
01377 }
01378 }
01379 else if (mU>0)
01380 {
01381 p=2002;
01382 if (mU==3 && mD==0) p+=222;
01383 else if (mU==2 && mD==1) p+=210;
01384 else if (mU==1 && mD==2) p+=110;
01385 else
01386 {
01387 #ifdef debug
01388 G4cout<<"***G4QC::SPDG:ExoBaryonU,U="<<mU<<",D="<<mD<<",S="<<mS<<GetThis()<<G4endl;
01389 #endif
01390 return GetZNSPDGCode();
01391 }
01392 }
01393 else if (mD==3) p=1114;
01394 else
01395 {
01396 #ifdef debug
01397 G4cout<<"**G4QC::SPDG:ExoticBaryonD,U="<<mU<<",D="<<mD<<",S="<<mS<<GetThis()<<G4endl;
01398 #endif
01399 return GetZNSPDGCode();
01400 }
01401 if (b<0) p=-p;
01402 }
01403 else
01404 {
01405 #ifdef debug
01406 G4cout<<"G4QC::SPDG:mDUS="<<mD<<","<<mU<<","<<mS<<",b,c,s="<<b<<","<<c<<","<<s_value<<G4endl;
01407 #endif
01408 if(n>4)
01409 {
01410 #ifdef debug
01411 G4cout<<"G4QC::SPDG:n>4 SEx:U="<<mU<<",D="<<mD<<",S="<<mS<<",QC="<<GetThis()<<G4endl;
01412 #endif
01413 return 0;
01414 }
01415 if(n==4) return 10;
01416 if(abs(s_value)>1)
01417 {
01418 #ifdef debug
01419 G4cout<<"**G4QC::SPDG:Stran="<<s_value<<",QC="<<GetThis()<<" - Superstrange Meson"<<G4endl;
01420 #endif
01421 return 0;
01422 }
01423
01424 if(mS>0)
01425 {
01426 p=301;
01427 if (mS==2)
01428 {
01429
01430
01431 p=221;
01432 }
01433 else if (mU==1 && mD==0) p+=20;
01434 else if (mU==0 && mD==1) p+=10;
01435 else
01436 {
01437 #ifdef debug
01438 G4cout<<"*G4QC::SPDG:ExMS U="<<mU<<",D="<<mD<<",S="<<mS<<",QC="<<GetThis()<<G4endl;
01439 #endif
01440 return 0;
01441 }
01442 }
01443 else if (mU>0)
01444 {
01445 p=201;
01446
01447 if (mU==2 && mD==0) p=111;
01448 else if (mU==1 && mD==1) p+=10;
01449 else
01450 {
01451 #ifdef debug
01452 G4cout<<"*G4QC::SPDG:ExMU U="<<mU<<",D="<<mD<<",S="<<mS<<",QC="<<GetThis()<<G4endl;
01453 #endif
01454 return 0;
01455 }
01456 }
01457 else if (mD==2) p=111;
01458 else
01459 {
01460 #ifdef debug
01461 G4cout<<"***G4QC::SPDG:ExMD U="<<mU<<",D="<<mD<<",S="<<mS<<",QC="<<GetThis()<<G4endl;
01462 #endif
01463 return 0;
01464 }
01465 if (c<0 || (c==0 && mS>0 && s_value>0)) p=-p;
01466 }
01467 #ifdef debug
01468 G4cout<<"G4QC::GetSPDG:output SPDGcode="<<p<<" for the QuarkContent="<<GetThis()<<G4endl;
01469 #endif
01470 return p;
01471 }
01472
01473
01474 G4int G4QContent::NOfCombinations(const G4QContent& rhs) const
01475 {
01476 G4int c=1;
01477 #ifdef ppdebug
01478 G4cout<<"G4QContent::NOfComb: This="<<GetThis()<<", selectQC="<<rhs<<G4endl;
01479 #endif
01480 G4int mD=rhs.GetD();
01481 G4int mU=rhs.GetU();
01482 G4int mS=rhs.GetS();
01483 G4int mAD=rhs.GetAD();
01484 G4int mAU=rhs.GetAU();
01485 G4int mAS=rhs.GetAS();
01486 G4int mN=mD+mU+mS-mAD-mAU-mAS;
01488 if (( ((nD < mD || nAD < mAD) && !(mD-mAD)) ||
01489 ((nU < mU || nAU < mAU) && !(mU-mAU)) ||
01490 ((nS < mS || nAS < mAS) && !(mS-mAS)) ) && !mN) return 1;
01491 if(mD>0)
01492 {
01493 int j=nD;
01494 if (j<=0) return 0;
01495 if(mD>1||j>1) for (int i=1; i<=mD; i++)
01496 {
01497 if(!j) return 0;
01498 c*=j/i;
01499 j--;
01500 }
01501 }
01502 if(mU>0)
01503 {
01504 int j=nU;
01505 if (j<=0) return 0;
01506 if(mU>1||j>1) for (int i=1; i<=mU; i++)
01507 {
01508 if(!j) return 0;
01509 c*=j/i;
01510 j--;
01511 }
01512 }
01513 if(mS>0)
01514 {
01515 int j=nS;
01516 if (j<=0) return 0;
01517 if(mS>1||j>1) for (int i=1; i<=mS; i++)
01518 {
01519 if(!j) return 0;
01520 c*=j/i;
01521 j--;
01522 }
01523 }
01524 if(mAD>0)
01525 {
01526 int j=nAD;
01527 if (j<=0) return 0;
01528 if(mAD>1||j>1) for (int i=1; i<=mAD; i++)
01529 {
01530 if(!j) return 0;
01531 c*=j/i;
01532 j--;
01533 }
01534 }
01535 if(mAU>0)
01536 {
01537 int j=nAU;
01538 if (j<=0) return 0;
01539 if(mAU>1||j>1) for (int i=1; i<=mAU; i++)
01540 {
01541 if(!j) return 0;
01542 c*=j/i;
01543 j--;
01544 }
01545 }
01546 if(mAS>0)
01547 {
01548 int j=nAS;
01549 if (j<=0) return 0;
01550 if(mAS>1||j>1) for (int i=1; i<=mAS; i++)
01551 {
01552 if(!j) return 0;
01553 c*=j/i;
01554 j--;
01555 }
01556 }
01557 return c;
01558 }
01559
01560
01561 std::pair<G4int,G4int> G4QContent::MakePartonPair() const
01562 {
01563 G4double S=0.;
01564 S+=nD;
01565 G4double dP=S;
01566 S+=nU;
01567 G4double uP=S;
01568 S+=nS;
01569 G4double sP=S;
01570 S+=nAD;
01571 G4double dA=S;
01572 S+=nAU;
01573 G4double uA=S;
01574 S+=nAS;
01575 if(!S)
01576 {
01577 G4int f= static_cast<int>(1.+2.3*G4UniformRand());
01578 return std::make_pair(f,-f);
01579 }
01580 G4int f=0;
01581 G4double R=S*G4UniformRand();
01582 if (R<dP) f=1;
01583 else if(R<uP) f=2;
01584 else if(R<sP) f=3;
01585 else if(R<dA) f=-1;
01586 else if(R<uA) f=-2;
01587 else f=-3;
01588
01589 if (f < 0) {
01590 if(nD || nU || nS)
01591 {
01592 if (nD) return std::make_pair(1,f);
01593 else if(nU) return std::make_pair(2,f);
01594 else return std::make_pair(3,f);
01595 }
01596 else
01597 {
01598
01599 G4int AD=nAD;
01600 if(f==-1) AD--;
01601 G4int AU=nAU;
01602 if(f==-1) AU--;
01603 G4int AS=nAS;
01604 if(f==-1) AS--;
01605 if (AS)
01606 {
01607 if (AS==2) return std::make_pair(-3303,f);
01608 else if(AU) return std::make_pair(-3201,f);
01609 else return std::make_pair(-3101,f);
01610 }
01611 else if(AU)
01612 {
01613 if (AU==2) return std::make_pair(-2203,f);
01614 else return std::make_pair(-2101,f);
01615 }
01616 else return std::make_pair(-1103,f);
01617 }
01618
01619 } else {
01620
01621 if(nAD || nAU || nAS)
01622 {
01623 if (nAD) return std::make_pair(f,-1);
01624 else if(nAU) return std::make_pair(f,-2);
01625 else return std::make_pair(f,-3);
01626 }
01627 else
01628 {
01629
01630
01631
01632 G4int AU=nU;
01633
01634 G4int AS=nS;
01635
01636
01637 if (AS) {
01638 if (AS==2) return std::make_pair(f,3303);
01639 else if(AU) return std::make_pair(f,3201);
01640 else return std::make_pair(f,3101);
01641
01642 } else if(AU) {
01643 if (AU==2) return std::make_pair(f,2203);
01644 else return std::make_pair(f,2101);
01645 }
01646 else return std::make_pair(f,1103);
01647 }
01648 }
01649 }
01650
01651
01652
01653
01654
01655 G4int G4QContent::AddParton(G4int pPDG) const
01656 {
01657 #ifdef debug
01658 G4cout<<"G4QContent::AddParton: This="<<GetThis()<<", pPDG="<<pPDG<<G4endl;
01659 #endif
01660 if(!pPDG || pPDG==9 || pPDG==21)
01661 {
01662 #ifdef debug
01663 G4cout<<"-Warning-G4QContent::AddParton: ImpossibleToAdd PartonWithPDG="<<pPDG<<G4endl;
01664 #endif
01665 return 0;
01666 }
01667 G4int aPDG = std::abs(pPDG);
01668 if( (aPDG>3 && aPDG<1101) || pPDG>3303)
01669 {
01670 #ifdef debug
01671 G4cout<<"-Warning-G4QContent::AddParton: Impossible Parton with PDG="<<pPDG<<G4endl;
01672 #endif
01673 return 0;
01674 }
01675 G4int HBN = GetBaryonNumber();
01676 if( HBN > 1 || HBN <-1)
01677 {
01678 #ifdef debug
01679 G4cout<<"-Warning-G4QContent::AddParton: Impossible Hadron with BaryonN="<<HBN<<G4endl;
01680 #endif
01681 return 0;
01682 }
01683 G4int AD=nAD;
01684 G4int AU=nAU;
01685 G4int AS=nAS;
01686 G4int QD=nD;
01687 G4int QU=nU;
01688 G4int QS=nS;
01689 if(aPDG>99)
01690 {
01691 G4int rPDG=aPDG/100;
01692 G4int P1=rPDG/10;
01693 G4int P2=rPDG%10;
01694 #ifdef debug
01695 G4cout<<"G4QContent::AddParton: DiQuark/AntiDiQuark, P1="<<P1<<", P2="<<P2<<G4endl;
01696 #endif
01697 if(pPDG>0)
01698 {
01699 #ifdef debug
01700 G4cout<<"G4QContent::AddParton: DiQuark, P1="<<P1<<", P2="<<P2<<",HBN="<<HBN<<G4endl;
01701 #endif
01702 if (P1==3 && P2==3)
01703 {
01704 if(HBN<0 && AS>1) AS-=2;
01705 else if(!HBN && AS==1)
01706 {
01707 AS=0;
01708 ++QS;
01709 }
01710 else if(HBN || (!HBN && !AS)) return 0;
01711 }
01712 else if(P1==3 && P2==2)
01713 {
01714 if(HBN<0 && AS && AU)
01715 {
01716 --AS;
01717 --AU;
01718 }
01719 else if(!HBN && (AS || AU))
01720 {
01721 if(AS)
01722 {
01723 --AS;
01724 ++QU;
01725 }
01726 else
01727 {
01728 --AU;
01729 ++QS;
01730 }
01731 }
01732 else if(HBN || (!HBN && !AS && !AU)) return 0;
01733 }
01734 else if(P1==3 && P2==1)
01735 {
01736 if(HBN<0 && AS && AD)
01737 {
01738 --AS;
01739 --AD;
01740 }
01741 else if(!HBN && (AS || AD))
01742 {
01743 if(AS)
01744 {
01745 --AS;
01746 ++QD;
01747 }
01748 else
01749 {
01750 --AD;
01751 ++QS;
01752 }
01753 }
01754 else if(HBN || (!HBN && !AS && !AD)) return 0;
01755 }
01756 else if(P1==2 && P2==2)
01757 {
01758 if(HBN<0 && AU>1) AU-=2;
01759 else if(!HBN && AU==1)
01760 {
01761 AU=0;
01762 ++QU;
01763 }
01764 else if(HBN || (!HBN && !AU)) return 0;
01765 }
01766 else if(P1==2 && P2==1)
01767 {
01768 if(HBN<0 && AD && AU)
01769 {
01770 --AD;
01771 --AU;
01772 }
01773 else if(!HBN && (AD || AU))
01774 {
01775 if(AD)
01776 {
01777 --AD;
01778 ++QU;
01779 }
01780 else
01781 {
01782 --AU;
01783 ++QD;
01784 }
01785 }
01786 else if(HBN || (!HBN && !AU && !AD)) return 0;
01787 }
01788 else
01789 {
01790 if(HBN<0 && AD>1) AD-=2;
01791 else if(!HBN && AD==1)
01792 {
01793 AD=0;
01794 ++QD;
01795 }
01796 else if(HBN || (!HBN && !AD)) return 0;
01797 }
01798 #ifdef debug
01799 G4cout<<"G4QContent::AddParton: DQ, QC="<<QD<<","<<QU<<","<<QS<<","<<AD<<","<<AU<<","
01800 <<AS<<G4endl;
01801 #endif
01802 if (HBN<0)
01803 {
01804 if (AD) return -1;
01805 else if(AU) return -2;
01806 else return -3;
01807 }
01808 else
01809 {
01810 if (QS)
01811 {
01812 if (QS==2) return 3303;
01813 else if(QU) return 3201;
01814 else return 3101;
01815 }
01816 else if(QU)
01817 {
01818 if (QU==2) return 2203;
01819 else return 2101;
01820 }
01821 else return 1103;
01822 }
01823 }
01824 else
01825 {
01826 #ifdef debug
01827 G4cout<<"G4QContent::AddParton: AntiDiQuark,P1="<<P1<<",P2="<<P2<<",B="<<HBN<<G4endl;
01828 #endif
01829 if (P1==3 && P2==3)
01830 {
01831 if(HBN>0 && QS>1) QS-=2;
01832 else if(!HBN && QS==1)
01833 {
01834 QS=0;
01835 ++AS;
01836 }
01837 else if(HBN || (!HBN && !QS)) return 0;
01838 }
01839 else if(P1==3 && P2==2)
01840 {
01841 if(HBN>0 && QS && QU)
01842 {
01843 --QS;
01844 --QU;
01845 }
01846 else if(!HBN && (QS || QU))
01847 {
01848 if(QS)
01849 {
01850 --QS;
01851 ++AU;
01852 }
01853 else
01854 {
01855 --QU;
01856 ++AS;
01857 }
01858 }
01859 else if(HBN || (!HBN && !QS && !QU)) return 0;
01860 }
01861 else if(P1==3 && P2==1)
01862 {
01863 if(HBN>0 && QS && QD)
01864 {
01865 --QS;
01866 --QD;
01867 }
01868 else if(!HBN && (QS || QD))
01869 {
01870 if(QS)
01871 {
01872 --QS;
01873 ++AD;
01874 }
01875 else
01876 {
01877 --QD;
01878 ++AS;
01879 }
01880 }
01881 else if(HBN || (!HBN && !QS && !QD)) return 0;
01882 }
01883 else if(P1==2 && P2==2)
01884 {
01885 if(HBN>0 && QU>1) QU-=2;
01886 else if(!HBN && QU==1)
01887 {
01888 QU=0;
01889 ++AU;
01890 }
01891 else if(HBN || (!HBN && !QU)) return 0;
01892 }
01893 else if(P1==2 && P2==1)
01894 {
01895 if(HBN>0 && QU && QD)
01896 {
01897 --QU;
01898 --QD;
01899 }
01900 else if(!HBN && (QU || QD))
01901 {
01902 if(QU)
01903 {
01904 --QU;
01905 ++AD;
01906 }
01907 else
01908 {
01909 --QD;
01910 ++AU;
01911 }
01912 }
01913 else if(HBN || (!HBN && !QU && !QD)) return 0;
01914 }
01915 else
01916 {
01917 if(HBN>0 && QD>1) QD-=2;
01918 else if(!HBN && QD==1)
01919 {
01920 QD=0;
01921 ++AD;
01922 }
01923 else if(HBN || (!HBN && !QD)) return 0;
01924 }
01925 #ifdef debug
01926 G4cout<<"G4QContent::AddParton:ADQ, QC="<<QD<<","<<QU<<","<<QS<<","<<AD<<","<<AU<<","
01927 <<AS<<G4endl;
01928 #endif
01929 if (HBN>0)
01930 {
01931 if (QD) return 1;
01932 else if(QU) return 2;
01933 else return 3;
01934 }
01935 else
01936 {
01937 if (AS)
01938 {
01939 if (AS==2) return -3303;
01940 else if(AU) return -3201;
01941 else return -3101;
01942 }
01943 else if(AU)
01944 {
01945 if (AU==2) return -2203;
01946 else return -2101;
01947 }
01948 else return -1103;
01949 }
01950 }
01951 }
01952 else
01953 {
01954 if(pPDG>0)
01955 {
01956 #ifdef debug
01957 G4cout<<"G4QContent::AddParton: Quark, A="<<AD<<","<<AU<<","<<AS<<",B="<<HBN<<G4endl;
01958 #endif
01959 if (aPDG==1)
01960 {
01961 if(HBN<0 && AD) AD--;
01962 else if(HBN || (!HBN && !AD)) return 0;
01963 }
01964 else if(aPDG==2)
01965 {
01966 if(HBN<0 && AU) AU--;
01967 else if(HBN || (!HBN && !AU)) return 0;
01968 }
01969 else
01970 {
01971 if(HBN<0 && AS) AS--;
01972 else if(HBN || (!HBN && !AS)) return 0;
01973 }
01974 #ifdef debug
01975 G4cout<<"G4QContent::AddParton: Q, QC="<<QD<<","<<QU<<","<<QS<<","<<AD<<","<<AU<<","
01976 <<AS<<G4endl;
01977 #endif
01978 if (!HBN)
01979 {
01980 if (QD) return 1;
01981 else if(QU) return 2;
01982 else return 3;
01983 }
01984 else
01985 {
01986 if (AS)
01987 {
01988 if (AS==2) return -3303;
01989 else if(AU) return -3201;
01990 else return -3101;
01991 }
01992 else if(AU)
01993 {
01994 if (AU==2) return -2203;
01995 else return -2101;
01996 }
01997 else return -1103;
01998 }
01999 }
02000 else
02001 {
02002 #ifdef debug
02003 G4cout<<"G4QContent::AddParton: antiQ, Q="<<QD<<","<<QU<<","<<QS<<",B="<<HBN<<G4endl;
02004 #endif
02005 if (aPDG==1)
02006 {
02007 if(HBN>0 && QD) QD--;
02008 else if(HBN || (!HBN && !QD)) return 0;
02009 }
02010 else if(aPDG==2)
02011 {
02012 if(HBN>0 && QU) QU--;
02013 else if(HBN || (!HBN && !QU)) return 0;
02014 }
02015 else
02016 {
02017 if(HBN>0 && QS) QS--;
02018 else if(HBN || (!HBN && !QS)) return 0;
02019 }
02020 #ifdef debug
02021 G4cout<<"G4QContent::AddParton: AQ, QC="<<QD<<","<<QU<<","<<QS<<","<<AD<<","<<AU<<","
02022 <<AS<<G4endl;
02023 #endif
02024 if (!HBN)
02025 {
02026 if (AD) return -1;
02027 else if(AU) return -2;
02028 else return -3;
02029 }
02030 else
02031 {
02032 if (QS)
02033 {
02034 if (QS==2) return 3303;
02035 else if(QU) return 3201;
02036 else return 3101;
02037 }
02038 else if(QU)
02039 {
02040 if (QU==2) return 2203;
02041 else return 2101;
02042 }
02043 else return 1103;
02044 }
02045 }
02046 }
02047 }