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
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066 #include <cmath>
00067 #include <cstdlib>
00068
00069 #include "G4Quasmon.hh"
00070 #include "G4SystemOfUnits.hh"
00071
00072 using namespace std;
00073
00074 G4Quasmon::G4Quasmon(G4QContent qQCont, G4LorentzVector q4M, G4LorentzVector ph4M)
00075 : q4Mom(q4M), valQ(qQCont), theWorld(0), phot4M(ph4M), f2all(0), rEP(0.), rMo(0.)
00076 {
00077 #ifdef debug
00078 G4cout<<"G4Quasmon:Constructor:QC="<<qQCont<<",Q4M="<<q4M<<",photonE="<<ph4M.e()<<G4endl;
00079 #endif
00080 #ifdef pardeb
00081 G4cout<<"**>G4Q:Con:(1),T="<<Temperature<<",S="<<SSin2Gluons<<",E="<<EtaEtaprime<<G4endl;
00082 #endif
00083 if(phot4M.e()>0.) q4Mom+=phot4M;
00084 valQ.DecQAQ(-1);
00085 status=4;
00086 }
00087
00088 G4Quasmon::G4Quasmon(const G4Quasmon& right): totMass(0), bEn(0), mbEn(0)
00089 {
00090 q4Mom = right.q4Mom;
00091 valQ = right.valQ;
00092
00093 status = right.status;
00094
00095 G4int nQH = right.theQHadrons.size();
00096 if(nQH) for(G4int ih=0; ih<nQH; ih++)
00097 {
00098 G4QHadron* curQH = new G4QHadron(right.theQHadrons[ih]);
00099 theQHadrons.push_back(curQH);
00100 }
00101 theWorld = right.theWorld;
00102 phot4M = right.phot4M;
00103 nBarClust = right.nBarClust;
00104 nOfQ = right.nOfQ;
00105
00106 G4int nQC = right.theQCandidates.size();
00107 if(nQC) for(G4int iq=0; iq<nQC; iq++)
00108 {
00109 G4QCandidate* curQC = new G4QCandidate(right.theQCandidates[iq]);
00110 theQCandidates.push_back(curQC);
00111 }
00112 f2all = right.f2all;
00113 rEP = right.rEP;
00114 rMo = right.rMo;
00115 }
00116
00117 G4Quasmon::G4Quasmon(G4Quasmon* right): totMass(0), bEn(0), mbEn(0)
00118 {
00119 #ifdef sdebug
00120 G4cout<<"G4Quasmon::Copy-Constructor: ***CALLED*** E="<<right->theEnvironment<<G4endl;
00121 #endif
00122 q4Mom = right->q4Mom;
00123 valQ = right->valQ;
00124
00125 status = right->status;
00126
00127 G4int nQH = right->theQHadrons.size();
00128 #ifdef sdebug
00129 G4cout<<"G4Quasmon::Copy-Constructor:nQH="<<nQH<<G4endl;
00130 #endif
00131 if(nQH) for(G4int ih=0; ih<nQH; ih++)
00132 {
00133 #ifdef debug
00134 G4cout<<"G4Quasmon:Copy-Constructor:H#"<<ih<<",QH="<<right->theQHadrons[ih]<<G4endl;
00135 #endif
00136 G4QHadron* curQH = new G4QHadron(right->theQHadrons[ih]);
00137 theQHadrons.push_back(curQH);
00138 }
00139 theWorld = right->theWorld;
00140 phot4M = right->phot4M;
00141 nBarClust = right->nBarClust;
00142 nOfQ = right->nOfQ;
00143
00144 G4int nQC = right->theQCandidates.size();
00145 #ifdef sdebug
00146 G4cout<<"G4Quasmon:Copy-Constructor: nCand="<<nQC<<G4endl;
00147 #endif
00148 if(nQC) for(G4int iq=0; iq<nQC; iq++)
00149 {
00150 #ifdef sdebug
00151 G4cout<<"G4Quasmon:Copy-Constructor:C#"<<iq<<",QC="<<right->theQCandidates[iq]<<G4endl;
00152 #endif
00153 G4QCandidate* curQC = new G4QCandidate(right->theQCandidates[iq]);
00154 theQCandidates.push_back(curQC);
00155 }
00156 f2all = right->f2all;
00157 rEP = right->rEP;
00158 rMo = right->rMo;
00159 #ifdef sdebug
00160 G4cout<<"G4Quasmon:Copy-Constructor: >->-> DONE <-<-<"<<G4endl;
00161 #endif
00162 }
00163
00164 G4Quasmon::~G4Quasmon()
00165 {
00166 #ifdef sdebug
00167 G4cout<<"G4Quasmon::Destructor before theQCandidates delete"<<G4endl;
00168 #endif
00169 for_each(theQCandidates.begin(), theQCandidates.end(), DeleteQCandidate());
00170 #ifdef sdebug
00171 G4cout<<"G4Quasmon::Destructor before theQHadrons"<<G4endl;
00172 #endif
00173 for_each(theQHadrons.begin(), theQHadrons.end(), DeleteQHadron());
00174 #ifdef sdebug
00175 G4cout<<"G4Quasmon::Destructor === DONE ==="<<G4endl;
00176 #endif
00177 }
00178
00179 G4double G4Quasmon::Temperature=180.;
00180 G4double G4Quasmon::SSin2Gluons=0.1;
00181 G4double G4Quasmon::EtaEtaprime=0.3;
00182 G4bool G4Quasmon::WeakDecays=false;
00183 G4bool G4Quasmon::ElMaDecays=true;
00184
00185
00186 void G4Quasmon::OpenElectromagneticDecays(){ElMaDecays=true;}
00187
00188
00189 void G4Quasmon::CloseElectromagneticDecays(){ElMaDecays=false;}
00190
00191
00192 void G4Quasmon::SetParameters(G4double temperature, G4double ssin2g, G4double etaetap)
00193 {
00194 Temperature=temperature;
00195 SSin2Gluons=ssin2g;
00196 EtaEtaprime=etaetap;
00197 }
00198 void G4Quasmon::SetTemper(G4double temperature) {Temperature=temperature;}
00199 void G4Quasmon::SetSOverU(G4double ssin2g) {SSin2Gluons=ssin2g;}
00200 void G4Quasmon::SetEtaSup(G4double etaetap) {EtaEtaprime=etaetap;}
00201
00202 const G4Quasmon& G4Quasmon::operator=(const G4Quasmon& right)
00203 {
00204 if(this != &right)
00205 {
00206 q4Mom = right.q4Mom;
00207 valQ = right.valQ;
00208
00209 status = right.status;
00210
00211 G4int iQH = theQHadrons.size();
00212 if(iQH) for(G4int jh=0; jh<iQH; jh++) delete theQHadrons[jh];
00213 theQHadrons.clear();
00214 G4int nQH = right.theQHadrons.size();
00215 if(nQH) for(G4int ih=0; ih<nQH; ih++)
00216 {
00217 G4QHadron* curQH = new G4QHadron(right.theQHadrons[ih]);
00218 theQHadrons.push_back(curQH);
00219 }
00220 theWorld = right.theWorld;
00221 phot4M = right.phot4M;
00222 nBarClust = right.nBarClust;
00223 nOfQ = right.nOfQ;
00224
00225 G4int iQC = theQCandidates.size();
00226 if(iQC) for(G4int jq=0; jq<iQC; jq++) delete theQCandidates[jq];
00227 theQCandidates.clear();
00228 G4int nQC = right.theQCandidates.size();
00229 if(nQC) for(G4int iq=0; iq<nQC; iq++)
00230 {
00231 G4QCandidate* curQC = new G4QCandidate(right.theQCandidates[iq]);
00232 theQCandidates.push_back(curQC);
00233 }
00234 f2all = right.f2all;
00235 rEP = right.rEP;
00236 rMo = right.rMo;
00237 }
00238 return *this;
00239 }
00240
00241
00242 G4QHadronVector G4Quasmon::HadronizeQuasmon(G4QNucleus& qEnv, G4int nQuasms)
00243 {
00245 static const G4int NUCPDG = 90000000;
00246 static const G4int MINPDG = 80000000;
00247 static const G4double np = 1877.9;
00248 static const G4double BIG = 1000000.;
00249 static const G4double BIG2 = BIG*BIG;
00250 static const G4QContent neutQC(2,1,0,0,0,0);
00251 static const G4QContent protQC(1,2,0,0,0,0);
00252 static const G4QContent lambQC(1,1,1,0,0,0);
00253 static const G4QContent PiQC(0,1,0,1,0,0);
00254 static const G4QContent PiMQC(1,0,0,0,1,0);
00255 static const G4QContent Pi0QC(0,1,0,0,1,0);
00256 static const G4QContent K0QC(1,0,0,0,0,1);
00257 static const G4QContent KpQC(0,1,0,0,0,1);
00258 static const G4QContent zeroQC(0,0,0,0,0,0);
00259
00260 static const G4LorentzVector zeroLV(0.,0.,0.,0.);
00261 static const G4QNucleus vacuum= G4QNucleus(zeroLV,NUCPDG);
00262 static const G4double mDeut= G4QPDGCode(2112).GetNuclMass(1,1,0);
00263 static const G4double mTrit= G4QPDGCode(2112).GetNuclMass(1,2,0);
00264 static const G4double mHel3= G4QPDGCode(2112).GetNuclMass(2,1,0);
00265 static const G4double mAlph= G4QPDGCode(2112).GetNuclMass(2,2,0);
00267 static const G4double mNeut= G4QPDGCode(2112).GetMass();
00269 static const G4double mProt= G4QPDGCode(2212).GetMass();
00270 static const G4double mLamb= G4QPDGCode(3122).GetMass();
00271 static const G4double mPi = G4QPDGCode(211).GetMass();
00272 static const G4double mPi0 = G4QPDGCode(111).GetMass();
00273 static const G4double mK = G4QPDGCode(321).GetMass();
00274 static const G4double mK0 = G4QPDGCode(311).GetMass();
00275 static const G4double mEta = G4QPDGCode(221).GetMass();
00276 static const G4double mEtaP= G4QPDGCode(331).GetMass();
00277 static const G4double diPiM= mPi0 + mPi0;
00278 static const G4double PiNM = mPi + mNeut;
00280 static const G4double mPi02= mPi0* mPi0;
00286 static const G4double mP2 = mProt*mProt;
00292 static const G4double eps = 0.003;
00293 G4double momPhoton=phot4M.rho();
00294 G4double addPhoton=phot4M.e();
00295 #ifdef debug
00296 G4cout<<"G4Quasmon::HadrQ:*==>>START QUASMON HADRONIZATION<<==*, aP="<<addPhoton<<",Env="
00297 <<qEnv<<qEnv.GetProbability()<<", #ofQuasms="<<nQuasms<<G4endl;
00298 #endif
00299 G4bool first=false;
00300 if(nQuasms<0)
00301 {
00302 if(nQuasms==-1) first=true;
00303 else G4cout<<"G4Quasmon::HadrQ: Negative #of Quasmons n="<<nQuasms<<G4endl;
00304 nQuasms=-nQuasms;
00305 }
00306 G4bool piF = false;
00307 G4bool gaF = false;
00308 if(addPhoton<0.)
00309 {
00310 #ifdef debug
00311 G4cout<<"G4Q::HQ: PionAtRest, addP="<<addPhoton<<", momP="<<momPhoton<<G4endl;
00312 #endif
00313 addPhoton=0.;
00314 momPhoton=0.;
00315 piF=true;
00316 }
00317 else if(addPhoton>0.) gaF=true;
00318 theEnvironment=qEnv;
00319 #ifdef debug
00320 G4cout<<"G4Quasmon::HadrQ:Env="<<theEnvironment<<theEnvironment.GetProbability()<<G4endl;
00321 #endif
00322
00323 theWorld= G4QCHIPSWorld::Get();
00324 G4int nP= theWorld->GetQPEntries();
00325
00326 G4int nMesons = G4QNucleus().GetNDefMesonC();
00327 #ifdef debug
00328 G4cout<<"G4Quasmon::HadrQ:CHIPSWorld initialized with nP="<<nP<<",nM="<<nMesons<<G4endl;
00329 #endif
00330 if (nP<34) nMesons = 9;
00331 else if(nP<51) nMesons = 18;
00332 else if(nP<65) nMesons = 27;
00333 else if(nP<82) nMesons = 36;
00334 G4int nBaryons = G4QNucleus().GetNDefBaryonC();
00335 if (nP<45) nBaryons = 16;
00336 else if(nP<59) nBaryons = 36;
00337 else if(nP<76) nBaryons = 52;
00338 G4int nClusters= nP-G4QPDGCode().GetNQHadr();
00339 #ifdef debug
00340 G4cout<<"G4Quasmon:HadrQ: Init Candidates:"<<theEnvironment<<",n="<<theQCandidates.size()
00341 <<",nMesons="<<nMesons<<",nBaryons="<<nBaryons<<",nClusters="<<nClusters<<G4endl;
00342 #endif
00343 theEnvironment.InitCandidateVector(theQCandidates,nMesons,nBaryons,nClusters);
00344 #ifdef debug
00345 G4cout<<"G4Quasmon:HadrQ:CandidatesAreInitialized,n="<<theQCandidates.size()<<",nMesons="
00346 <<nMesons<<", nBaryons="<<nBaryons<<", nClusters="<<nClusters<<G4endl;
00347 #endif
00348 if(!status||q4Mom==zeroLV)
00349 {
00350 #ifdef debug
00351 G4cout<<"G4Q::HQ:NOTHING-TO-DO: Q4M="<<q4Mom<<", QEnv="<<theEnvironment<<G4endl;
00352 #endif
00353 if(addPhoton)
00354 {
00355 G4ExceptionDescription ed;
00356 ed << "OverheadPhoton for theZeroQuasmon: Q4M=" << q4Mom << ",status="
00357 << status << ", phE=" << addPhoton << G4endl;
00358 G4Exception("G4Quasmon::HadronizeQuasmon()", "HAD_CHPS_0000", FatalException, ed);
00359 }
00360 KillQuasmon();
00361 qEnv=theEnvironment;
00362 return theQHadrons;
00363 }
00364 status=2;
00365 G4int sPDG=0;
00366 G4int pPDG=NUCPDG;
00367
00369 G4double npqp2=0;
00370 G4double sMass=0.;
00371 G4double sM2=0.;
00372 G4int pBaryn=0;
00373 G4double dMass=0.;
00374 G4double pMass=0.;
00375 G4double pNMass=0.;
00376 G4double delta=0.;
00377 G4double deltaN=0.;
00378 G4double minT=0.;
00379 G4double minSqT=0.;
00380 G4double minSqB=0.;
00381 G4double minSqN=0.;
00382 G4double hili=0.;
00383 G4double loli=0.;
00384 G4double tmpTM2=BIG2;
00385 G4double reTNM2=0.;
00386 G4QContent curQ=zeroQC;
00387 G4QContent memQ=zeroQC;
00388 G4QContent pQC=zeroQC;
00389 G4QContent sQC=zeroQC;
00390 G4QContent transQC=zeroQC;
00391 G4LorentzVector m4Mom=zeroLV;
00392 G4LorentzVector kp4Mom=zeroLV;
00393 G4LorentzVector check=-theEnvironment.Get4Momentum()-q4Mom;
00394 G4int ccheck=-theEnvironment.GetZ()-valQ.GetCharge();
00395 #ifdef chdebug
00396 G4int cSum=-ccheck;
00397 G4QNucleus oldEnv(theEnvironment);
00398 G4QContent oldCQC(valQ);
00399 G4int oldNH=theQHadrons.size();
00400 #endif
00401 G4bool start=true;
00402 #ifdef debug
00403 G4cout<<"G4Q::HQ: Before the loop EnvPDG="<<theEnvironment.GetPDG()<<G4endl;
00404 #endif
00405 while(theEnvironment.GetPDG()==NUCPDG || start)
00406 {
00407 #ifdef chdebug
00408 G4int ccSum=theEnvironment.GetZ()+valQ.GetCharge();
00409 G4int nHd=theQHadrons.size();
00410 if(nHd) for(int ih=0; ih<nHd; ih++) ccSum+=theQHadrons[ih]->GetCharge();
00411 if(ccSum!=cSum)
00412 {
00413 G4cerr<<"*G4Q::HQ:C"<<cSum<<",c="<<ccSum<<",E="<<theEnvironment<<",Q="<<valQ<<G4endl;
00414 G4cerr<<":G4Q::HQ:oldE="<<oldEnv<<"oldQ="<<oldCQC<<",oN="<<oldNH<<",N="<<nHd<<G4endl;
00415 if(nHd) for(int h=0; h<nHd; h++)
00416 {
00417 G4QHadron* cH = theQHadrons[h];
00418 G4cerr<<"::G4Q::HQ:#h"<<h<<",C="<<cH->GetCharge()<<",P="<<cH->GetPDGCode()<<G4endl;
00419 }
00420 }
00421 oldEnv=G4QNucleus(theEnvironment);
00422 oldCQC=G4QContent(valQ);
00423 oldNH=nHd;
00424 #endif
00425 start=false;
00426
00427 G4double qM2 = q4Mom.m2();
00428 G4double tmpEq=q4Mom.e();
00429 G4double tmpPq=q4Mom.rho();
00430 if(fabs(qM2)<.0001 || tmpEq<=tmpPq)
00431 {
00432 qM2=0.;
00433 if(!valQ.GetCharge() && !valQ.GetBaryonNumber() && !valQ.GetStrangeness())
00434 {
00435 if(fabs(qM2)<.001)
00436 {
00437 q4Mom.setE(tmpPq);
00438 #ifdef debug
00439 G4cout<<"G4Q::HQ:Quasmon is gamma, Q4M="<<q4Mom<<",E="<<theEnvironment<<G4endl;
00440 #endif
00441 G4QHadron* gamH = new G4QHadron(22,q4Mom);
00442 FillHadronVector(gamH);
00443 KillQuasmon();
00444 qEnv=theEnvironment;
00445 return theQHadrons;
00446 }
00447 else if(tmpPq<.001)
00448 {
00449 #ifdef debug
00450 G4cout<<"G4Q::HQ:Quasmon is nothing, Q4M="<<q4Mom<<",E="<<theEnvironment<<G4endl;
00451 #endif
00452 KillQuasmon();
00453 qEnv=theEnvironment;
00454 return theQHadrons;
00455 }
00456 }
00457 else q4Mom.setE(tmpPq*1.00001);
00458 }
00459 G4double quasM= sqrt(qM2);
00460 G4double qurF = quasM/(tmpEq-tmpPq);
00461 G4ThreeVector qltb=q4Mom.boostVector();
00463 #ifdef debug
00464 G4cout<<"G4Q::HQ: Quasm="<<q4Mom<<",qM="<<quasM<<",qQC="<<valQ<<G4endl;
00465 #endif
00466 CalculateNumberOfQPartons(quasM);
00467
00468 G4int envPDG=theEnvironment.GetPDG();
00469 G4int envN =theEnvironment.GetN();
00470 G4int envZ =theEnvironment.GetZ();
00471 G4int envS =theEnvironment.GetS();
00472 G4int envA =theEnvironment.GetA();
00473
00474 G4int maxActEnv=256;
00475 G4int dmaxActEnv=maxActEnv+maxActEnv;
00476
00477 if(envA>dmaxActEnv)
00478
00479 {
00480
00481 G4int zEn=maxActEnv;
00482 G4int nEn=zEn;
00483 bEn = zEn+nEn;
00484 mbEn = G4QPDGCode(2112).GetNuclMass(zEn,nEn,0);
00485 bEn4M=G4LorentzVector(0.,0.,0.,mbEn);
00486 bEnQC=G4QContent(bEn+nEn,bEn+zEn,0,0,0,0);
00487 }
00488 else
00489 {
00490 bEn = 256;
00491 mbEn = G4QPDGCode(2112).GetNuclMass(128,128,0);
00492 bEn4M = G4LorentzVector(0.,0.,0.,mbEn);
00493 bEnQC = G4QContent(384,384,0,0,0,0);
00494 }
00495 G4double envM=theEnvironment.GetMass();
00496 G4QContent envQC=theEnvironment.GetQCZNS();
00497 #ifdef debug
00498 G4cout<<"G4Q::HQ: ePDG="<<envPDG<<",eM="<<envM<<",eQC="<<envQC<<G4endl;
00499 #endif
00500 G4QContent totQC=valQ+envQC;
00501 G4int totBN=totQC.GetBaryonNumber();
00502 G4int totS=totQC.GetStrangeness();
00503 G4int totZ=totQC.GetCharge();
00504 G4QNucleus totN(totQC);
00505 G4int totNeut=totN.GetN();
00507 G4double totM =totN.GetMZNS();
00508 #ifdef debug
00509 G4cout<<"G4Q::HQ: tN="<<totN<<",tGSM="<<totM<<",tQC="<<totQC<<G4endl;
00510 #endif
00512 #ifdef debug
00513 G4int resNPDG=0;
00514 G4double resNM =10000000.;
00515 #endif
00516 if(totNeut>0)
00517 {
00518 G4QContent resNQC=totQC-G4QContent(2,1,0,0,0,0);
00519 G4QNucleus resNN(resNQC);
00520 #ifdef debug
00521 resNM = resNN.GetMZNS();
00522 resNPDG= resNN.GetPDG();
00523 #endif
00524 }
00525 G4LorentzVector env4M =G4LorentzVector(0.,0.,0.,envM);
00526 G4LorentzVector tot4M =q4Mom+env4M;
00527 totMass=tot4M.m();
00528 G4int totPDG=totN.GetPDG();
00529 #ifdef debug
00530 G4cout<<"G4Q::HQ: totPDG="<<totPDG<<",totM="<<totMass<<",rPDG="<<resNPDG<<G4endl;
00531 #endif
00532 G4double totEn=tot4M.e();
00533 G4double totMo=tot4M.rho();
00534 if(totEn<totMo)
00535 {
00536 G4cerr<<"---Warning---G4Q::HQ: *Boost* tot4M="<<tot4M<<", E-p="<<totEn-totMo<<G4endl;
00537 G4double accuracy=.000001*totMo;
00538 G4double emodif=fabs(totEn-totMo);
00539
00540
00541 G4cerr<<"G4Q::HQ: *Boost* E-p shift is corrected to "<<emodif<<G4endl;
00542 tot4M.setE(totMo+emodif+.01*accuracy);
00543
00544 }
00545 G4ThreeVector totBoost = tot4M.boostVector();
00546 G4ThreeVector totRBoost= -totBoost;
00547 G4int iniPDG =valQ.GetSPDGCode();
00548 G4int iniBN =valQ.GetBaryonNumber();
00549 G4int iniQChg=valQ.GetCharge();
00550 G4int iniN =valQ.GetN();
00551 G4int iniP =valQ.GetP();
00552 G4int iniS =valQ.GetL();
00553 #ifdef debug
00554 G4cout<<"G4Q::HQ: iniPDG="<<iniPDG<<", Z="<<iniP<<", N="<<iniN<<", S="<<iniS<<G4endl;
00555 #endif
00556 G4QNucleus iniRN(iniP,iniN-1,iniS);
00557 G4double iniQM =G4QPDGCode(iniPDG).GetMass();
00558 #ifdef debug
00559 G4double iniRM = iniRN.GetMZNS();
00560 if(iniBN<2||envA>0) iniRM=0.;
00561 G4cout<<"G4Q::HQ: iniRN="<<iniRN<<", iniRM="<<iniRM<<", iniQM="<<iniQM<<G4endl;
00562 #endif
00564 G4double excE = totMass-totM;
00566 #ifdef debug
00567 G4double bndQM = totM-envM;
00568 if(envPDG==NUCPDG) bndQM=iniQM;
00569 G4double bndQM2= bndQM*bndQM;
00570 G4double quen = iniQM+envM;
00571 G4cout<<"G4Q::HQ:mQ="<<quasM<<valQ<<bndQM2<<",nQ="<<nOfQ<<",Env="<<envM
00572 <<envQC<<",Q+E="<<quen<<",tM="<<totPDG<<totQC<<totM<<"<"<<totMass<<G4endl;
00573 #endif
00574 G4int tQ = valQ.GetTot();
00575 G4int bQ = abs(valQ.GetBaryonNumber());
00576 G4QContent cQ = valQ;
00577 G4int s_value = 4;
00578 if (bQ) s_value = 3*bQ + 2;
00579 if (tQ> s_value) cQ.DecQAQ((tQ-s_value)/2);
00580 #ifdef debug
00581 G4int rsPDG = cQ.GetSPDGCode();
00582 G4cout<<"G4Q::HQ:eN="<<envN<<",eZ="<<envZ<<",Q="<<rsPDG<<cQ<<",piF="<<piF<<",gaF="<<gaF
00583 <<G4endl;
00584 #endif
00585 theEnvironment.UpdateClusters(false);
00586
00587
00588
00589 theEnvironment.PrepareCandidates(theQCandidates,piF,piF);
00590 ModifyInMatterCandidates();
00591 G4double kMom = 0.;
00592 G4double minK = 0.;
00593 G4double maxK = quasM/2.;
00594 G4double kLS = 0;
00595 G4double cost = 0.;
00596 G4bool kCond = true;
00597 G4bool qCond = true;
00598 G4bool pCond = true;
00599 G4bool fskip=false;
00600
00601
00602 G4LorentzVector k4Mom=zeroLV;
00603 G4LorentzVector cr4Mom=zeroLV;
00604 G4int kCount =0;
00605
00606
00607
00608
00609
00610 G4int qCountMax=1;
00611 if(excE > diPiM) qCountMax=(G4int)(excE/mPi0);
00612
00613
00614
00615
00616 G4int kCountMax=1;
00617
00618
00619
00620
00621
00622
00623 G4int pCountMax=1;
00624
00625 if(envA>0&&envA<19) pCountMax=36/envA;
00626
00627
00628 G4bool gintFlag=false;
00629 while(kCount<kCountMax&&kCond)
00630 {
00631 kCond=true;
00632 G4double miM2=0.;
00633 if(envPDG==NUCPDG)
00634 {
00635 if(excE>mPi0) miM2=mPi02;
00636 else miM2=mP2;
00637 }
00638 else
00639 {
00640 minK=100000.;
00641
00642
00643
00644
00645
00646
00647
00648
00649
00650
00651
00652
00653
00654
00655
00656
00657
00658
00659
00660
00661
00662 #ifdef debug
00663
00664
00665
00666
00667 #endif
00668
00669
00670
00671
00672
00673
00674
00675
00676
00677
00678
00679
00680
00681 envA=envZ+envN;
00682 G4int totN_value=totBN-totZ;
00683 if(totN_value>0&&totBN>1)
00684 {
00685 G4QNucleus tmpNN(totZ,totN_value-1,0);
00686 G4double delN=tmpNN.GetMZNS()+mNeut-totM;
00687 if(envN>0&&envA>1)
00688 {
00689 G4QNucleus envNN(envZ,envN-1,0);
00690 G4double delEN=envNN.GetMZNS()+mNeut-envM;
00691 if(delEN>delN) delN=delEN;
00692 }
00693 delN*=qurF;
00694 if(delN<minK) minK=delN;
00695 }
00696 if(totZ>0&&totBN>1)
00697 {
00698 G4double proCB=theEnvironment.CoulombBarrier(1,1);
00699 G4QNucleus tmpPN(totZ-1,totN_value,0);
00700 G4double delP=tmpPN.GetMZNS()+mProt-totM+proCB;
00701 if(envZ>0&&envA>1)
00702 {
00703 G4QNucleus envPN(envZ-1,envN,0);
00704 G4double delEP=envPN.GetMZNS()+mProt-envM+proCB;
00705 if(delEP>delP) delP=delEP;
00706 }
00707 delP*=qurF;
00708 if(delP<minK) minK=delP;
00709 }
00710 if(totN_value>0&&totZ>0&&totBN>2)
00711 {
00712 G4double proCB=theEnvironment.CoulombBarrier(1,2);
00713 G4QNucleus tmpDN(totZ-1,totN_value-1,0);
00714 G4double delD=tmpDN.GetMZNS()+mDeut-totM+proCB;
00715 if(envN>0&&envZ>0&&envA>2)
00716 {
00717 G4QNucleus envDN(envZ-1,envN-1,0);
00718 G4double delED=envDN.GetMZNS()+mDeut-envM+proCB;
00719 if(delED>delD) delD=delED;
00720 }
00721 delD*=qurF;
00722 if(delD<minK) minK=delD;
00723 }
00724 if(totN_value>1&&totZ>0&&totBN>3)
00725 {
00726 G4double proCB=theEnvironment.CoulombBarrier(1,3);
00727 G4QNucleus tmpTN(totZ-1,totN_value-2,0);
00728 G4double delT=tmpTN.GetMZNS()+mTrit-totM+proCB;
00729 if(envN>1&&envZ>0&&envA>3)
00730 {
00731 G4QNucleus envTN(envZ-1,envN-2,0);
00732 G4double delET=envTN.GetMZNS()+mTrit-envM+proCB;
00733 if(delET>delT) delT=delET;
00734 }
00735 delT*=qurF;
00736 if(delT<minK) minK=delT;
00737 }
00738 if(totN_value>0&&totZ>1&&totBN>3)
00739 {
00740 G4double proCB=theEnvironment.CoulombBarrier(2,3);
00741 G4QNucleus tmpRN(totZ-2,totN_value-1,0);
00742 G4double delR=tmpRN.GetMZNS()+mHel3-totM+proCB;
00743 if(envN>0&&envZ>1&&envA>3)
00744 {
00745 G4QNucleus envRN(envZ-2,envN-1,0);
00746 G4double delER=envRN.GetMZNS()+mHel3-envM+proCB;
00747 if(delER>delR) delR=delER;
00748 }
00749 delR*=qurF;
00750 if(delR<minK) minK=delR;
00751 }
00752 if(totN_value>1&&totZ>1&&totBN>4)
00753 {
00754 G4double proCB=theEnvironment.CoulombBarrier(2,4);
00755 G4QNucleus tmpAN(totZ-2,totN_value-2,0);
00756 G4double delA=tmpAN.GetMZNS()+mAlph-totM+proCB;
00757 if(envN>1&&envZ>1&&envA>4)
00758 {
00759 G4QNucleus envAN(envZ-2,envN-2,0);
00760 G4double delEA=envAN.GetMZNS()+mAlph-envM+proCB;
00761 if(delEA>delA) delA=delEA;
00762 }
00763 delA*=qurF;
00764 if(delA*qurF<minK) minK=delA*qurF;
00765 }
00766
00767
00768
00769
00770
00771
00772
00773
00774
00775
00776
00777
00778
00779
00780
00781 if(minK<0. || minK+minK > quasM) minK=0.;
00782 }
00783
00784
00785
00786
00787
00788
00789
00790
00791
00793
00794 if(addPhoton>0.)
00795
00796 {
00797
00798
00799 gintFlag=true;
00800 q4Mom-= phot4M;
00801 qM2 = q4Mom.m2();
00802 quasM = sqrt(qM2);
00803 G4double kpow=static_cast<double>(nOfQ-2);
00804 if(kpow<1.) kpow=1.;
00805 G4double xmi=(momPhoton-addPhoton)*quasM;
00806 if(xmi<0.) xmi=0.;
00807
00808
00809
00810
00811
00812
00813
00814 G4double rn=G4UniformRand();
00815 kMom=(1.-(1.-xmi)*pow(rn,1./kpow))*quasM/2;
00816
00818 G4double dkM=kMom+kMom;
00819
00820
00821 G4double cor=200./dkM/addPhoton;
00822 G4double bas=std::log(2+cor);
00823 cost=1.-std::exp(bas-(bas-std::log(cor))*G4UniformRand());
00824
00825
00826
00827
00828
00829
00830
00831
00832 if(cost>1.) cost=1.;
00833 if(cost<-1.) cost=-1.;
00834 #ifdef debug
00835 G4cout<<"G4Q::HQ:**PHOTON out of Q**k="<<kMom<<",ct="<<cost<<",QM="<<quasM<<G4endl;
00836 #endif
00837
00838
00839
00840
00841
00842
00843
00844
00845
00846
00847
00848
00849
00850
00851
00852
00853
00854
00855
00856
00857
00858
00861
00862
00863
00864
00865
00866 }
00867 else
00868 {
00869 gaF=false;
00870 gintFlag=false;
00871
00872 if(!miM2) miM2=(minK+minK)*quasM;
00873 if(qM2<.0001) kMom=0.;
00874 else kMom = GetQPartonMomentum(maxK,miM2);
00875
00876
00877
00878
00879
00880
00881
00882
00883
00884
00885
00886 G4double rnc=G4UniformRand();
00887 cost = 1.-rnc-rnc;
00888 }
00889 G4double cQM2=qM2-(kMom+kMom)*quasM;
00890 if(cQM2<0.)
00891 {
00892
00893 cQM2=0.;
00894 }
00895 G4double cQM=sqrt(cQM2);
00896 k4Mom=zeroLV;
00897 cr4Mom=G4LorentzVector(0.,0.,0.,cQM);
00898 G4LorentzVector dir4M=q4Mom-G4LorentzVector(0.,0.,0.,q4Mom.e()*.01);
00899 if(!G4QHadron(q4Mom).RelDecayIn2(k4Mom,cr4Mom,dir4M,cost,cost))
00900 {
00901 #ifdef debug
00902 G4cerr<<"*G4Q::HQ:PB,M="<<quasM<<",cM="<<cQM<<",c="<<cost<<",F="<<gintFlag<<G4endl;
00903 #endif
00904 kCond=true;
00905 if(addPhoton&&gintFlag)
00906 {
00907 q4Mom+= phot4M;
00908 qM2 = q4Mom.m2();
00909 if(qM2>0.) quasM = sqrt(qM2);
00910 else
00911 {
00912 if(qM2<-.0001)G4cerr<<"--Warning-- G4Q::HQ:Phot.M2="<<qM2<<" Cor to 0"<<G4endl;
00913 quasM=0.;
00914 }
00915 gintFlag=false;
00916 }
00917 }
00918 else
00919 {
00920 if(addPhoton&&gintFlag)
00921 {
00922 q4Mom+= phot4M;
00923 k4Mom+= phot4M;
00924 qM2 = q4Mom.m2();
00925 quasM = sqrt(qM2);
00926 kMom=k4Mom.e();
00927 gintFlag=false;
00928 }
00929 #ifdef debug
00930 G4cout<<"G4Q::HQ:(PhBack) k="<<kMom<<",k4M="<<k4Mom<<",ct="<<cost<<",gF="<<gintFlag
00931 <<G4endl;
00932 #endif
00933 kLS=k4Mom.e();
00934 G4double rEn=cr4Mom.e();
00935 rMo=cr4Mom.rho();
00936 rEP=rEn+rMo;
00937 G4int totCand = theQCandidates.size();
00938 #ifdef sdebug
00939 G4cout<<"G4Q::HQ: ****>>K-ITERATION #"<<kCount<<", k="<<kMom<<k4Mom<<G4endl;
00940 #endif
00941 for (G4int index=0; index<totCand; index++)
00942 {
00943 G4QCandidate* curCand=theQCandidates[index];
00944 G4int cPDG = curCand->GetPDGCode();
00945 if(cPDG==90000001||cPDG==90001000||cPDG==91000000||cPDG<MINPDG)
00946 {
00947 G4bool poss= curCand->GetPossibility();
00948 #ifdef debug
00949 if(cPDG==90000001 || cPDG==90001000 || cPDG==90000002 || cPDG==90001001)
00950 G4cout<<"G4Q::HQ:pos="<<poss<<",cPDG="<<cPDG<<",iQC="<<iniQChg<<G4endl;
00951 #endif
00952 if(poss)
00953 {
00954 G4double cMs=curCand->GetEBMass();
00955 G4QContent cQC=curCand->GetQC();
00956 G4double cfM=curCand->GetQPDG().GetMass();
00957 G4QContent rtQC=curQ+envQC-cQC;
00958 G4QNucleus rtN(rtQC);
00962 #ifdef debug
00963 if(cPDG==90000001 || cPDG==90001000 || cPDG==90000002 || cPDG==90001001)
00964 G4cout<<"G4Q::HQ:cfM="<<cfM<<",cMs="<<cMs<<",ind="<<index<<G4endl;
00965 #endif
00966
00967
00968 G4double kMin=0.;
00969 if(cMs) kMin=(cfM*cfM-cMs*cMs)/(cMs+cMs);
00970 if(kMin<0.) kMin=0.;
00971 #ifdef debug
00972 G4double totr = rtN.GetMZNS();
00973 G4double bnM = totr-envM+cMs;
00974
00975
00976
00977 G4double dR=0.;
00978 if(cPDG==90000001 || cPDG==90001000 || cPDG==90000002 || cPDG==90001001)
00979 G4cout<<"G4Q::HQ:i="<<index<<",cPDG="<<cPDG<<",k="<<kMom<<","<<kLS<<">kMin="
00980 <<kMin<<",bM="<<bnM<<",rEP="<<rEP<<",dR="<<dR<<",kCo="<<kCond<<G4endl;
00981 #endif
00982 if(kLS>kMin)
00983 {
00984 kCond=false;
00985 break;
00986 }
00987 }
00988 }
00989 }
00990 }
00991 kCount++;
00992 }
00993
00994 #ifdef debug
00995 if(addPhoton)
00996 G4cout<<"G4Q::HQ:***PHOTON OK***k="<<k4Mom<<",Q="<<q4Mom<<",PhE="<<addPhoton<<G4endl;
00997 #endif
00998
00999 #ifdef debug
01000 G4cout<<"G4Q::HQ:Select="<<kMom<<",ki="<<minK<<",ka="<<maxK<<",k="<<k4Mom<<kLS<<G4endl;
01001 #endif
01002 CalculateHadronizationProbabilities(excE,kMom,k4Mom,piF,gaF,first);
01003
01004
01005 addPhoton=0.;
01006 momPhoton=0.;
01007
01008 G4double dk = kMom+kMom;
01009 G4int nCandid = theQCandidates.size();
01010 G4bool fprob = true;
01012 G4bool fdul = false;
01013 int i=0;
01014 G4double maxP = 0.;
01015 if(nCandid) maxP = theQCandidates[nCandid-1]->GetIntegProbability();
01016 #ifdef debug
01017 G4cout<<"G4Q::HQ:***RANDOMIZE CANDIDATEs***a#OfCand="<<nCandid<<",maxP="<<maxP<<G4endl;
01018 #endif
01019 if (maxP<=0.)
01020 {
01021 #ifdef debug
01022 if(status == 4) G4Exception("G4Quasmon::HadronizeQuasmon()", "HAD_CHPS_0001",
01023 FatalException, "HQ:*TMP EXCEPTION - NoChanalsOfFragmentation");
01024 #endif
01025 #ifdef debug
01026 G4cout<<"G4Q::HQ:Z="<<iniP<<",B="<<iniBN<<G4endl;
01027 #endif
01028 G4double qCB=theEnvironment.CoulombBarrier(iniP,iniBN);
01029 #ifdef debug
01030 G4cout<<"G4Q::HQ:qCB="<<qCB<<",Z="<<iniP<<",B="<<iniBN<<",eE="<<excE<<G4endl;
01031 #endif
01032 G4double pCB=theEnvironment.CoulBarPenProb(qCB,excE,iniP,iniBN);
01033 G4bool qenv=envPDG!=90000000&&envPDG!=90000002&&envPDG!=90002000&&iniPDG!=90002000
01034 &&iniPDG!=90000002;
01035 #ifdef debug
01036 G4cout<<"G4Q::HQ: qCB="<<qCB<<",pCB="<<pCB<<",cond="<<qenv<<",N="<<nQuasms<<G4endl;
01037 #endif
01038
01039 if(!first&&excE>qCB&&envM+iniQM<totMass&&nQuasms==1&&qenv&&G4UniformRand()<pCB)
01040
01041 {
01042 G4QHadron* resNuc = new G4QHadron(valQ,q4Mom);
01043 resNuc->CorMDecayIn2(iniQM,env4M);
01044
01045 FillHadronVector(resNuc);
01046 #ifdef debug
01047 G4LorentzVector oQ4Mom=resNuc->Get4Momentum();
01048 G4cout<<"G4Q::HQ: outQM="<<oQ4Mom.m()<<oQ4Mom<<",GSQM="<<iniQM<<G4endl;
01049 #endif
01050 if(nQuasms==1) qEnv = G4QNucleus(env4M,envPDG);
01051 else
01052 {
01053 G4QHadron* envH = new G4QHadron(envPDG,env4M);
01054 FillHadronVector(envH);
01055 qEnv = vacuum;
01056 }
01057 #ifdef debug
01058 G4double eM=env4M.m();
01059 G4LorentzVector dif=tot4M-oQ4Mom-env4M;
01060 G4cout<<"G4Q::HQ: envM="<<envM<<"=="<<eM<<", envT="<<env4M.e()-eM<<dif<<G4endl;
01061 #endif
01062 ClearQuasmon();
01063 return theQHadrons;
01064 }
01065 else
01066 {
01067 #ifdef debug
01068 G4cout<<"G4Q::HQ:Q2E:E="<<theEnvironment<<",valQ="<<valQ<<",tot4M="<<tot4M<<G4endl;
01069 #endif
01070
01071 if(CheckGroundState(true)) ClearQuasmon();
01072
01073
01074
01075 #ifdef debug
01076
01077 #endif
01078
01079
01080
01081 else if(envPDG==NUCPDG && quasM>iniQM)
01082 {
01083 #ifdef debug
01084 G4cout<<"***G4Q::HQ: Emergency Decay in Gamma/Pi0 + Residual GSQuasmon"<<G4endl;
01085
01086 #endif
01087 G4int gamPDG=22;
01088 G4double gamM=0.;
01089 if(quasM>mPi0+iniQM)
01090 {
01091 gamPDG=111;
01092 gamM=mPi0;
01093 }
01094 G4LorentzVector r4Mom(0.,0.,0.,gamM);
01095 G4LorentzVector s4Mom(0.,0.,0.,iniQM);
01096 G4double sum=gamM+iniQM;
01097 if(sum>0. && fabs(quasM-sum)<eps)
01098 {
01099 r4Mom=q4Mom*(gamM/sum);
01100 s4Mom=q4Mom*(iniQM/sum);
01101 }
01102 else if(quasM<sum || !G4QHadron(q4Mom).DecayIn2(r4Mom, s4Mom))
01103 {
01104
01105
01106
01107 G4ExceptionDescription ed;
01108 ed << "(E=0)G/Pi0+GSQ decay error: Q=" << q4Mom << quasM
01109 << "->g/pi0(M=" << gamM << ")+GSQ=" << iniPDG << "(M="
01110 << iniQM << ")=" << sum << ", d=" << sum-quasM << G4endl;
01111 G4Exception("G4Quasmon::HadronizeQuasmon()", "HAD_CHPS_0002",
01112 FatalException, ed);
01113 }
01114 #ifdef debug
01115 G4cout<<"G4Q::HQ:=== 0 ===>HadrVec, Q="<<q4Mom<<quasM<<"->g/pi0("<<gamPDG<<")="
01116 <<r4Mom<<gamM<<"+GSQ="<<iniPDG<<r4Mom<<iniQM<<G4endl;
01117 #endif
01118 G4QHadron* curHadr2 = new G4QHadron(gamPDG,r4Mom);
01119 FillHadronVector(curHadr2);
01120 G4QHadron* curHadr1 = new G4QHadron(iniPDG,s4Mom);
01121 FillHadronVector(curHadr1);
01122 ClearQuasmon();
01123 }
01124 #ifdef debug
01125 G4cout<<"***G4Q::HQ:dE="<<excE<<">minK="<<minK<<",Env="<<theEnvironment<<",k="<<kLS
01126 <<",Q="<<valQ<<quasM<<", nQ="<<nQuasms<<G4endl;
01127 #endif
01128 qEnv=theEnvironment;
01129 return theQHadrons;
01130 }
01131 }
01132 G4bool nucflag=false;
01133 G4bool hsflag=false;
01134
01135 G4LorentzVector rQ4Mom(0.,0.,0.,0.);
01136 G4LorentzVector fQ4Mom(0.,0.,0.,0.);
01137 G4LorentzVector fr4Mom(0.,0.,0.,0.);
01138 G4double rMass=0.;
01139 G4double kt=0.;
01140
01141 G4double ku=0.;
01142
01143 G4double kn=0.;
01144 G4double sCBE=0.;
01145 G4int rPDG=0;
01146 pCond=true;
01147 G4int pCount=0;
01148 G4LorentzVector PMEMfr4M(0.,0.,0.,0.);
01149 G4LorentzVector PMEMrQ4M(0.,0.,0.,0.);
01150 G4QContent PMEMpQC(0,0,0,0,0,0);
01151 G4QContent PMEMsQC(0,0,0,0,0,0);
01152 G4QContent PMEMtQC(0,0,0,0,0,0);
01153 G4QContent PMEMcQC(0,0,0,0,0,0);
01154 G4double PMEMktM2=0.;
01155 G4double PMEMknM2=0.;
01156 G4double PMEMreM2=0.;
01157 G4double PMEMrMas=0.;
01158 G4double PMEMpMas=0.;
01159 G4double PMEMsMas=0.;
01160 G4double PMEMdMas=0.;
01161 G4double PMEMmiSN=0.;
01162 G4double PMEMmiST=0.;
01163 G4double PMEMmiSB=0.;
01164 G4int PMEMrPDG=0;
01165 G4int PMEMsPDG=0;
01166 G4int PMEMpPDG=0;
01167 G4bool PMEMhsfl=hsflag;
01168 G4bool PMEMnucf=nucflag;
01169 #ifdef debug
01170 G4cout<<"G4Q::HQ: fp="<<fprob<<",QM="<<quasM<<",QQC="<<valQ<<",k="<<kMom<<G4endl;
01171 #endif
01172 while(pCount<pCountMax&&pCond)
01173 {
01174 hsflag=false;
01175 #ifdef debug
01176 G4cout<<"G4Q::HQ:***>New p-Attempt#"<<pCount<<",pMax="<<pCountMax<<",hsfl=0"<<G4endl;
01177 #endif
01178 G4double totP = maxP * G4UniformRand();
01179 while(theQCandidates[i]->GetIntegProbability() < totP) i++;
01180 if (i>=nCandid)
01181 {
01182 G4ExceptionDescription ed;
01183 ed <<"Too big number of the candidate: Cand#"<< i <<" >= Tot#"<< nCandid << G4endl;
01184 G4Exception("G4Quasmon::HadronizeQuasmon()", "HAD_CHPS_0003", FatalException, ed);
01185 }
01186 curQ = valQ;
01187 if (!fprob)
01188 {
01189 memQ=curQ;
01190 sPDG=curQ.GetSPDGCode();
01191 if(!sPDG&&theEnvironment.GetPDG()!=NUCPDG&&totBN>1&&totMass>totM&&totS>=0)
01192 {
01193 #ifdef pdebug
01194 G4ExceptionDescription ed;
01195 ed << "Why Fail? (1): NEED-EVAP-1:Q=" << q4Mom << valQ << ",En="
01196 << theEnvironment << G4endl;
01197 G4Exception("G4Quasmon::HadronizeQuasmon()","HAD_CHPS_0004", FatalException, ed);
01198 #endif
01199 qEnv=theEnvironment;
01200 return theQHadrons;
01201 }
01202 else if(!sPDG)
01203 {
01204 G4ExceptionDescription ed;
01205 ed << "DecayPartSelection,Evaporation: sPDG=0,E=" << theEnvironment
01206 << ",B=" << totBN << ",S=" << totS << G4endl;
01207 G4Exception("G4Quasmon::HadronizeQuasmon()","HAD_CHPS_0005", FatalException, ed);
01208 }
01209 else if(sPDG==10)
01210 {
01211
01212 G4QChipolino chipQ(valQ);
01213 G4QPDGCode QPDG1=chipQ.GetQPDG1();
01214 sPDG = QPDG1.GetPDGCode();
01215 sMass= QPDG1.GetMass();
01216 G4QPDGCode QPDG2=chipQ.GetQPDG2();
01217 rPDG = QPDG2.GetPDGCode();
01218 rMass= QPDG2.GetMass();
01219 if(sMass+rMass>quasM)
01220 {
01221 if(totBN>1&&totMass>totM&&totS>=0)
01222 {
01223 #ifdef pdebug
01224 G4ExceptionDescription ed;
01225 ed << " Why Fail? (2): NEED-EVAP-2:Q=" << q4Mom << valQ << ",E="
01226 << theEnvironment << G4endl;
01227 G4Exception("G4Quasmon::HadronizeQuasmon()", "HAD_CHPS_0006",
01228 FatalException, ed);
01229 #endif
01230 qEnv=theEnvironment;
01231 return theQHadrons;
01232 }
01233 else
01234 {
01235 G4ExceptionDescription ed;
01236 ed << "VirtChipo Can't EvapNucleus: QM=" << quasM << "<S=" << sMass
01237 << "+R=" << rMass << "=" << sMass+rMass << ",tB=" << totBN
01238 << ",tS=" << totS << ",tM=" << totMass << ">minM=" << totM << G4endl;
01239 G4Exception("G4Quasmon::HadronizeQuasmon()", "HAD_CHPS_0007",
01240 FatalException, ed);
01241 }
01242 }
01243 }
01244 else if(nQuasms>1)
01245 {
01246 #ifdef pdebug
01247 G4ExceptionDescription ed;
01248 ed << "Why Fail ? (0): NEED-EVAP-0:Q=" << q4Mom << valQ << ",En="
01249 << theEnvironment << G4endl;
01250 G4Exception("G4Quasmon::HadronizeQuasmon()","HAD_CHPS_0008", FatalException, ed);
01251 #endif
01252 qEnv=theEnvironment;
01253 return theQHadrons;
01254 }
01255 else
01256 {
01257 sMass=G4QPDGCode(sPDG).GetMass();
01258 rPDG=envPDG;
01259 if (rPDG>MINPDG&&rPDG!=NUCPDG)
01260 {
01261 rMass=theEnvironment.GetMZNS();
01262 q4Mom+=G4LorentzVector(0.,0.,0.,rMass);
01263 valQ +=theEnvironment.GetQC();
01264 quasM=q4Mom.m();
01265 KillEnvironment();
01266
01267 }
01268 else if(rPDG!=NUCPDG&&totBN>1&&totMass>totM&&totS>=0)
01269 {
01270 #ifdef pdebug
01271 G4QContent nTotQC=totQC-neutQC;
01272 G4QNucleus nTotN(nTotQC);
01273 G4double nTotM =nTotN.GetMZNS();
01274 G4double dMnT=totMass-nTotM-mNeut;
01275 G4QContent pTotQC=totQC-protQC;
01276 G4QNucleus pTotN(pTotQC);
01277 G4double pTotM =pTotN.GetMZNS();
01278 G4double dMpT=totMass-pTotM-mProt;
01279 if(dMpT>dMnT)dMnT=dMpT;
01280 G4ExceptionDescription ed;
01281 ed << "Why Fail ? (3): NEED-EVAP3:s=" << sPDG << ",Q=" << q4Mom
01282 << valQ << ",r=" << rPDG << G4endl;
01283 G4Exception("G4Quasmon::HadronizeQuasmon()","HAD_CHPS_0009",FatalException,ed);
01284 #endif
01285 qEnv=theEnvironment;
01286 return theQHadrons;
01287 }
01288 else if(rPDG==NUCPDG)
01289 {
01290 #ifdef debug
01291 G4cout<<"G4Quasmon::HadronizeQuasm:SafatyDecayIn PI0/GAM, rPDG="<<rPDG<<G4endl;
01292 #endif
01293 if(totMass-totM>mPi0)
01294 {
01295 rMass=mPi0;
01296 rPDG=111;
01297 }
01298 else
01299 {
01300 rMass=0.;
01301 rPDG=22;
01302 }
01303 }
01304 }
01305 hsflag=true;
01306 #ifdef debug
01307 G4cout<<"G4Q::HQ:hsflagTRUE,s="<<sPDG<<","<<sMass<<",r="<<rPDG<<","<<rMass<<G4endl;
01308 #endif
01309 }
01310 else
01311 {
01312 G4QCandidate* curCand = theQCandidates[i];
01313 sPDG = curCand->GetPDGCode();
01315 #ifdef debug
01316 G4cout<<"G4Q::HQ:hsfl="<<hsflag<<", sPDG="<<sPDG<<", i="<<i<<G4endl;
01317 #endif
01318
01319 if(sPDG>MINPDG&&sPDG!=NUCPDG)
01320 {
01321 G4int ip=0;
01322 G4int nParCandid = curCand->GetPClustEntries();
01323 G4double sppm = curCand->TakeParClust(nParCandid-1)->GetProbability();
01324 if (sppm<=0)
01325 {
01326 G4cerr<<"***G4Quasmon::HadronizeQ:P="<<theQCandidates[i]->GetIntegProbability()
01327 <<",nC="<<nParCandid<<",pP="<<sppm<<",QM="<<quasM<<",QC="<<valQ;
01328 for(int ipp=0; ipp<nParCandid; ipp++)
01329 G4cerr<<", "<<ipp<<": "<<curCand->TakeParClust(ip)->GetProbability();
01330 G4Exception("G4Quasmon::HadronizeQuasmon()", "HAD_CHPS_0010", FatalException,
01331 "NoParentClust forTheFragment");
01332 }
01333 else
01334 {
01335 G4double totPP = sppm * G4UniformRand();
01336 while(curCand->TakeParClust(ip)->GetProbability() < totPP) ip++;
01337 #ifdef debug
01338 G4cout<<"G4Q::HQ:p#ip="<<ip<<",f#i="<<i<<",tP="<<totPP<<",sP="<<sppm<<G4endl;
01339 #endif
01340 }
01341 G4QParentCluster* parCluster=curCand->TakeParClust(ip);
01342 pPDG = parCluster->GetPDGCode();
01343 G4QPDGCode pQPDG(pPDG);
01344 pQC = pQPDG.GetQuarkContent();
01345 pBaryn= pQC.GetBaryonNumber();
01346 pMass = parCluster->GetEBMass();
01347 pNMass = parCluster->GetNBMass();
01348 transQC = parCluster->GetTransQC();
01349 delta = parCluster->GetEBind();
01350 deltaN = parCluster->GetNBind();
01351 loli = parCluster->GetLow();
01352 hili = parCluster->GetHigh();
01353
01354
01355
01356 npqp2 = parCluster->GetNQPart2();
01357
01358 G4QPDGCode sQPDG(curCand->GetPDGCode());
01359 sQC = sQPDG.GetQuarkContent();
01360
01361 if(sPDG==90001001 && G4UniformRand()>1.0) sMass=np;
01362 else sMass = sQPDG.GetMass();
01363 sM2 = sMass*sMass;
01364 curQ += transQC;
01365 #ifdef debug
01366 G4cout<<"G4Q::HQ:valQ="<<valQ<<"+transQ="<<transQC<<"("<<pPDG<<" to "<<sPDG
01367 <<") = curQ="<<curQ<<",InvBinding="<<delta<<",pM="<<pMass<<pQC<<G4endl;
01368 #endif
01369 }
01370 else
01371 {
01372 pBaryn=0;
01373 sQC=theQCandidates[i]->GetQC();
01374 sMass = theQCandidates[i]->GetNBMass();
01375 sM2=sMass*sMass;
01376 curQ-= sQC;
01377 #ifdef debug
01378 G4cout<<"G4Q::HQ: hsfl="<<hsflag<<", valQ="<<valQ<<"-sQ="<<sQC<<",sM="<<sMass
01379 <<",C="<<theQCandidates[i]->GetPDGCode()<<",Q="<<curQ<<",M2="<<sM2<<G4endl;
01380 #endif
01381 }
01382 G4QContent resNQC=totQC-sQC;
01383 G4QNucleus resTN(resNQC);
01384 G4double resTNM=resTN.GetMZNS();
01385 G4double sCB=0;
01386 if(resTN.GetA()>0) sCB=totN.CoulombBarrier(sQC.GetCharge(),sQC.GetBaryonNumber());
01387 sCBE=sCB+sMass;
01388 #ifdef debug
01389 G4cout<<"G4Q::HQ:rQC="<<resNQC<<",rM="<<resTNM<<",sM="<<sMass<<",CB="<<sCB<<G4endl;
01390 #endif
01391 memQ=curQ;
01392 G4double rtM=0.;
01393 G4double reM=0.;
01394
01395
01396 if(envA>=pBaryn)
01397 {
01398
01399 G4QContent RNQC=curQ+envQC;
01400 if(sPDG>MINPDG&&sPDG!=NUCPDG) RNQC-=pQC;
01401 if(envA-pBaryn>bEn) RNQC=curQ+bEnQC;
01402 G4int RNPDG = RNQC.GetSPDGCode();
01403 if(RNPDG==10) minSqN=G4QChipolino(RNQC).GetMass2();
01404 else if(!RNPDG)
01405 {
01406
01407 G4cout<<"**G4Q::HQ:*KinematicTotal*, PDG="<<RNPDG<<curQ<<", QC="<<RNQC<<G4endl;
01408
01409 minSqN=1000000.;
01410 }
01411 else
01412 {
01413 G4double minN=G4QPDGCode(RNPDG).GetMass();
01414 minSqN=minN*minN;
01415 #ifdef debug
01416 G4cout<<"G4Q::HQ:M="<<bEn<<",A="<<envA<<",B="<<pBaryn<<",N="<<minN<<G4endl;
01417 #endif
01418 }
01419 }
01420 else
01421 {
01422 #ifdef debug
01423 G4cout<<"*G4Q::HQ:EnvironmentA="<<envA<<" < SecondaryFragmentA="<<pBaryn<<G4endl;
01424 #endif
01425 }
01426
01427 G4int rqPDG = curQ.GetSPDGCode();
01428 if(rqPDG==111&&sPDG!=111&&G4UniformRand()>.5) rqPDG=221;
01429
01430 G4int rQQ=G4QPDGCode(curQ).GetQCode();
01431 if(rqPDG==10) {
01432 minSqT=G4QChipolino(curQ).GetMass2();
01433 minSqB=minSqT;
01434 minT=sqrt(minSqT);
01435
01436 } else if(!rqPDG||rQQ<-1) {
01437 #ifdef debug
01438 G4cerr<<"*G4Q::HQ:*** ResidualQuasmon *** PDG="<<rqPDG<<curQ<<",Q="<<rQQ<<G4endl;
01439 #endif
01440 minT=100000.;
01441 minSqT=10000000000.;
01442 minSqB=10000000000.;
01443
01444 } else {
01446 minT=G4QPDGCode(rqPDG).GetMass();
01447 if(sPDG<MINPDG&&envPDG>MINPDG&&envPDG!=NUCPDG)
01448 {
01449 G4int rqZ=curQ.GetCharge();
01450 G4int rqS=curQ.GetStrangeness();
01451 G4int rqN=curQ.GetBaryonNumber()-rqS-rqZ;
01452 G4double qpeM=G4QNucleus(envZ+rqZ,envN+rqN,envS+rqS).GetGSMass();
01453 #ifdef debug
01454 G4cout<<"G4Q::HadQ:Z="<<rqZ<<",N="<<rqN<<",S="<<rqS<<",eZ="<<envZ<<",eN="<<envN
01455 <<",eS="<<envS<<",ePDG="<<envPDG<<",eM="<<envM<<",tM="<<qpeM<<G4endl;
01456 #endif
01457 minT=qpeM-envM;
01458 }
01459 minSqT=minT*minT;
01460 #ifdef debug
01461 G4cout<<"G4Q::HQ:rPDG="<<rqPDG<<curQ<<",minT="<<minT<<",minSqT="<<minSqT
01462 <<",hsfl="<<hsflag<<G4endl;
01463 #endif
01464 G4double newT=0.;
01465
01466
01467
01468 if ( ( (sPDG < MINPDG && envPDG > MINPDG && envPDG != NUCPDG) ||
01469 (sPDG > MINPDG && sPDG != NUCPDG && envPDG > pPDG)
01470 ) && ( (rqPDG > MINPDG && rqPDG != NUCPDG) ||
01471 rqPDG==2112 || rqPDG==2212 || rqPDG==3122
01472 )
01473 )
01474 {
01475 if(sPDG<MINPDG)
01476 {
01477
01478 G4QContent rtQC=curQ;
01479 if (envA > bEn)
01480 {
01481 reM=mbEn;
01482 rtQC+=bEnQC;
01483 }
01484 else
01485 {
01486 reM=envM;
01487 rtQC+=envQC;
01488 }
01489 G4QNucleus rtN(rtQC);
01490 rtM=rtN.GetMZNS();
01491 newT=rtM-reM;
01492 #ifdef debug
01493 G4cout<<"G4Q::HQ:***VacuumFragmentation** M="<<newT<<",rM="<<rtM<<rtQC
01494 <<",eM="<<envM<<",mM="<<minT<<G4endl;
01495 #endif
01496 }
01497 else
01498 {
01499
01500 G4QContent reQC=envQC-pQC;
01501 if(envA-pBaryn>bEn) reQC=bEnQC;
01502 G4QNucleus reN(reQC);
01503 reM=reN.GetMZNS();
01504
01505 G4QContent rtQC=curQ;
01506 #ifdef debug
01507 G4cout<<"G4Q::HQ:reQC="<<reQC<<",rtQC="<<rtQC<<",eA="<<envA<<",pB="<<pBaryn
01508 <<",bE="<<bEn<<bEnQC<<G4endl;
01509 #endif
01510 rtQC+=reQC;
01511 G4QNucleus rtN(rtQC);
01512 rtM=rtN.GetMZNS();
01513
01514 if (envA-pBaryn > bEn) newT=rtM-mbEn;
01515 else newT=rtM-reM;
01516 #ifdef debug
01517 G4cout<<"G4Q::HQ:NuclFrM="<<newT<<",r="<<rtM<<rtQC<<",e="<<envM<<envQC<<",p="
01518 <<pMass<<pQC<<",re="<<reM<<reQC<<",exEn="<<totMass-rtM-sMass<<G4endl;
01519 #endif
01520 }
01521 if(minT<newT) newT=minT;
01522 }
01523 minSqB=newT*newT;
01524 }
01525 #ifdef debug
01526 G4cout<<"G4Q::HQ:rq="<<rqPDG<<",miT="<<minSqT<<",miB="<<minSqB<<",M="<<rtM<<G4endl;
01527 #endif
01528 if(!minSqT)
01529 {
01530 G4ExceptionDescription ed;
01531 ed << "MinResMass can't be calculated: minSqT=0(!), curQ=" << curQ << G4endl;
01532 G4Exception("G4Quasmon::HadronizeQuasmon()","HAD_CHPS_0011", FatalException, ed);
01533 }
01534 G4double m2_value = BIG2;
01535
01536 if (sPDG > MINPDG && sPDG != NUCPDG) {
01537 #ifdef debug
01538 G4cout<<"G4Q::HQ: BoundM="<<pMass<<",FreeM="<<sMass<<",QM="<<quasM<<G4endl;
01539 #endif
01540
01541
01542 G4LorentzVector cl4Mom(0.,0.,0.,pMass);
01543 G4LorentzVector tot4Mom=q4Mom+cl4Mom;
01544 #ifdef debug
01545 G4cout<<"G4Q::HQ:Q("<<quasM<<")->k("<<k4Mom<<")+CRQ("<<cr4Mom.m()<<")"<<G4endl;
01546 #endif
01547 G4LorentzVector cc4Mom=k4Mom+cl4Mom;
01548 G4double ccM2=cc4Mom.m2();
01549 G4double frM2=sMass*sMass;
01550 if (ccM2 <= frM2)
01551 {
01552 #ifdef debug
01553 G4cout<<"***G4Q::HQ:FailedToFind FragmM:"<<ccM2<<"<"<<frM2<<",M="<<pMass<<"+k="
01554 <<k4Mom<<"="<<sqrt(ccM2)<<cc4Mom<<" < fM="<<sMass<<",miK="<<minK<<G4endl;
01555 #endif
01556 dMass=pMass-pNMass;
01557 pMass=pNMass;
01558 delta=deltaN;
01559 cl4Mom=G4LorentzVector(0.,0.,0.,pMass);
01560 tot4Mom=q4Mom+cl4Mom;
01561 cc4Mom=k4Mom+cl4Mom;
01562 ccM2=cc4Mom.m2();
01563 if (ccM2 <= frM2)
01564 {
01565 #ifdef debug
01566 G4cout<<"G4Q::HQ:hsflagTRUE*NuclBINDING,ccM2="<<ccM2<<"<frM2="<<frM2<<G4endl;
01567 #endif
01568 hsflag=true;
01569 }
01570 else
01571 {
01572 #ifdef debug
01573 G4cout<<"G4Q::HQ:***NUCLEAR BINDING***ccM2="<<ccM2<<" > frM2="<<frM2<<G4endl;
01574 #endif
01575 nucflag=true;
01576 }
01577 }
01578 else
01579 {
01580 #ifdef debug
01581 G4double crMass2 = cr4Mom.m2();
01582 G4cout<<"G4Q::HQ:cM2="<<crMass2<<"="<<rEP*(rEP-rMo-rMo)<<",h="<<hili<<",l="
01583 <<loli<<G4endl;
01584 #endif
01586 if(hili<loli) hili=loli;
01587 G4double fpqp2=static_cast<double>(npqp2);
01588 G4double pw=1./fpqp2;
01589
01590 qCond=true;
01591 G4int qCount=0;
01592 #ifdef pdebug
01593 G4double dM=0.;
01594 if(sPDG==90001001) dM=2.25;
01595 G4cout<<"G4Q::HQ:Is xE="<<excE<<" > sM="<<sMass<<"-pM="<<pMass<<"-dM="<<dM
01596 <<" = "<<sMass-pMass-dM<<G4endl;
01597 #endif
01598 #ifdef debug
01599 G4cout<<"G4Q::HQ: must totM="<<totMass<<" > rTM="<<resTNM<<"+sM="<<sMass<<" = "
01600 <<sMass+resTNM<<G4endl;
01601 #endif
01602 if(resTNM && totMass<resTNM+sMass)
01603 {
01604 #ifdef pdebug
01605 G4cout<<"***G4Quasmon::HadronizeQuasmon:***PANIC#1***TotalDE="<<excE<<"< bE="
01606 <<sMass-pMass-dM<<", dM="<<dM<<", sM="<<sMass<<", bM="<<pMass<<G4endl;
01607
01608 #endif
01609 status =-1;
01610 qEnv=theEnvironment;
01611 return theQHadrons;
01612 }
01613 G4double ex=kLS-delta;
01614 G4double dex=ex+ex;
01615 G4QContent tmpEQ=envQC-pQC;
01616 if(envA-pBaryn>bEn) tmpEQ=bEnQC;
01617 G4QNucleus tmpN(tmpEQ);
01618 G4double tmpNM=tmpN.GetMZNS();
01619 #ifdef debug
01620 G4cout<<"G4Q::HQ:eQC="<<envQC<<",pQC="<<pQC<<",rEnvM="<<tmpNM<<",hsfl="<<hsflag
01621 <<G4endl;
01622 #endif
01623 G4QContent tmpRQ=valQ+transQC;
01624 G4QContent tmpTQ=tmpRQ+tmpEQ;
01625 G4QNucleus tmpT(tmpTQ);
01626 G4double tmpTM=tmpT.GetMZNS();
01627 tmpTM2=tmpTM*tmpTM;
01628 G4LorentzVector ResEnv4Mom(0.,0.,0.,tmpNM);
01629 G4LorentzVector tCRN=cr4Mom+ResEnv4Mom;
01630
01631 G4double tmpBE=minT+tmpNM-tmpTM;
01632 #ifdef debug
01633
01634 G4double tcEP=tCRN.e()+tCRN.rho();
01635
01636 G4double cta=1.-(dex/(1.-pow(loli,pw))-pMass)/kLS;
01637 if(cta>1.0001)G4cerr<<"Warn-G4Q::HQ: cost_max="<<cta<<">1.CorHadrProb"<<G4endl;
01638 G4double kap_a=ex/(1.+kLS*(1.-cta)/pMass);
01639 G4double cti=1.-(dex/(1.-pow(hili,pw))-pMass)/kLS;
01640 if(cti<-1.0001)G4cerr<<"Warn-G4Q::HQ: cost_i="<<cti<<"<-1.CorHadrProb"<<G4endl;
01641 G4double kap_i=ex/(1.+kLS*(1.-cti)/pMass);
01642 G4double q_lim=(tmpTM2-tCRN.m2())/(tcEP+tcEP);
01643 if(cti>cta+.0001)G4cerr<<"**G4Q::HQ:ci="<<cti<<">ca="<<cta<<".CorHPro"<<G4endl;
01644 G4cout<<"G4Q::HQ:qi="<<kap_i<<",ci="<<cti<<",a="<<kap_a<<",ca="<<cta<<",e="<<ex
01645 <<",q="<<q_lim<<",S="<<tmpTM*tmpTM<<",R2="<<tCRN.m2()<<","<<tcEP<<G4endl;
01646 #endif
01647
01648 #ifdef debug
01649
01650 #endif
01651
01652
01653 #ifdef debug
01654
01655 #endif
01656
01657 G4LorentzVector MEMkp4M(0.,0.,0.,0.);
01658 G4LorentzVector MEMfr4M(0.,0.,0.,0.);
01659 G4LorentzVector MEMrQ4M(0.,0.,0.,0.);
01660 G4double MEMrQM2=0.;
01661 G4double MEMsCBE=0.;
01662 G4double MEMreM2=0.;
01663
01664
01665
01666
01667
01668 while(qCount<qCountMax&&qCond)
01669 {
01670 G4double z = pow(loli+(hili-loli)*G4UniformRand(),pw);
01671
01672
01673
01674 G4double ctkk=1.-(dex/(1.-z)-pMass)/kLS;
01675 #ifdef pdebug
01676 if(qCount) G4cout<<"G4Q::HQ:qC="<<qCount<<",ct="<<ctkk<<",M="<<pMass<<",z="
01677 <<z<<",zl="<<pow(loli,pw)<<",zh="<<pow(hili,pw)<<",dE="
01678 <<totMass-totM<<",bE="<<sMass-pMass<<G4endl;
01679 #endif
01680 #ifdef debug
01681 G4cout<<"G4Q::HQ:ct="<<ctkk<<",pM="<<pMass<<",z="<<z<<",zl="<<pow(loli,pw)
01682 <<",zh="<<pow(hili,pw)<<",ex="<<ex<<",li="<<loli<<",hi="<<hili<<G4endl;
01683 #endif
01684 if(abs(ctkk)>1.00001)
01685 {
01686 #ifdef debug
01687 G4cerr<<"***G4Q:HQ:ctkk="<<ctkk<<",ex="<<ex<<",z="<<z<<",pM="<<pMass
01688 <<",kLS="<<kLS<<",hi="<<hili<<",lo="<<loli<<",n="<<npqp2<<G4endl;
01689
01690 #endif
01691 if(ctkk> 1.)ctkk= 1.;
01692 if(ctkk<-1.)ctkk=-1.;
01693 }
01694 G4double cen=kLS+pMass;
01695 G4double ctc=(cen*ctkk-kLS)/(cen-kLS*ctkk);
01696 if(abs(ctc)>1.)
01697 {
01698
01699 if(ctc>1.) ctc=1.;
01700 else if(ctc<-1.) ctc=-1.;
01701 }
01702 kp4Mom=zeroLV;
01703 fr4Mom=G4LorentzVector(0.,0.,0.,sMass);
01704 if(!G4QHadron(cc4Mom).RelDecayIn2(kp4Mom, fr4Mom, k4Mom, ctc, ctc))
01705 {
01706 G4ExceptionDescription ed;
01707 ed << "Can't dec ColClust(Fr+kap): c4M=" << cc4Mom << ",sM="
01708 << sMass << ",ct=" << ctc << G4endl;
01709 G4Exception("G4Quasmon::HadronizeQuasmon()", "HAD_CHPS_0012",
01710 FatalException, ed);
01711 }
01712
01713
01714
01715
01716
01717 #ifdef debug
01718 G4double ccM=sqrt(ccM2);
01719 G4double kappa=ex/(1.+kLS*(1.-ctkk)/pMass);
01720 G4cout<<"G4Q::HQ:>>ColDec>>CF("<<ccM<<")->F("<<sMass<<")+q"<<kp4Mom<<"="
01721 <<kappa<<",hsfl="<<hsflag<<G4endl;
01722 #endif
01723
01724
01725 rQ4Mom=cr4Mom+kp4Mom;
01726 G4LorentzVector retN4Mom=rQ4Mom+ResEnv4Mom;
01727 reTNM2=retN4Mom.m2();
01728
01729
01730
01731
01732
01733
01734
01735
01736
01737
01738 fQ4Mom=rQ4Mom+G4LorentzVector(0.,0.,0.,tmpBE);
01739 G4double fQM2=fQ4Mom.m2();
01740 if(fQM2>=minSqT && reTNM2>=tmpTM2 && fr4Mom.e()>=sCBE)
01741 {
01742 qCond=false;
01743
01744 #ifdef debug
01745
01746
01747
01748
01749
01750
01751 G4cout<<"G4Q::HQ:Attemp#"<<qCount<<".Yes.M2="<<fQM2<<">"<<minSqT<<G4endl;
01752 #endif
01753 }
01754 else
01755 {
01756 #ifdef debug
01757
01758
01759
01760
01761
01762
01763 G4cout<<"G4Q::HQ:Attempt#"<<qCount<<",NO.M2="<<fQM2<<"<"<<minSqT<<G4endl;
01764
01765
01766 #endif
01767
01768
01769
01770
01771
01772 if(reTNM2<tmpTM2 && fQM2>MEMrQM2 && fr4Mom.e()>=sCBE)
01773 {
01774
01775 MEMrQM2=fQM2;
01776
01777
01778
01779
01780
01781 MEMkp4M=kp4Mom;
01782 MEMfr4M=fr4Mom;
01783 MEMrQ4M=rQ4Mom;
01784 MEMreM2=reTNM2;
01785 }
01786
01787 else if(fr4Mom.e()<sCBE&&fr4Mom.e()>MEMsCBE&&reTNM2>=tmpTM2&&fQM2>MEMrQM2)
01788
01789
01790
01791
01792 {
01793 MEMsCBE=fr4Mom.e();
01794 MEMkp4M=kp4Mom;
01795 MEMfr4M=fr4Mom;
01796 MEMrQ4M=rQ4Mom;
01797 MEMreM2=reTNM2;
01798 }
01799 else if(!qCount)
01800 {
01801
01802 MEMrQM2=fQM2;
01803
01804
01805
01806
01807
01808 MEMsCBE=fr4Mom.e();
01809 MEMkp4M=kp4Mom;
01810 MEMfr4M=fr4Mom;
01811 MEMrQ4M=rQ4Mom;
01812 MEMreM2=reTNM2;
01813 }
01814 else
01815 {
01816
01817 fQM2=MEMrQM2;
01818
01819
01820
01821
01822
01823 kp4Mom=MEMkp4M;
01824 fr4Mom=MEMfr4M;
01825 rQ4Mom=MEMrQ4M;
01826 reTNM2=MEMreM2;
01827 }
01828 }
01829 qCount++;
01830 }
01831
01832
01833 #ifdef debug
01834 G4cout<<"G4Q::HadQ:RQ("<<rQ4Mom.m()<<")=C("<<cr4Mom.m()<<")+q"<<kp4Mom<<G4endl;
01835 #endif
01836
01837
01838 ku=fQ4Mom.m2();
01839
01840
01841
01842 G4LorentzVector totC4Mom=rQ4Mom;
01843 if(envA-pBaryn>bEn) totC4Mom+=G4LorentzVector(0.,0.,0.,mbEn);
01844 else totC4Mom+=G4LorentzVector(0.,0.,0.,envM-pMass);
01845 kn=totC4Mom.m2();
01846 #ifdef debug
01847 G4cout<<"G4Q::HQ:A="<<envA<<",B="<<pBaryn<<",Q="<<rQ4Mom.m()<<","<<piF<<G4endl;
01848 #endif
01849 m2_value=kt;
01850 tot4Mom-=rQ4Mom+fr4Mom;
01851 #ifdef debug
01852 G4cout<<"G4Q::HQ:t4M="<<tot4Mom<<",hsfl="<<hsflag<<".Is kt="<<kt<<">"<<minSqB
01853 <<" or kn="<<kn<<">"<<minSqN<<"? m2="<<m2_value<<", sPDG="<<sPDG<<G4endl;
01854 #endif
01855
01856
01857
01858
01859 if(kn<minSqN && ku<minSqT)
01860
01861
01862
01863
01864
01865
01866
01867
01868
01869
01870
01871
01872
01873 {
01874 hsflag=true;
01875 #ifdef debug
01876 G4cout<<"G4Q::HQ:**hsflag=1** No, sPDG="<<sPDG<<", kt="<<kt<<"<"<<minSqB
01877 <<" or kn="<<kn<<"<"<<minSqN<<G4endl;
01878 #endif
01879 }
01880 #ifdef debug
01881 else G4cout<<"G4Q::HQ:YES,t="<<kt<<">"<<minSqB<<",n="<<kn<<">"<<minSqN<<G4endl;
01882 #endif
01883 }
01884 }
01885 else
01886 {
01887 kt = (quasM-dk)*(quasM-sM2/dk);
01888 G4double rQM=0.;
01889 if(kt>0.) rQM=sqrt(kt);
01890 fr4Mom=G4LorentzVector(0.,0.,0.,sMass);
01891 rQ4Mom=G4LorentzVector(0.,0.,0.,rQM);
01892 if(!G4QHadron(q4Mom).DecayIn2(fr4Mom, rQ4Mom))
01893 {
01894 G4ExceptionDescription ed;
01895 ed << "Can't dec Quasmon in Fr+rQuas: q4M=" << q4Mom << ", sM="
01896 << sMass << ", rQM=" << rQM << G4endl;
01897 G4Exception("G4Quasmon::HadronizeQuasmon()","HAD_CHPS_0013",FatalException,ed);
01898 }
01899
01900 if(envPDG>MINPDG&&envPDG!=NUCPDG)
01901 {
01902
01903 G4LorentzVector TCRN=rQ4Mom;
01904 if(envA>bEn) TCRN+=bEn4M;
01905 else TCRN+=env4M;
01906 kn=TCRN.m2();
01907 }
01908 else kn=kt;
01909 }
01910
01911 G4LorentzVector tL=rQ4Mom;
01912 tL+=G4LorentzVector(0.,0.,0.,reM);
01913 G4double tM=tL.m();
01914
01915 #ifdef debug
01916 G4cout<<"G4Q::HQ:k="<<kMom<<".F:"<<kt<<">"<<minSqB<<",N:"<<kn<<">"<<minSqN<<" &tM="
01917 <<tM<<">rtM="<<rtM<<" & hsfl="<<hsflag<<" to avoid decay R+S="<<sPDG<<G4endl;
01918 #endif
01919
01920
01921
01922
01923
01924 if(kt>minSqB+.01 && tM>rtM && !hsflag)
01925
01926
01927
01928
01929
01930
01931 {
01932
01933
01934
01935
01936
01937
01938
01939
01940
01941 #ifdef debug
01942
01943
01944 #endif
01945
01946
01947
01948
01949 #ifdef debug
01950 G4cout<<"G4Q::HQ:YES forFragment it's big enough:kn="<<kn<<">"<<minSqN<<G4endl;
01951 #endif
01952 m2_value = kt;
01953
01954 }
01955 else
01956 {
01957 hsflag=true;
01958 #ifdef debug
01959 G4cout<<"G4Q::HQ:NO,hsfl=1,kt="<<kt<<"<"<<minSqB<<" or M="<<tM<<"<"<<rtM<<G4endl;
01960 #endif
01961 }
01962 #ifdef debug
01963 G4cout<<"G4Q::HQ:******>>rM="<<rMass<<",sqM2="<<sqrt(m2_value)<<",hsfl="<<hsflag<<G4endl;
01964 #endif
01965 rMass=sqrt(m2_value);
01966 G4double m3_value=m2_value*rMass;
01967 G4int cB = abs(curQ.GetBaryonNumber());
01969 rPDG = curQ.GetSPDGCode();
01970 G4double rrn=G4UniformRand();
01971 if(rPDG==111&&sPDG!=111&&rrn>.5) rPDG=221;
01972
01973 G4int aPDG = abs(rPDG);
01974 G4int rb = abs(curQ.GetBaryonNumber());
01975 G4double rcMass=-BIG;
01976 if (!rPDG)
01977 {
01978 G4ExceptionDescription ed;
01979 ed << "Unidentifiable residual Hadron: rQ =" << curQ << ", rPDG=" << rPDG
01980 << "(b=" << rb << ") + sPDG=" << sPDG << "(sM=" << sMass << ")" << G4endl;
01981 G4Exception("G4Quasmon::HadronizeQuasmon()", "HAD_CHPS_0014",
01982 FatalException, ed);
01983 }
01984 G4double sB = 1473.;
01985 if(rPDG!=10) rcMass=G4QPDGCode(rPDG).GetMass();
01986 else sB=0.;
01987 if(rPDG==221 || rPDG==331) rcMass=mPi0;
01988 G4double bs = rcMass+mPi0;
01989 G4bool rFl=false;
01990
01991
01992 if(sPDG<MINPDG&&envPDG==NUCPDG)rFl=G4UniformRand()<bs*bs*bs/m3_value;
01993 #ifdef debug
01994 G4cout<<"G4Q::HQ: sPDG="<<sPDG<<", hsflag="<<hsflag<<", rPDG="<<rPDG<<curQ<<",rM="
01995 <<rMass<<",rb="<<rb<<",F="<<rFl<<",v="<<m2_value<<", bs="<<bs<<G4endl;
01996 #endif
01997 if(!hsflag&&rFl&&rPDG&& rPDG!=10 && rb<2 && aPDG!=1114 && aPDG!=2224 && aPDG!=3334)
01998 {
01999 G4int regPDG = 0;
02000 G4int refPDG = 0;
02001 G4int redPDG = 0;
02002 G4int repPDG = 0;
02003 if(rPDG && rPDG!=10)
02004 {
02005 if (rPDG== 3122) rPDG= 3212;
02006 else if(rPDG==-3122) rPDG=-3212;
02007 if(rPDG>0)repPDG=rPDG+2;
02008 else repPDG=rPDG-2;
02009 if(repPDG>0)redPDG=repPDG+2;
02010 else redPDG=repPDG-2;
02011 if(redPDG>0)refPDG=redPDG+2;
02012 else refPDG=redPDG-2;
02013 if(refPDG>0)regPDG=refPDG+2;
02014 else regPDG=refPDG-2;
02015
02016 #ifdef debug
02017 G4cout<<"G4Q::HQ:QuasM="<<quasM<<valQ<<")->H("<<sPDG<<")+R("<<rPDG<<")"<<",rp="
02018 <<repPDG<<",rd="<<redPDG<<",rf="<<refPDG<<",rg="<<regPDG<<G4endl;
02019 #endif
02020 G4double resM2 = G4QPDGCode( rPDG).GetMass2();
02021 G4double repM2 = G4QPDGCode(repPDG).GetMass2();
02022 G4double redM2 = G4QPDGCode(redPDG).GetMass2();
02023 G4double refM2 = G4QPDGCode(refPDG).GetMass2();
02024 sB = sqrt((resM2+repM2)/2.);
02025 G4double pB = sqrt((repM2+redM2)/2.);
02026 G4double dB = sqrt((redM2+refM2)/2.);
02027 G4double fB = sqrt(refM2)+150.;
02028 if(!cB) fB= sqrt((refM2+G4QPDGCode(regPDG).GetMass2())/2.);
02029 G4double dif=quasM-sMass;
02030 G4double rM = GetRandomMass(repPDG,dif);
02031 G4double dM = GetRandomMass(redPDG,dif);
02032 G4double fM = GetRandomMass(refPDG,dif);
02033 #ifdef debug
02034 G4cout<<"G4Q::HQ: rM="<<rM<<",rMa="<<rMass<<",sB="<<sB<<"(bQ="<<bQ<<"),pB="<<pB
02035 <<",dM="<<dM<<",dB="<<dB<<",fM="<<fM<<",fB="<<fB<<G4endl;
02036 #endif
02037 if(((rM>0 && rMass<pB && rMass>sB) || (dM>0 && rMass>pB && rMass<dB) ||
02038 (fM>0 && rMass>dB && rMass<fB)) && theEnvironment.GetPDG()==NUCPDG)
02039 {
02040 if (rMass>pB && rMass<dB && dM>0)
02041 {
02042 repPDG=redPDG;
02043 rM=dM;
02044 }
02045 else if(rMass>dB && rMass<fB && dM>0)
02046 {
02047 repPDG=refPDG;
02048 rM=fM;
02049 }
02050 #ifdef debug
02051 G4cout<<"G4Q::HQ:s="<<sPDG<<",Q=>rM="<<rMass<<"(minQ="<<rPDG<<curQ<<")+sB="
02052 <<sB<<G4endl;
02053 #endif
02054 if(quasM<rM+sMass &&(sPDG==221||sPDG==331))
02055 {
02056 sPDG=111;
02057 sMass=mPi0;
02058 }
02059 G4LorentzVector r4Mom(0.,0.,0.,rM);
02060 G4LorentzVector s4Mom(0.,0.,0.,sMass);
02061 G4double sum=rM+sMass;
02062 if(fabs(quasM-sum)<eps)
02063 {
02064 r4Mom=q4Mom*(rM/sum);
02065 s4Mom=q4Mom*(sMass/sum);
02066 }
02067 else if(quasM<sum || !G4QHadron(q4Mom).DecayIn2(r4Mom, s4Mom))
02068 {
02069
02070
02071
02072 G4ExceptionDescription ed;
02073 ed << "H+Res Decay failed: rPD=" << repPDG << "(rM=" << rMass
02074 << ")+sPD=" << sPDG << "(sM=" << sMass << "), Env="
02075 << theEnvironment << G4endl;
02076 G4Exception("G4Quasmon::HadronizeQuasmon()", "HAD_CHPS_0015",
02077 FatalException, ed);
02078 }
02079 #ifdef debug
02080 G4cout<<"G4Q::HQ:=== 1 ===> HadronVec, Q="<<q4Mom<<" -> s4M="<<s4Mom<<"("
02081 <<sPDG<<"), r4M="<<r4Mom<<"("<<repPDG<<")"<<G4endl;
02082 #endif
02083
02084 G4QHadron* curHadr1 = new G4QHadron(repPDG,r4Mom,curQ);
02085 FillHadronVector(curHadr1);
02086
02087 if (sPDG==2112) sPDG=90000001;
02088 else if(sPDG==2212) sPDG=90001000;
02089 else if(sPDG==3122) sPDG=91000000;
02090
02091 G4QHadron* curHadr2 = new G4QHadron(sPDG,s4Mom);
02092 FillHadronVector(curHadr2);
02093 ClearQuasmon();
02094 qEnv=theEnvironment;
02095 return theQHadrons;
02096 }
02097 }
02098 }
02099 curQ = memQ;
02100
02101 fdul = rFl && rPDG!=10;
02102 }
02103
02104
02105
02106 if (kn > minSqN && ku > minSqT)
02107
02108
02109
02110
02111 {
02112 pCond=false;
02113 #ifdef debug
02114
02115
02116
02117
02118 G4cout<<"G4Q::HQ:P-Attempt#"<<pCount<<" *Yes* sPDG="<<sPDG<<",ku="<<ku<<">"
02119 <<minSqT<<" || kn="<<kn<<">"<<minSqN<<G4endl;
02120
02121
02122
02123 #endif
02124 }
02125 else
02126 {
02127 #ifdef debug
02128
02129
02130
02131
02132 G4cout<<"G4Q::HQ:P-Attempt#"<<pCount<<",No. ku="<<ku<<"<"<<minSqT<<" or kn="<<kn
02133 <<"<"<<minSqN<<" or E="<<fr4Mom.e()<<"<"<<sCBE<<G4endl;
02134
02135
02136
02137 #endif
02138
02139
02140
02141
02142
02143
02144
02145 if (kn < minSqN && ku < minSqT)
02146 {
02147
02148
02149
02150
02151
02152
02153
02154
02155 if(ku < minSqT && ku > PMEMktM2)
02156 {
02157
02158 PMEMktM2=ku;
02159
02160
02161
02162 PMEMknM2=kn;
02163 PMEMfr4M=fr4Mom;
02164 PMEMrQ4M=rQ4Mom;
02165 PMEMreM2=reTNM2;
02166 PMEMrMas=rMass;
02167 PMEMpMas=pMass;
02168 PMEMsMas=sMass;
02169 PMEMdMas=dMass;
02170 PMEMmiSN=minSqN;
02171 PMEMmiST=minSqT;
02172 PMEMmiSB=minSqB;
02173 PMEMrPDG=rPDG;
02174 PMEMsPDG=sPDG;
02175 PMEMpPDG=pPDG;
02176 PMEMpQC =pQC;
02177 PMEMsQC =sQC;
02178 PMEMtQC =transQC;
02179 PMEMcQC =curQ;
02180 PMEMhsfl=hsflag;
02181 PMEMnucf=nucflag;
02182 #ifdef debug
02183 G4cout<<"G4Q::HQ:RemTheBest rPDG="<<rPDG<<",sPDG="<<sPDG<<",kt="<<kt<<G4endl;
02184 #endif
02185 }
02186 else if(!pCount)
02187 {
02188
02189 PMEMktM2=ku;
02190
02191
02192
02193 PMEMknM2=kn;
02194 PMEMfr4M=fr4Mom;
02195 PMEMrQ4M=rQ4Mom;
02196 PMEMreM2=reTNM2;
02197 PMEMrMas=rMass;
02198 PMEMpMas=pMass;
02199 PMEMsMas=sMass;
02200 PMEMdMas=dMass;
02201 PMEMmiSN=minSqN;
02202 PMEMmiST=minSqT;
02203 PMEMmiSB=minSqB;
02204 PMEMrPDG=rPDG;
02205 PMEMsPDG=sPDG;
02206 PMEMpPDG=pPDG;
02207 PMEMpQC =pQC;
02208 PMEMsQC =sQC;
02209 PMEMtQC =transQC;
02210 PMEMcQC =curQ;
02211 PMEMhsfl=hsflag;
02212 PMEMnucf=nucflag;
02213 #ifdef debug
02214 G4cout<<"G4Q::HQ:RemTheFirst rPDG="<<rPDG<<",sPDG="<<sPDG<<",kt="<<kt<<G4endl;
02215 #endif
02216 }
02217 else
02218 {
02219 fr4Mom=PMEMfr4M;
02220 rQ4Mom=PMEMrQ4M;
02221 reTNM2=PMEMreM2;
02222 rMass =PMEMrMas;
02223 pMass =PMEMpMas;
02224 sMass =PMEMsMas;
02225 dMass =PMEMdMas;
02226 minSqN=PMEMmiSN;
02227 minSqT=PMEMmiST;
02228 minSqB=PMEMmiSB;
02229 rPDG =PMEMrPDG;
02230 sPDG =PMEMsPDG;
02231
02232 ku=PMEMktM2;
02233
02234
02235
02236 kn =PMEMknM2;
02237 pPDG=PMEMpPDG;
02238 pQC=PMEMpQC;
02239 sQC=PMEMsQC;
02240 transQC=PMEMtQC;
02241 curQ=PMEMcQC;
02242 hsflag=PMEMhsfl;
02243 nucflag=PMEMnucf;
02244 }
02245 }
02246 }
02247 pCount++;
02248 }
02249 #ifdef debug
02250 G4cout<<"G4Q::HQ:>rPDG="<<rPDG<<curQ<<",sPDG="<<sPDG<<",kt="<<kt<<",F="<<fprob
02251 <<",totQC="<<totQC<<",sQC="<<sQC<<G4endl;
02252 #endif
02253 if(fprob)
02254 {
02255 rPDG=curQ.GetSPDGCode();
02256 G4double rrr=G4UniformRand();
02257 if(rPDG==111&&sPDG!=111&&rrr>.5) rPDG=221;
02258 if(rPDG==221&&sPDG!=221&&sPDG!=331&&rrr<.5) rPDG=111;
02259 }
02260
02261 G4double reMass=sqrt(minSqB);
02262 if (!rPDG)
02263 {
02264 G4ExceptionDescription ed;
02265 ed << "Unidentifiable residual Hadron: Q=" << curQ << ",r=" << rPDG
02266 << "+s=" << sPDG << "(sM=" << sMass << ")" << G4endl;
02267 G4Exception("G4Quasmon::HadronizeQuasmon()", "HAD_CHPS_0016", FatalException, ed);
02268 }
02269 if(rPDG==221||rPDG==331) reMass=mPi0;
02270 G4double aMass=0.;
02271
02272
02273
02274 if ( ( ( (sPDG < MINPDG && envPDG > MINPDG && envPDG != NUCPDG) ||
02275 (sPDG > MINPDG && sPDG!=NUCPDG && envPDG > pPDG)
02276 ) && iniBN > 0
02277 ) || iniBN > 1 || rPDG == 10
02278 ) aMass=0.;
02279
02280 #ifdef debug
02281 G4cout <<"G4Q::HQ:Is hsfl="<<hsflag<<" or fdul="<<fdul<<" or [rM="<<rMass<<"<"<<reMass
02282 <<" + "<<aMass<<" or rM2="<<reTNM2<<" < miM2="<<tmpTM2<<" and ePDG="<<envPDG
02283 <<">pPDG="<<pPDG<<"] to fail?"<<G4endl;
02284 #endif
02285
02286 if ( hsflag ||
02287 (sPDG < MINPDG && rMass < reMass+aMass) ||
02288 (sPDG > MINPDG && envPDG > pPDG && reTNM2 < tmpTM2) ||
02289 fdul )
02290 {
02291
02292
02293 #ifdef debug
02294 G4cout<<"G4Q::HQ: Yes(No), hsf="<<hsflag<<",sPDG="<<sPDG<<",pM="<<pMass<<",Env="
02295 <<envPDG<<",QM="<<quasM<<valQ<<", fpr="<<fprob<<G4endl;
02296 #endif
02297 G4QPDGCode rQPDG=G4QPDGCode(rPDG);
02298 if(hsflag) rMass=rQPDG.GetMass();
02299 if(sPDG>MINPDG&&sPDG!=NUCPDG)
02300 {
02301 G4QContent tmpEQ=envQC-pQC;
02302 G4QNucleus tmpN(tmpEQ);
02303 G4double tmpNM=tmpN.GetMZNS();
02304 G4QContent tmpRQ=valQ+transQC;
02305
02306
02307 G4LorentzVector ResEnv4Mom(0.,0.,0.,tmpNM);
02308 if(rQ4Mom==zeroLV)
02309 {
02310 #ifdef pdebug
02311 G4ExceptionDescription ed;
02312 ed << "Why Fail? (5): NEEDS-EVAP-5,Q=" << q4Mom << valQ << ",QEnv="
02313 << theEnvironment << G4endl;
02314 G4Exception("G4Quasmon::HadronizeQuasmon()","HAD_CHPS_0017", FatalException, ed);
02315 #endif
02316 qEnv=theEnvironment;
02317 return theQHadrons;
02318 }
02319 G4LorentzVector retN4Mom=rQ4Mom+ResEnv4Mom;
02320 G4double retNM2=retN4Mom.m2();
02321 G4double retNM=sqrt(retNM2);
02322 G4QContent tmpTQ=tmpEQ+tmpRQ;
02323 G4QNucleus tmpT(tmpTQ);
02324
02325 G4double tmpTM=tmpT.GetMZNS();
02326 if(tmpTM>retNM) tmpT=G4QNucleus(tmpTQ,retN4Mom);
02327 G4QPDGCode sQPDG(sPDG);
02328
02329
02331
02332
02334 rMass=rQPDG.GetMass();
02335
02336 G4int rB=rQPDG.GetBaryNum();
02337 G4double rCB=theEnvironment.CoulombBarrier(rQPDG.GetCharge(),rB);
02338 if(rCB < 0.) rCB=0.;
02339 G4int sB=sQPDG.GetBaryNum();
02340 G4double sCB=theEnvironment.CoulombBarrier(sQPDG.GetCharge(),sB);
02341 if(sCB < 0.) sCB=0.;
02342
02343 #ifdef debug
02344
02345
02346
02347
02348
02349
02350
02351
02352
02353
02354 G4cout<<"G4Q::HQ:tM="<<totMass<<",totQC="<<totQC<<",rtQC="<<tmpTQ<<",pQC="<<pQC
02355 <<",sB="<<sB<<",resB="<<tmpT.GetA()<<G4endl;
02356 #endif
02357
02358 if(nQuasms==1 && tmpNM+rMass+rCB+sMass+sCB < totMass &&
02359 (sB==1 || (sB==2 && G4UniformRand()<.2) || (sB==3 && G4UniformRand()<.1)) &&
02360 (rB==1 || (rB==2 && G4UniformRand()<.2) || (rB==3 && G4UniformRand()<.1)) )
02361
02362 {
02363 G4LorentzVector fr4M = G4LorentzVector(0.,0.,0.,sMass);
02364 G4LorentzVector re4M = G4LorentzVector(0.,0.,0.,tmpNM);
02365 G4LorentzVector rq4M = G4LorentzVector(0.,0.,0.,rMass);
02366 #ifdef debug
02367 G4double cfM=fr4Mom.m();
02368 G4double ctM=tot4M.m();
02369 G4cout<<"G4Q::HQ: *YES*,tM="<<ctM<<"="<<totMass<<",fM="<<cfM<<"="<<sMass<<G4endl;
02370 #endif
02371 G4double sum=tmpNM+sMass+rMass;
02372 if(fabs(totMass-sum)<eps)
02373 {
02374 re4M=tot4M*(tmpNM/sum);
02375 rq4M=tot4M*(rMass/sum);
02376 fr4M=tot4M*(sMass/sum);
02377 }
02378 else if(totMass<sum || !G4QHadron(tot4M).DecayIn3(rq4M,re4M,fr4M))
02379 {
02380 G4ExceptionDescription ed;
02381 ed << "DecayIn Frag+ResQ+ResE failed: Decay (" << totMass << ") in Fragm("
02382 << sMass <<")+ResQ("<< rMass <<")+ResEnv("<< tmpNM <<")="<< sum << G4endl;
02383 G4Exception("G4Quasmon::HadrQuasmon()", "HAD_CHPS_0018", FatalException, ed);
02384 }
02385 G4QHadron* resQH = new G4QHadron(tmpRQ,rq4M);
02386 FillHadronVector(resQH);
02387 if(nQuasms==1)
02388 {
02389 G4QHadron* envH = new G4QHadron(tmpEQ,re4M);
02390 FillHadronVector(envH);
02391 qEnv = vacuum;
02392 }
02393 else
02394 {
02395 qEnv=G4QNucleus(tmpEQ,re4M);
02396 #ifdef debug
02397 G4cout<<"**G4Q::HQ:(3)**KeepEnvironmentMoving**, nQ="<<nQuasms<<G4endl;
02398 #endif
02399 }
02400 G4QHadron* candH = new G4QHadron(sPDG,fr4M);
02401 FillHadronVector(candH);
02402 ClearQuasmon();
02403 return theQHadrons;
02404 }
02405
02406
02407 else if(nQuasms==1 && tmpTM+sMass+sCB < totMass)
02408
02409 {
02410 qEnv = G4QNucleus(tmpTQ,retN4Mom);
02411 #ifdef debug
02412 G4cout<<"G4Q::HQ:(2)*KeepEnvironmentMoving*,nQ="<<nQuasms<<",Env="<<qEnv<<G4endl;
02413 #endif
02414
02415
02416
02417 #ifdef debug
02418 G4LorentzVector d4M=tot4M-retN4Mom-fr4Mom;
02419 G4QContent dQC=totQC-tmpTQ-sQC;
02420 G4cout<<"G4Q::HQ: rTotM="<<retN4Mom.m()<<" >? GSM="<<tmpTM<<",d4M="<<d4M<<",dQC="
02421 <<dQC<<G4endl;
02422 #endif
02423 G4QHadron* candHadr = new G4QHadron(sPDG,fr4Mom);
02424 FillHadronVector(candHadr);
02425 #ifdef debug
02426 G4double frM=fr4Mom.m();
02427 G4LorentzVector dif2=tot4M-retN4Mom-fr4Mom;
02428 G4cout<<"G4Q::HQ:sM="<<sMass<<"="<<frM<<", fT="<<fr4Mom.e()-frM<<",dif24M="<<dif2
02429 <<G4endl;
02430 #endif
02431 ClearQuasmon();
02432 return theQHadrons;
02433 }
02434 else if(totBN>1 &&totMass>totM &&totS>=0&&envPDG>MINPDG&&envPDG!=NUCPDG)
02435
02436 {
02437 #ifdef ppdebug
02438
02439 G4double fraM=fr4Mom.m();
02440 G4double kinE=fr4Mom.e()-fraM;
02441 G4double sumM=tmpTM+fraM;
02442 G4ExceptionDescription ed;
02443 ed << "Why Fail?(6): ProductMasses>totalMass: EV-6: TotEVAPORATION:s=" << sPDG
02444 << ",T=" << kinE << ",RM=" << retN4Mom.m() << "<" << tmpTM << ",tQC="
02445 << transQC << ",E=" << excE << ",sM=" << sumM << ">tM=" << totMass << ",nQ="
02446 << nQuasms << G4endl;
02447 G4Exception("G4Quasmon::HadronizeQuasmon()","HAD_CHPS_0019", FatalException, ed);
02448 #endif
02449 #ifdef debug
02450 G4cout<<"G4Q::HQ:Q="<<q4Mom<<quasM<<",E="<<theEnvironment<<",P="<<phot4M<<G4endl;
02451 #endif
02452 qEnv=theEnvironment;
02453 return theQHadrons;
02454 }
02455 else if(totBN==1 && nQuasms==1)
02456 {
02457 #ifdef debug
02458 G4cout<<"G4Q::HQ:tB=1,nQ=1,Z="<<totZ<<",S="<<totS<<totQC<<",M="<<totMass<<G4endl;
02459 #endif
02460 G4double nucM= mProt;
02461 G4double piM = 0.;
02462 G4int nucPDG = 2212;
02463 G4int piPDG = 22;
02464 if(abs(totS)==1)
02465 {
02466 if(totS==1)
02467 {
02468 if(!totZ&&totMass>mLamb+mPi0)
02469 {
02470 nucM = mLamb;
02471 nucPDG= 3122;
02472 piM = mPi0;
02473 piPDG = 111;
02474 }
02475 else if(abs(totZ)==1&&totMass>mLamb+mPi)
02476 {
02477 nucM = mLamb;
02478 nucPDG= 3122;
02479 piM = mPi;
02480 if(totZ>0) piPDG = 211;
02481 else piPDG =-211;
02482 }
02483 else
02484 {
02485 G4ExceptionDescription ed;
02486 ed << "Pi + Lambda decay error:Z=" << totZ << ",S=" << totS
02487 << totQC << ",tM=" << totMass << G4endl;
02488 G4Exception("G4Quasmon::HadronizeQuasmon()", "HAD_CHPS_0020",
02489 FatalException, ed);
02490 }
02491 }
02492 else
02493 {
02494 if(!totZ&&totMass>mNeut+mK0)
02495 {
02496 nucM = mNeut;
02497 nucPDG= 2112;
02498 piM = mK0;
02499 piPDG = 311;
02500 }
02501 else if(totZ==2&&totMass>mProt+mK)
02502 {
02503 piM = mK;
02504 piPDG = 321;
02505 }
02506 else if(totZ==1&&totMass>mProt+mK0&&G4UniformRand()>0.5)
02507 {
02508 piM = mK0;
02509 piPDG = 311;
02510 }
02511 else if(totZ==1&&totMass>=mNeut+mK)
02512 {
02513 nucM = mNeut;
02514 nucPDG= 2112;
02515 piM = mK;
02516 piPDG = 321;
02517 }
02518 else
02519 {
02520 G4ExceptionDescription ed;
02521 ed << "K + Nucleon decay error: Z=" << totZ << ",S=" << totS
02522 << totQC << ",tM=" << totMass << G4endl;
02523 G4Exception("G4Quasmon::HadronizeQuasmon", "HAD_CHPS_0021",
02524 FatalException, ed);
02525 }
02526 }
02527 }
02528 else if(totMass>PiNM&&!totS)
02529 {
02530 if(!totZ&&totMass>mProt+mPi&&G4UniformRand()<0.5)
02531 {
02532 piM = mPi;
02533 piPDG = -211;
02534 }
02535 else if(!totZ&&totMass>mNeut+mPi0)
02536 {
02537 nucM = mNeut;
02538 nucPDG= 2112;
02539 piM = mPi0;
02540 piPDG = 111;
02541 }
02542 else if(totZ==1&&totMass>mNeut+mPi&&G4UniformRand()<0.5)
02543 {
02544 nucM = mNeut;
02545 nucPDG= 2112;
02546 piM = mPi;
02547 piPDG = 211;
02548 }
02549 else if(totZ==1&&totMass>mProt+mPi0)
02550 {
02551 piM = mPi0;
02552 piPDG = 111;
02553 }
02554 else if(totZ==-1)
02555 {
02556 nucM = mNeut;
02557 nucPDG= 2112;
02558 piM = mPi;
02559 piPDG = -211;
02560 }
02561 else if(totZ==2)
02562 {
02563 piM = mPi;
02564 piPDG = 211;
02565 }
02566 else
02567 {
02568 G4ExceptionDescription ed;
02569 ed << "Pi + Nucleon decay error: Z=" << totZ << ",B=" << totBN
02570 << ",E=" << envQC << ",Q=" << valQ << G4endl;
02571 G4Exception("G4Quasmon::HadronizeQuasmon()", "HAD_CHPS_0022",
02572 FatalException, ed);
02573 }
02574 }
02575 else if(!totS)
02576 {
02577 if(!totZ)
02578 {
02579 nucM=mNeut;
02580 nucPDG=2112;
02581 }
02582 else if(totZ<0||totZ>1)
02583 {
02584 G4ExceptionDescription ed;
02585 ed << "Photon+Nucleon decay error: Z=" << totZ << ",B=" << totBN
02586 << ",E=" << envQC <<",Q=" << valQ << G4endl;
02587 G4Exception("G4Quasmon::HadronizeQuasmon()", "HAD_CHPS_0023",
02588 FatalException, ed);
02589 }
02590 }
02591 G4LorentzVector pi4M(0.,0.,0.,piM);
02592 G4LorentzVector nuc4M(0.,0.,0.,nucM);
02593 G4double sum=piM+nucM;
02594 if(fabs(totMass-sum)<eps)
02595 {
02596 pi4M=tot4M*(piM/sum);
02597 nuc4M=tot4M*(nucM/sum);
02598 }
02599 else if(totMass<sum || !G4QHadron(tot4M).DecayIn2(pi4M, nuc4M))
02600 {
02601 G4ExceptionDescription ed;
02602 ed << "Gam/Pi/K+N decay error: T=" << tot4M << totMass
02603 << "->gam/pi/K(" << piM << ")+N=" << nucPDG << "(" << nucM
02604 << ")=" << sum << G4endl;
02605 G4Exception("G4Quasmon::HadronizeQuasmon()", "HAD_CHPS_0024",
02606 FatalException, ed);
02607 }
02608 #ifdef debug
02609 G4cout<<"G4Q::HQ:T="<<tot4M<<totMass<<"->GPK="<<piPDG<<pi4M<<"+B="<<nucPDG<<nuc4M
02610 <<G4endl;
02611 #endif
02612 G4QHadron* piH = new G4QHadron(piPDG,pi4M);
02613 FillHadronVector(piH);
02614 G4QHadron* nucH = new G4QHadron(nucPDG,nuc4M);
02615 FillHadronVector(nucH);
02616 ClearQuasmon();
02617 qEnv=vacuum;
02618 return theQHadrons;
02619 }
02620 #ifdef debug
02621 else G4cout<<"***G4Q::HQ: B="<<totBN<<",tM="<<totMass<<" > M="<<totM<<",S="<<totS
02622 <<", envPDG="<<envPDG<<G4endl;
02623 #endif
02624 }
02625 G4double dm=quasM-sMass;
02626 #ifdef debug
02627 G4cout<<"G4Q::HQ:f="<<fprob<<",d="<<dm<<",rPDG="<<rPDG<<",rM="<<rMass<<",M="<<reMass
02628 <<",sM="<<sMass<<G4endl;
02629 #endif
02630 if(abs(dm)<.000001)
02631 {
02632 if(sPDG==iniPDG)
02633 {
02634 G4QHadron* quasH = new G4QHadron(iniPDG,q4Mom);
02635 FillHadronVector(quasH);
02636 ClearQuasmon();
02637 qEnv=theEnvironment;
02638 return theQHadrons;
02639 }
02640 else G4cerr<<"---Warning---G4Q::HQ:Q=H,q="<<iniPDG<<",s="<<sPDG<<",d="<<dm<<G4endl;
02641 }
02642 G4double rWi=0.;
02643 if(rPDG!=10) rWi=G4QPDGCode(rPDG).GetWidth();
02644 if(rPDG!=10&&rMass>dm&&!rWi)
02645 {
02646 G4double sWi=G4QPDGCode(sPDG).GetWidth();
02647 G4double sMM=G4QPDGCode(sPDG).GetMass();
02648 if(sWi)
02649 {
02650 G4double mmm=theWorld->GetQParticle(G4QPDGCode(sPDG))->MinMassOfFragm();
02651 G4double ddm=quasM-rMass;
02652 if(fabs(sMM-ddm)<1.5*sWi-.001 && ddm>mmm)
02653 {
02654 #ifdef debug
02655 G4double msm=sMass;
02656 #endif
02657 sMass=GetRandomMass(sPDG,ddm);
02658 if(fabs(sMass)<.001)
02659 {
02660 #ifdef debug
02661 G4cerr<<"***G4Q::HQ:ChangeToM=0, "<<sPDG<<",new="<<ddm<<",old="<<msm<<G4endl;
02662 #endif
02663 sMass=ddm;
02664 }
02665 if(sMass<ddm) sMass=ddm;
02666 #ifdef debug
02667 G4cout<<"G4Q::HQ: sPDG="<<sPDG<<",sM="<<sMass<<",d="<<ddm<<",isM="<<msm<<",W="
02668 <<sWi<<G4endl;
02669 #endif
02670 }
02671
02672
02673
02674
02675
02676
02677
02678
02679
02680
02681 }
02682 }
02683
02684 G4double rnd=G4UniformRand();
02685
02686 #ifdef debug
02687 G4cout<<"G4Q::HQ:BEFrPDGcor,d="<<dm<<",R="<<rnd<<",r="<<rPDG<<",rM="<<rMass<<G4endl;
02688 #endif
02689
02690
02691 if(rPDG==111 && sPDG!=111 && dm>548.)
02692 {
02693
02694 if(dm>958.) rPDG=331;
02695 else rPDG=221;
02696 }
02697 if(rPDG==221 && dm>958. && rnd>.5 ) rPDG=331;
02698 if(rPDG==331 &&(dm<958. || rnd<.5)) rPDG=221;
02699
02700 if(rPDG==221 && dm<548.) rPDG=111;
02701
02702
02703 if ( ( (rPDG == 111 && sPDG!= 111) || rPDG == 221) &&
02704 rMass > 544. && dm > 544. && rnd > .5) rPDG=113;
02705
02706 if ( ( (rPDG == 111 && sPDG != 111) || rPDG == 221) &&
02707 rMass > 782. && dm > 782. && rnd < .5) rPDG = 223;
02708
02709 if ( rPDG == 331 && rMass > 1020. && dm > 1020. && rnd < .5) rPDG=333;
02710
02711 if(rPDG== 211 && dm>544. && rnd>.5) rPDG= 213;
02712 if(rPDG==-211 && dm>544. && rnd>.5) rPDG=-213;
02713 #ifdef debug
02714 G4cout<<"G4Q::HQ:rCor,Q="<<quasM<<",sM="<<sMass<<",r="<<rPDG<<",rM="<<rMass<<G4endl;
02715 #endif
02716 if (rPDG < MINPDG && rPDG != 2212 && rPDG != 2112 && rPDG != 3122 && rPDG != 10)
02717 {
02718 reMass=GetRandomMass(rPDG,dm);
02719 #ifdef debug
02720 G4cout<<"G4Q::HQ:dm="<<dm<<", ResQM="<<reMass<<" is changed to PDG="<<rPDG<<G4endl;
02721 #endif
02722 if(reMass==0.)
02723 {
02724 if(sPDG==221 || sPDG==331)
02725 {
02726 if (sPDG==221) dm+=mEta-mPi0;
02727 else if(sPDG==331) dm+=mEtaP-mPi0;
02728 if(dm<0)
02729 {
02730 dm+=mPi0;
02731 sPDG=22;
02732 sMass=0.;
02733 }
02734 else
02735 {
02736 sPDG=111;
02737 sMass=mPi0;
02738 }
02739 if(dm<mPi0-.00001&&rPDG==111)
02740 {
02741 rPDG=22;
02742 reMass=0.;
02743 }
02744 else reMass=GetRandomMass(rPDG,dm);
02745 if(reMass==0.)G4cerr<<"-W-G4Q::HQ:2,M="<<quasM<<",r="<<rPDG<<",d="<<dm<<G4endl;
02746 }
02747 else if(rPDG==111)
02748 {
02749 rPDG=22;
02750 reMass=0.;
02751 }
02752 else
02753 {
02754 if(CheckGroundState()) ClearQuasmon();
02755
02756 #ifdef pdebug
02757 G4ExceptionDescription ed;
02758 ed << "Why Fail? (7): NeedsEvap7:s=" << sPDG << ",Q=" << q4Mom
02759 << valQ << ",r=" << rPDG << G4endl;
02760 G4Exception("G4Quasmon::HadronizeQuasmon()","HAD_CHPS_0025",FatalException,ed);
02761 #endif
02762 qEnv=theEnvironment;
02763 return theQHadrons;
02764 }
02765 }
02766 }
02767 else if(rPDG==NUCPDG)
02768 {
02769 if(dm>mPi0)
02770 {
02771 rPDG=111;
02772 reMass=mPi0;
02773 }
02774 else
02775 {
02776 rPDG=22;
02777 reMass=0.;
02778 }
02779 }
02780 G4double freeRQM=rQPDG.GetMass();
02781 G4int RQB = rQPDG.GetBaryNum();
02782 G4double fRQW= 3*rQPDG.GetWidth();
02783 if(fRQW<.001) fRQW=.001;
02784 G4QPDGCode sQPDG(sPDG);
02785 G4int sChg=sQPDG.GetCharge();
02786 G4int sBaryn=sQPDG.GetBaryNum();
02787 G4double sCB=theEnvironment.CoulombBarrier(sChg,sBaryn);
02788 #ifdef debug
02789 G4cout<<"G4Q::HQ:h="<<sCB<<",C="<<sChg<<",B="<<sBaryn<<",E="<<theEnvironment<<G4endl;
02790 #endif
02791 G4int rChg=rQPDG.GetCharge();
02792 G4int rBaryn=rQPDG.GetBaryNum();
02793 G4double rCB=theEnvironment.CoulombBarrier(rChg,rBaryn);
02794 #ifdef debug
02795 G4cout<<"G4Q::HQ:rqCB="<<rCB<<",rqC="<<rChg<<",rqB="<<sBaryn<<",rM="<<rQPDG<<",reM="
02796 <<reMass<<G4endl;
02797 #endif
02798 if ( totBN > 1 && totS >= 0 && envPDG > MINPDG && envPDG != NUCPDG &&
02799 (reMass+sMass > quasM || sCB+rCB+reMass+sMass+envM > totMass ||
02800 (!RQB && quasM < diPiM)
02801 )
02802 )
02803
02804 {
02805 #ifdef pdebug
02806 G4ExceptionDescription ed;
02807 ed << "Why Fail? (8): RQM+SM=" << reMass+sMass << ">QM=" << quasM << ", sCB="
02808 << sCB << " + rCB=" << rCB << " + rM=" << reMass << " + sMass=" << sMass
02809 << " + eM=" << envM << " = " << sCB+rCB+reMass+sMass+envM << ">tM=" << totMass
02810 << "," << reMass+sMass+envM << G4endl;
02811 G4Exception("G4Quasmon::HadronizeQuasmon()", "HAD_CHPS_0026", FatalException, ed);
02812 #endif
02813 qEnv=theEnvironment;
02814 return theQHadrons;
02815 }
02816 if(rPDG==NUCPDG)
02817 {
02818 G4ExceptionDescription ed;
02819 ed << "Residual Particle is Vacuum: rPDG=90000000, MV=" << reMass << G4endl;
02820 G4Exception("G4Quasmon::HadronizeQuasmon()", "HAD_CHPS_0027", FatalException, ed);
02821 }
02822 if(rPDG==2212&&sPDG==311&&reMass+sMass>quasM)
02823 {
02824 if(mNeut+mK<=quasM+.001)
02825 {
02826 reMass=mNeut;
02827 rPDG =2112;
02828 rQPDG=G4QPDGCode(rPDG);
02829 rChg=rQPDG.GetCharge();
02830 rBaryn=rQPDG.GetBaryNum();
02831 rCB=theEnvironment.CoulombBarrier(rChg,rBaryn);
02832 #ifdef debug
02833 G4cout<<"G4Q::HQ:NCB="<<rCB<<",NC="<<rChg<<",sB="<<sBaryn<<",r="<<rQPDG<<G4endl;
02834 #endif
02835 freeRQM=mNeut;
02836 RQB=1;
02837 fRQW=0.;
02838 sMass =mK;
02839 if(mNeut+mK<=quasM) sMass=quasM-mNeut;
02840 sPDG =321;
02841 sQPDG=G4QPDGCode(sPDG);
02842 sChg=sQPDG.GetCharge();
02843 sBaryn=sQPDG.GetBaryNum();
02844 sCB=theEnvironment.CoulombBarrier(sChg,sBaryn);
02845 #ifdef debug
02846 G4cout<<"G4Q::HQ:KCB="<<sCB<<",KC="<<sChg<<",frB="<<sBaryn<<",E="<<theEnvironment
02847 <<G4endl;
02848 #endif
02849 curQ=neutQC;
02850 }
02851 else
02852 {
02853 G4ExceptionDescription ed;
02854 ed<<"Can't decay Q in N and K: (NK) QM="<< quasM<<",d="<< quasM-mNeut-mK<<G4endl;
02855 G4Exception("G4Quasmon::HadronizeQuasmon()","HAD_CHPS_0028", FatalException, ed);
02856 }
02857 }
02858 #ifdef debug
02859 G4cout<<"G4Q::HQ: ****** Before reM="<<reMass<<", rM="<<rMass<<G4endl;
02860 #endif
02861 G4QPDGCode tmpQPDG(rPDG);
02862 if(tmpQPDG.GetWidth()<.000001) reMass=tmpQPDG.GetMass();
02863 if(!reMass) reMass=rMass;
02864 #ifdef debug
02865 G4cout<<"G4Q::HQ: Decay in sM="<<sMass<<" + reM="<<reMass<<" (rM="<<rMass<<G4endl;
02866 #endif
02867 G4LorentzVector r4Mom(0.,0.,0.,reMass);
02868 G4LorentzVector s4Mom(0.,0.,0.,sMass);
02869 if(sPDG>MINPDG)
02870 {
02871 #ifdef debug
02872 G4cout<<"G4Q::HQ:Q->RQ+QEX s="<<sPDG<<",pM="<<pMass<<",E="<<theEnvironment<<G4endl;
02873 #endif
02874 q4Mom+=G4LorentzVector(0.,0.,0.,pMass);
02875 }
02876 G4double tmM=q4Mom.m()+.001;;
02877 G4double sum=reMass+sMass;
02878 if(fabs(tmM-sum)<eps)
02879 {
02880 r4Mom=q4Mom*(reMass/sum);
02881 s4Mom=q4Mom*(sMass/sum);
02882 }
02883 else if(tmM<sum || !G4QHadron(q4Mom).DecayIn2(r4Mom, s4Mom))
02884 {
02885 G4QContent resNQC=totQC-sQC;
02886 #ifdef debug
02887 G4cerr<<"---Warning---G4Q::HQ:M="<<tmM<<"=>rPDG="<<rPDG<<"(rM="<<reMass<<")+sPDG="
02888 <<sPDG<<"(sM="<<sMass<<")="<<sum<<",resNQC="<<resNQC<<G4endl;
02889 #endif
02890 G4QNucleus resTN(resNQC);
02891 G4double resTNM=resTN.GetMZNS();
02892 if(sPDG==311 && tmpQPDG.GetCharge()>0)
02893 {
02894 G4QContent crQC=tmpQPDG.GetQuarkContent()-KpQC+K0QC;
02895 G4QNucleus nNuc(crQC);
02896 G4double nreM=nNuc.GetGSMass();
02897 if(tmM>mK+nreM)
02898 {
02899 sMass=mK;
02900 sPDG=321;
02901 s4Mom=G4LorentzVector(0.,0.,0.,sMass);
02902 curQ+=K0QC-KpQC;
02903 reMass=nreM;
02904 rPDG=nNuc.GetPDG();
02905 r4Mom=G4LorentzVector(0.,0.,0.,reMass);
02906 sum=reMass+sMass;
02907 if(fabs(tmM-sum)<eps)
02908 {
02909 r4Mom=q4Mom*(reMass/sum);
02910 s4Mom=q4Mom*(sMass/sum);
02911 }
02912 else if(tmM<sum || !G4QHadron(q4Mom).DecayIn2(r4Mom, s4Mom))
02913 {
02914 G4ExceptionDescription ed;
02915 ed << "Hadron+K+ DecayIn2: (I) KCor M=" << tmM << "=>rPDG=" << rPDG << "(rM="
02916 << reMass << ")+sPDG=" << sPDG << "(sM=" << sMass << ")=" << sum <<G4endl;
02917 G4Exception("G4Quasmon::HadronizeQuasmon()", "HAD_CHPS_0029",
02918 FatalException, ed);
02919 }
02920 }
02921 else
02922 {
02923 G4ExceptionDescription ed;
02924 ed << "Hadron+K+ DecayIn2:(O) KCor M=" << tmM << "=>rPDG=" << rPDG << "(rM="
02925 << reMass << ")+sPDG=" << sPDG << "(sM=" << sMass << ")=" << sum << G4endl;
02926 G4Exception("G4Quasmon::HadronizeQuasmon()","HAD_CHPS_0030",FatalException,ed);
02927 }
02928 }
02929 else if(sPDG==321 && tmpQPDG.GetCharge()<=tmpQPDG.GetBaryNum())
02930 {
02931 G4QContent crQC=tmpQPDG.GetQuarkContent()-K0QC+KpQC;
02932 G4QNucleus nNuc(crQC);
02933 G4double nreM=nNuc.GetGSMass();
02934 if(tmM>mK0+nreM)
02935 {
02936 sMass=mK0;
02937 sPDG=311;
02938 s4Mom=G4LorentzVector(0.,0.,0.,sMass);
02939 curQ+=KpQC-K0QC;
02940 reMass=nreM;
02941 rPDG=nNuc.GetPDG();
02942 r4Mom=G4LorentzVector(0.,0.,0.,reMass);
02943 sum=reMass+sMass;
02944 if(fabs(tmM-sum)<eps)
02945 {
02946 r4Mom=q4Mom*(reMass/sum);
02947 s4Mom=q4Mom*(sMass/sum);
02948 }
02949 else if(tmM<sum || !G4QHadron(q4Mom).DecayIn2(r4Mom, s4Mom))
02950 {
02951 G4ExceptionDescription ed;
02952 ed << "Hadron+K0 DecayIn2: (I) K0Cor M=" << tmM << "=>rPDG=" << rPDG
02953 << "(rM=" << reMass << ")+sPDG=" << sPDG << "(sM=" << sMass << ")="
02954 << sum << G4endl;
02955 G4Exception("G4Quasmon::HadronizeQuasmon()", "HAD_CHPS_0031",
02956 FatalException, ed);
02957 }
02958 }
02959 else
02960 {
02961 G4ExceptionDescription ed;
02962 ed << "Hadron+K0 DecayIn2: (O) K0Cor M=" << tmM << "=>rPDG=" << rPDG << "(rM="
02963 << reMass << ")+sPDG=" << sPDG << "(sM=" << sMass << ")=" << sum << G4endl;
02964 G4Exception("G4Quasmon::HadronizeQuasmon()","HAD_CHPS_0032",FatalException,ed);
02965 }
02966 }
02967 else if(sPDG==211 && tmpQPDG.GetCharge()<tmpQPDG.GetBaryNum())
02968 {
02969 G4QContent crQC=tmpQPDG.GetQuarkContent()-Pi0QC+PiQC;
02970 G4QNucleus nNuc(crQC);
02971 G4double nreM=nNuc.GetGSMass();
02972 if(tmM>mPi0+nreM)
02973 {
02974 sMass=mPi0;
02975 sPDG=111;
02976 curQ+=PiQC;
02977 s4Mom=G4LorentzVector(0.,0.,0.,sMass);
02978 reMass=nreM;
02979 rPDG=nNuc.GetPDG();
02980 r4Mom=G4LorentzVector(0.,0.,0.,reMass);
02981 sum=reMass+sMass;
02982 if(fabs(tmM-sum)<eps)
02983 {
02984 r4Mom=q4Mom*(reMass/sum);
02985 s4Mom=q4Mom*(sMass/sum);
02986 }
02987 else if(tmM<sum || !G4QHadron(q4Mom).DecayIn2(r4Mom, s4Mom))
02988 {
02989 G4ExceptionDescription ed;
02990 ed << "Hadron+Pi0 DecayIn2: (I) Pi+/Pi0Cor M=" << tmM << "=>rPDG=" << rPDG
02991 << "(rM=" << reMass << ")+sPDG=" << sPDG << "(sM=" << sMass << ")=" << sum
02992 << G4endl;
02993 G4Exception("G4Quasmon::HadronizeQuasmon()", "HAD_CHPS_0033",
02994 FatalException, ed);
02995 }
02996 }
02997 else if(tmM>nreM)
02998 {
02999 sMass=0.;
03000 sPDG=22;
03001 curQ+=PiQC;
03002 s4Mom=G4LorentzVector(0.,0.,0.,sMass);
03003 reMass=nreM;
03004 rPDG=nNuc.GetPDG();
03005 r4Mom=G4LorentzVector(0.,0.,0.,reMass);
03006 sum=reMass+sMass;
03007 if(fabs(tmM-sum)<eps)
03008 {
03009 r4Mom=q4Mom*(reMass/sum);
03010 s4Mom=q4Mom*(sMass/sum);
03011 }
03012 else if(tmM<sum || !G4QHadron(q4Mom).DecayIn2(r4Mom, s4Mom))
03013 {
03014 G4ExceptionDescription ed;
03015 ed << "Hadron+Gamma DecayIn2: (I) Pi+/GamCor M=" << tmM << "=>rPDG=" << rPDG
03016 << "(rM=" << reMass << ")+sPDG=" << sPDG << "(sM=" << sMass << ")=" << sum
03017 << G4endl;
03018 G4Exception("G4Quasmon::HadronizeQuasmon()", "HAD_CHPS_0034",
03019 FatalException, ed);
03020 }
03021 }
03022 else
03023 {
03024 G4ExceptionDescription ed;
03025 ed << "Hadron+Pi0/Gam DecayIn2: (O) Pi+/Pi0Cor M=" << tmM << "=>rPDG=" << rPDG
03026 <<"(rM="<< reMass <<")+sPDG="<< sPDG <<"(sM="<< sMass <<")="<< sum <<G4endl;
03027 G4Exception("G4Quasmon::HadronizeQuasmon()","HAD_CHPS_0035",FatalException,ed);
03028 }
03029 }
03030 else if(sPDG==-211 && tmpQPDG.GetCharge()>0)
03031 {
03032 G4QContent crQC=tmpQPDG.GetQuarkContent()-Pi0QC+PiMQC;
03033 G4QNucleus nNuc(crQC);
03034 G4double nreM=nNuc.GetGSMass();
03035 if(tmM>mPi0+nreM)
03036 {
03037 sMass=mPi0;
03038 sPDG=111;
03039 curQ+=PiMQC;
03040 s4Mom=G4LorentzVector(0.,0.,0.,sMass);
03041 reMass=nreM;
03042 rPDG=nNuc.GetPDG();
03043 r4Mom=G4LorentzVector(0.,0.,0.,reMass);
03044 sum=reMass+sMass;
03045 if(fabs(tmM-sum)<eps)
03046 {
03047 r4Mom=q4Mom*(reMass/sum);
03048 s4Mom=q4Mom*(sMass/sum);
03049 }
03050 else if(tmM<sum || !G4QHadron(q4Mom).DecayIn2(r4Mom, s4Mom))
03051 {
03052 G4ExceptionDescription ed;
03053 ed << "Hadron+Pi0 DecayIn2: (I) Pi-/Pi0Cor M=" << tmM << "=>rPDG=" << rPDG
03054 << "(rM=" << reMass << ")+sPDG=" << sPDG << "(sM=" << sMass << ")=" << sum
03055 << G4endl;
03056 G4Exception("G4Quasmon::HadronizeQuasmon()", "HAD_CHPS_0036",
03057 FatalException, ed);
03058 }
03059 }
03060 else if(tmM>nreM)
03061 {
03062 sMass=0.;
03063 sPDG=22;
03064 curQ+=PiMQC;
03065 s4Mom=G4LorentzVector(0.,0.,0.,sMass);
03066 reMass=nreM;
03067 rPDG=nNuc.GetPDG();
03068 r4Mom=G4LorentzVector(0.,0.,0.,reMass);
03069 sum=reMass+sMass;
03070 if(fabs(tmM-sum)<eps)
03071 {
03072 r4Mom=q4Mom*(reMass/sum);
03073 s4Mom=q4Mom*(sMass/sum);
03074 }
03075 else if(tmM<sum || !G4QHadron(q4Mom).DecayIn2(r4Mom, s4Mom))
03076 {
03077 G4ExceptionDescription ed;
03078 ed << "Hadron+Gamma DecayIn2: (I) Pi-/GamCor M=" << tmM << "=>rPDG=" << rPDG
03079 << "(rM=" << reMass << ")+sPDG=" << sPDG << "(sM=" << sMass << ")=" << sum
03080 << G4endl;
03081 G4Exception("G4Quasmon::HadronizeQuasmon()", "HAD_CHPS_0037",
03082 FatalException, ed);
03083 }
03084 }
03085 }
03086 else if((sPDG==221 || sPDG==331) && tmM>mPi0+reMass)
03087 {
03088 sMass=mPi0;
03089 sPDG=111;
03090 s4Mom=G4LorentzVector(0.,0.,0.,sMass);
03091 sum=reMass+sMass;
03092 if(fabs(tmM-sum)<eps)
03093 {
03094 r4Mom=q4Mom*(reMass/sum);
03095 s4Mom=q4Mom*(sMass/sum);
03096 }
03097 else if(tmM<sum || !G4QHadron(q4Mom).DecayIn2(r4Mom, s4Mom))
03098 {
03099 G4ExceptionDescription ed;
03100 ed <<"Hadron+Pi0 DecayIn2: Eta/Pi0Cor M="<< tmM <<"=>rPDG="<< rPDG << "(rM="
03101 << reMass << ")+sPDG=" << sPDG << "(sM=" << sMass << ")=" << sum << G4endl;
03102 G4Exception("G4Quasmon::HadronizeQuasmon()","HAD_CHPS_0038",FatalException,ed);
03103 }
03104 }
03105 else if((sPDG==111 || sPDG==221 || sPDG==331) && tmM>reMass)
03106 {
03107 sMass=0.;
03108 sPDG=22;
03109 s4Mom=G4LorentzVector(0.,0.,0.,sMass);
03110 sum=reMass+sMass;
03111 if(fabs(tmM-reMass)<eps)
03112 {
03113 r4Mom=q4Mom*(reMass/sum);
03114 s4Mom=q4Mom*(sMass/sum);
03115 }
03116 else if(!G4QHadron(q4Mom).DecayIn2(r4Mom, s4Mom))
03117 {
03118 G4ExceptionDescription ed;
03119 ed << "QHadron+Gamma DecayIn2: PiCor M=" << tmM << "=>rPDG=" << rPDG << "(rM="
03120 << reMass <<")+sPDG="<< sPDG << "(sM=" << sMass << ")=" << reMass << G4endl;
03121 G4Exception("G4Quasmon::HadronizeQuasmon()","HAD_CHPS_0039",FatalException,ed);
03122 }
03123 }
03124 else if(iniBN>0 && iniS>0)
03125 {
03126 G4QContent tmpSQC=G4QPDGCode(sPDG).GetQuarkContent();
03127 G4QContent lanQC=tmpQPDG.GetQuarkContent()+tmpSQC+K0QC;
03128 G4QNucleus nucM(lanQC-PiMQC);
03129 G4double nreM=nucM.GetGSMass();
03130 G4QNucleus nucZ(lanQC-Pi0QC);
03131 G4double nreZ=nucZ.GetGSMass();
03132 #ifdef debug
03133 G4cout<<"G4Q::HQ:LsPDG="<<sPDG<<",rPDG="<<rPDG<<",Z="<<nucZ<<",M="<<nucM<<G4endl;
03134 #endif
03135 if((G4UniformRand()<.33333 || mPi+nreM>tmM) && mPi0+nreZ<tmM)
03136 {
03137 sMass=mPi0;
03138 sPDG=111;
03139 curQ+=tmpSQC+K0QC;
03140 s4Mom=G4LorentzVector(0.,0.,0.,sMass);
03141 reMass=nreZ;
03142 rPDG=nucZ.GetPDG();
03143 r4Mom=G4LorentzVector(0.,0.,0.,reMass);
03144 sum=reMass+sMass;
03145 if(fabs(tmM-sum)<eps)
03146 {
03147 r4Mom=q4Mom*(reMass/sum);
03148 s4Mom=q4Mom*(sMass/sum);
03149 }
03150 else if(tmM<sum || !G4QHadron(q4Mom).DecayIn2(r4Mom, s4Mom))
03151 {
03152 G4ExceptionDescription ed;
03153 ed <<"(L->n)+Pi0 DecayIn2: LamPi0 Cor M="<< tmM <<"=>rPDG="<< rPDG << "(rM="
03154 << reMass << ")+sPDG=" << sPDG << "(sM=" << sMass << ")=" << sum <<G4endl;
03155 G4Exception("G4Quasmon::HadronizeQuasmon()", "HAD_CHPS_0040",
03156 FatalException, ed);
03157 }
03158 }
03159 else if(mPi+nreM<tmM)
03160 {
03161 sMass=mPi;
03162 sPDG=-211;
03163 curQ+=tmpSQC+K0QC-PiMQC;
03164 s4Mom=G4LorentzVector(0.,0.,0.,sMass);
03165 reMass=nreM;
03166 rPDG=nucM.GetPDG();
03167 r4Mom=G4LorentzVector(0.,0.,0.,reMass);
03168 sum=reMass+sMass;
03169 if(fabs(tmM-sum)<eps)
03170 {
03171 r4Mom=q4Mom*(reMass/sum);
03172 s4Mom=q4Mom*(sMass/sum);
03173 }
03174 else if(tmM<sum || !G4QHadron(q4Mom).DecayIn2(r4Mom, s4Mom))
03175 {
03176 G4ExceptionDescription ed;
03177 ed << "(L->n)+Pi- DecayIn2: LamPiM Cor M=" << tmM << "=>rPDG=" << rPDG
03178 << "(rM=" << reMass << ")+sPDG=" << sPDG << "(sM=" << sMass << ")=" << sum
03179 << G4endl;
03180 G4Exception("G4Quasmon::HadronizeQuasmon()", "HAD_CHPS_0041",
03181 FatalException, ed);
03182 }
03183 }
03184 else if(nreM<tmM)
03185 {
03186 sMass=0.;
03187 sPDG=22;
03188 curQ+=tmpSQC+K0QC;
03189 s4Mom=G4LorentzVector(0.,0.,0.,sMass);
03190 reMass=nreZ;
03191 rPDG=nucZ.GetPDG();
03192 r4Mom=G4LorentzVector(0.,0.,0.,reMass);
03193 sum=reMass+sMass;
03194 if(fabs(tmM-sum)<eps)
03195 {
03196 r4Mom=q4Mom*(reMass/sum);
03197 s4Mom=q4Mom*(sMass/sum);
03198 }
03199 else if(tmM<sum || !G4QHadron(q4Mom).DecayIn2(r4Mom, s4Mom))
03200 {
03201 G4ExceptionDescription ed;
03202 ed << "(L->n)+Gamma DecayIn2: LamNGam Cor M=" << tmM << "=>rPDG=" << rPDG
03203 << "(rM=" << reMass << ")+sPDG=" << sPDG << "(sM=" << sMass << ")=" << sum
03204 << G4endl;
03205 G4Exception("G4Quasmon::HadronizeQuasmon()", "HAD_CHPS_0042",
03206 FatalException, ed);
03207 }
03208 }
03209 else
03210 {
03211 G4ExceptionDescription ed;
03212 ed << "LamTo0N with Pi DecayIn2: LamToN M=" << tmM << totQC << "=>rM="
03213 << nucM.GetPDG() << "," << nucZ.GetPDG() << "(" << nreM << "," << nreZ
03214 << ")+PiM/PiZ=" << mPi+nreM << "," << mPi0+nreZ << G4endl;
03215 G4Exception("G4Quasmon::HadronizeQuasmon()", "HAD_CHPS_0043",
03216 FatalException, ed);
03217 }
03218 }
03219 else if(tmM>iniQM)
03220 {
03221 G4QContent tmpSQC=G4QPDGCode(sPDG).GetQuarkContent();
03222 sMass=0.;
03223 sPDG=22;
03224 s4Mom=G4LorentzVector(0.,0.,0.,sMass);
03225 curQ+=tmpSQC;
03226 reMass=iniQM;
03227 rPDG=iniPDG;
03228 r4Mom=G4LorentzVector(0.,0.,0.,reMass);
03229 sum=reMass+sMass;
03230 if(fabs(tmM-reMass)<eps)
03231 {
03232 r4Mom=q4Mom*(reMass/sum);
03233 s4Mom=q4Mom*(sMass/sum);
03234 }
03235 else if(!G4QHadron(q4Mom).DecayIn2(r4Mom, s4Mom))
03236 {
03237 G4ExceptionDescription ed;
03238 ed << "QHadron+Gamma DecayIn2: gam+TQ M=" << tmM << "=>rPDG=" << rPDG << "(rM="
03239 << reMass <<")+sPDG="<< sPDG << "(sM=" << sMass << ")=" << reMass << G4endl;
03240 G4Exception("G4Quasmon::HadronizeQuasmon()", "HAD_CHPS_0044",
03241 FatalException, ed);
03242 }
03243 }
03244 else if(totMass>resTNM+sMass)
03245 {
03246 G4LorentzVector re4M = G4LorentzVector(0.,0.,0.,resTNM);
03247 G4LorentzVector rs4M = G4LorentzVector(0.,0.,0.,sMass);
03248 #ifdef debug
03249 G4cout<<"G4Q::HQ:EMERGENCY,rEM="<<resTN<<resTNM<<",fM="<<sMass<<",tM="<<totMass
03250 <<",d="<<totMass-resTNM-sMass<<G4endl;
03251 #endif
03252 sum=resTNM+sMass;
03253 if(fabs(totMass-sum)<eps)
03254 {
03255 re4M=tot4M*(resTNM/sum);
03256 rs4M=tot4M*(sMass/sum);
03257 }
03258 else if(totMass<sum || !G4QHadron(tot4M).DecayIn2(re4M,rs4M))
03259 {
03260 G4ExceptionDescription ed;
03261 ed << "DecayIn2 Frag+ResE failed: HadrQ:Decay T=" << totMass
03262 << "->R=" << resTNM << "+S=" << sMass << ")=" << sum << G4endl;
03263 G4Exception("G4Quasmon::HadronizeQuasmon()","HAD_CHPS_0045",FatalException,ed);
03264 }
03265 else
03266 {
03267
03268 G4QHadron* fragH = new G4QHadron(sPDG,rs4M);
03269 FillHadronVector(fragH);
03270 if(nQuasms==1)
03271 {
03272 resTN.Set4Momentum(re4M);
03273 qEnv=resTN;
03274 }
03275 else
03276 {
03277 G4QHadron* envH = new G4QHadron(resNQC,re4M);
03278 FillHadronVector(envH);
03279 qEnv = vacuum;
03280 }
03281 ClearQuasmon();
03282 return theQHadrons;
03283 }
03284 }
03285 else if(totMass>totM)
03286 {
03287 G4LorentzVector re4M = G4LorentzVector(0.,0.,0.,totM);
03288 G4LorentzVector rs4M = G4LorentzVector(0.,0.,0.,0.);
03289 #ifdef debug
03290 G4cout<<"G4Q::HQ:EMERGENSY,minM="<<totM<<" < totM="<<totMass<<G4endl;
03291 #endif
03292 if(fabs(totMass-totM)<eps) re4M=tot4M*(resTNM/sum);
03293 else if(!G4QHadron(tot4M).DecayIn2(re4M,rs4M))
03294 {
03295 G4ExceptionDescription ed;
03296 ed<<"DecayIn2 gam+TotN failed:HadrQ:Decay,T="<<totMass<<"->g+M="<<totM<<G4endl;
03297 G4Exception("G4Quasmon::HadronizeQuasmon()", "HAD_CHPS_0046",
03298 FatalException, ed);
03299
03300 }
03301 else
03302 {
03303 G4QHadron* fragH = new G4QHadron(22,rs4M);
03304 FillHadronVector(fragH);
03305 if(nQuasms==1)
03306 {
03307 totN.Set4Momentum(re4M);
03308 qEnv=totN;
03309 }
03310 else
03311 {
03312 G4QHadron* envH = new G4QHadron(totPDG,re4M);
03313 FillHadronVector(envH);
03314 qEnv = vacuum;
03315 }
03316 ClearQuasmon();
03317 return theQHadrons;
03318 }
03319 }
03320 else
03321 {
03322 G4cerr<<"***G4Q::HQ:M="<<tmM<<"=>rPDG="<<rPDG<<"(rM="<<reMass<<")+sPDG="
03323 <<sPDG<<"(sM="<<sMass<<")="<<sum<<",QM="<<iniQM<<G4endl;
03324 if(fabs(tmM-sum)<1.)
03325
03326 {
03327 r4Mom=q4Mom*(reMass/sum);
03328 s4Mom=q4Mom*(sMass/sum);
03329 }
03330
03331 else G4Exception("G4Quasmon::HadronizeQuasmon()", "HAD_CHPS_0047",
03332 FatalException, "QHadr+SHadr DecayIn2");
03333 }
03334 }
03335 G4double sKE=s4Mom.e()-sMass;
03336 #ifdef rdebug
03337 G4cout<<"G4Q::HQ:=2.3=>QHVect s4M="<<s4Mom<<",sPDG="<<sPDG<<", r4M/M="<<r4Mom<<reMass
03338 <<",fR="<<freeRQM<<",fW="<<fRQW<<",PDG="<<rPDG<<",r="<<rCB<<",s="<<sCB<<G4endl;
03339 #endif
03341 //if(sKE<sCB||rKE<rCB) // => "KinEn is below CB, try once more" case
03342 if(sKE<sCB)
03343 {
03344 #ifdef pdebug
03345 G4cout<<"****G4Q::HQ:E-9: sKE="<<sKE<<"<sCB="<<sCB<<G4endl;
03346
03347 G4Exception("G4Quasmon::HadronizeQuasmon()", "HAD_CHPS_0048",
03348 FatalException, "Why Fail? (9)");
03349 #endif
03350 if(sPDG>MINPDG) q4Mom-=G4LorentzVector(0.,0.,0.,pMass);
03351 qEnv=theEnvironment;
03352 return theQHadrons;
03353 }
03354 else if(abs(reMass-freeRQM)<fRQW||envPDG==NUCPDG)
03355 {
03356 G4QHadron* curHadr1 = new G4QHadron(rPDG,r4Mom);
03357 FillHadronVector(curHadr1);
03358 G4QHadron* curHadr2 = new G4QHadron(sPDG,s4Mom);
03359 FillHadronVector(curHadr2);
03360 #ifdef rdebug
03361 G4cout<<"G4Q::HQ:DecayQuasmon "<<q4Mom<<" in 4M="<<r4Mom+s4Mom<<" RQ="<<rPDG<<r4Mom
03362 <<" + Fragment="<<sPDG<<s4Mom<<", Env="<<theEnvironment<<G4endl;
03363 #endif
03364 if(sPDG>MINPDG) theEnvironment.Reduce(pPDG);
03365 ClearQuasmon();
03366 qEnv=theEnvironment;
03367 return theQHadrons;
03368 }
03369 else
03370 {
03371 G4LorentzVector resTotN4Mom=r4Mom+G4LorentzVector(0.,0.,0.,envM);
03372 G4QContent resTotNQC=envQC+curQ;
03373 G4QNucleus resTotN(resTotNQC);
03375 G4double resTotNM=resTotN.GetMZNS();
03376 if(resTotN4Mom.m()<resTotNM)
03377 {
03378
03379 #ifdef pdebug
03380 G4ExceptionDescription ed;
03381 ed<<"Why Fail?(10):NEEDS-EVAP-10,M="<<resTotN4Mom.m()<<"<miM="<<resTotNM<<G4endl;
03382 G4Exception("G4Quasmon::HadronizeQuasmon()","HAD_CHPS_0049", FatalException, ed);
03383 #endif
03384 if(sPDG>MINPDG) q4Mom-=G4LorentzVector(0.,0.,0.,pMass);
03385 qEnv=theEnvironment;
03386 return theQHadrons;
03387 }
03388 else
03389 {
03390 G4QHadron* curHadr2 = new G4QHadron(sPDG,s4Mom);
03391 FillHadronVector(curHadr2);
03392 q4Mom = r4Mom;
03393 if(sPDG>MINPDG)
03394 {
03395 theEnvironment.Reduce(pPDG);
03396 valQ += transQC;
03397 }
03398 else valQ = curQ;
03399 #ifdef rdebug
03400 G4cout<<"OK***>G4Q::HQ:S="<<sPDG<<s4Mom<<",Env="<<theEnvironment<<",Q="<<q4Mom
03401 <<valQ<<curQ<<G4endl;
03402 #endif
03403 status=1;
03404 phot4M=zeroLV;
03405 piF=false;
03406 gaF=false;
03407 if(CheckGroundState()) ClearQuasmon();
03408
03409 #ifdef rdebug
03410 G4cout<<"***>G4Q::HQ:After,S="<<sPDG<<s4Mom<<",Env="<<theEnvironment<<",Q="
03411 <<q4Mom<<valQ<<curQ<<G4endl;
03412 #endif
03413 qEnv=theEnvironment;
03414 return theQHadrons;
03415 }
03416 }
03417 }
03418 #ifdef rdebug
03419 else G4cout<<"G4Q::HQ:NO-OK,h="<<hsflag<<",d="<<fdul<<",M="<<rMass<<"<"<<reMass<<",M2="
03420 <<reTNM2<<"<I="<<tmpTM2<<",sP="<<sPDG<<",eP="<<envPDG<<",pP="<<pPDG<<G4endl;
03421 #endif
03422 if(!fskip)
03423 {
03424
03425 #ifdef debug
03426 G4int ePDG=theEnvironment.GetPDG();
03427 G4double frKin=fr4Mom.e()-sMass;
03428 G4cout<<"G4Q::HQ:>>"<<sPDG<<fr4Mom<<fr4Mom.m()<<"="<<sMass<<",T="<<frKin<<",E="<<ePDG
03429 <<G4endl;
03430 #endif
03431
03432
03433 if(sPDG<MINPDG)
03434 {
03435 G4int SQ=totQC.GetStrangeness();
03436 #ifdef debug
03437 G4cout<<"G4Q::HQ: sPDG="<<sPDG<<", sM="<<sMass<<", SQ="<<SQ<<G4endl;
03438 #endif
03439 if(!sPDG&&SQ<0&&nQuasms==1)
03440
03441 {
03442 sPDG=321;
03443 sMass=mK;
03444 G4QContent resKPQC=totQC-G4QContent(0,1,0,0,0,1);
03445 G4QNucleus rKPN(resKPQC);
03446 G4double rKPM = rKPN.GetMZNS();
03447 G4int rKPPDG = rKPN.GetPDG();
03448 G4QContent resK0QC=totQC-G4QContent(1,0,0,0,0,1);
03449 G4QNucleus rK0N(resK0QC);
03450 G4int rK0PDG = rK0N.GetPDG();
03451 G4double rK0M = rK0N.GetMZNS();
03452 if ( (rKPM+mK > totMass && rK0M+mK0 > totMass) ||
03453 rKPPDG == NUCPDG ||
03454 rK0PDG == NUCPDG )
03455 {
03456 #ifdef pdebug
03457 G4ExceptionDescription ed;
03458 ed << "Why PANIC? (2): ***PANIC#2***tM=" << totMass << "<KM=" << mK << ","
03459 << mK0 << ",rM=" << rKPM << "," << rK0M << ",d=" << mK+rKPM-totMass << ","
03460 << mK0+rK0M-totMass << G4endl;
03461 G4Exception("G4Quasmon::HadronizeQuasmon()","HAD_CHPS_0050",FatalException,ed);
03462 #endif
03463 status =-1;
03464 qEnv=theEnvironment;
03465 return theQHadrons;
03466 }
03467 if(rKPM + mK > rK0M + mK0)
03468 {
03469 rPDG = rK0PDG;
03470 rMass = rK0M;
03471 sPDG = 311;
03472 sMass = mK0;
03473 }
03474 else
03475 {
03476 rPDG = rKPPDG;
03477 rMass = rKPM;
03478 sPDG = 321;
03479 sMass = mK;
03480 }
03481 G4double ctM=tot4M.m();
03482 G4LorentzVector r4Mom(0.,0.,0.,rMass);
03483 G4LorentzVector s4Mom(0.,0.,0.,sMass);
03484 G4double sum=rMass+sMass;
03485 if(fabs(ctM-sum)<eps)
03486 {
03487 r4Mom=tot4M*(rMass/sum);
03488 s4Mom=tot4M*(sMass/sum);
03489 }
03490 else if(ctM<sum || !G4QHadron(tot4M).DecayIn2(r4Mom, s4Mom))
03491 {
03492 G4ExceptionDescription ed;
03493 ed << "HadrQuasm:K+ResNuc DecayIn2 didn't succeed: tM=" << ctM
03494 << totQC << " => rPDG=" << rPDG << "(rM=" << rMass << ") + sPDG="
03495 << sPDG << "(sM=" << sMass << ")=" << sum << G4endl;
03496 G4Exception("G4Quasmon::HadronizeQuasmon()","HAD_CHPS_0051",FatalException,ed);
03497 }
03498 #ifdef debug
03499 G4cout<<"G4Q::HQ:===2.4===>HadrVec s="<<sPDG<<s4Mom<<",r="<<rPDG<<r4Mom<<G4endl;
03500 #endif
03501
03502 G4QHadron* curHadr1 = new G4QHadron(rPDG,r4Mom);
03503 FillHadronVector(curHadr1);
03504 G4QHadron* curHadr2 = new G4QHadron(sPDG,s4Mom);
03505 FillHadronVector(curHadr2);
03506 ClearQuasmon();
03507 qEnv=vacuum;
03508 return theQHadrons;
03509 }
03510 G4bool ffin=false;
03511 if(quasM<rMass+sMass&&(sPDG==221||sPDG==331))
03512 {
03513 sPDG = 111;
03514 sMass=mPi0;
03515 }
03516 else if(!sPDG)
03517 {
03518 if (iniS<0&&iniQChg+iniQChg>=iniBN)
03519 {
03520 sPDG = 321;
03521 sMass= mK;
03522 G4QNucleus totQN(valQ+KpQC);
03523 rPDG = totQN.GetPDG();
03524 rMass= totQN.GetMZNS();
03525 }
03526 else if(iniS<0)
03527 {
03528 sPDG = 311;
03529 sMass= mK0;
03530 G4QNucleus totQN(valQ+K0QC);
03531 rPDG = totQN.GetPDG();
03532 rMass= totQN.GetMZNS();
03533 }
03534 else if(iniQChg>iniBN-iniS)
03535 {
03536 sPDG = 211;
03537 sMass= mPi;
03538 G4QNucleus totQN(valQ-PiQC);
03539 rPDG = totQN.GetPDG();
03540 rMass= totQN.GetMZNS();
03541 }
03542 else if(iniQChg<0)
03543 {
03544 sPDG = -211;
03545 sMass= mPi;
03546 G4QNucleus totQN(valQ+PiQC);
03547 rPDG = totQN.GetPDG();
03548 rMass= totQN.GetMZNS();
03549 }
03550 else if(quasM>iniQM+mPi0)
03551 {
03552 sPDG = 111;
03553 sMass= mPi0;
03554 rPDG = iniPDG;
03555 rMass= iniQM;
03556 }
03557 else
03558 {
03559 sPDG = 22;
03560 sMass= 0.;
03561 rPDG = iniPDG;
03562 rMass= iniQM;
03563 }
03564 ffin = true;
03565 }
03566 #ifdef debug
03567 G4cout<<"G4Q::HQ:MQ="<<q4Mom.m()<<"->sPDG="<<sPDG<<"(M="<<sMass<<") + rPDG="<<rPDG
03568 <<"(M="<<rMass<<")"<<",S="<<rMass+sMass<<G4endl;
03569 #endif
03570 if(q4Mom.m()+.003<rMass+sMass)
03571 {
03572 #ifdef debug
03573 G4cerr<<"G4Q::HQ:***PANIC#3***tM="<<q4Mom.m()<<"<rM="<<rMass<<",sM="<<sMass
03574 <<",d="<<rMass+sMass-q4Mom.m()<<G4endl;
03575 #endif
03576 #ifdef pdebug
03577
03578 G4Exception("G4Quasmon::HadronizeQuasmon()", "HAD_CHPS_0052",
03579 FatalException, "Why PANIC? (3)");
03580 #endif
03581 status =-1;
03582 qEnv=theEnvironment;
03583 return theQHadrons;
03584 }
03585 G4double cqM=q4Mom.m();
03586 G4LorentzVector resQ4Mom(0.,0.,0.,rMass);
03587 G4LorentzVector s4Mom(0.,0.,0.,sMass);
03588 G4double sum=rMass+sMass;
03589 if(fabs(cqM-sum)<eps)
03590 {
03591 resQ4Mom=q4Mom*(rMass/sum);
03592 s4Mom=q4Mom*(sMass/sum);
03593 }
03594 else if(cqM<sum || !G4QHadron(q4Mom).DecayIn2(resQ4Mom, s4Mom))
03595 {
03596
03597
03598
03599 G4ExceptionDescription ed;
03600 ed << "Quasm+Hadr DecayIn2 error: MQ=" << cqM << "-> rPDG=" << rPDG
03601 << ", (M=" << rMass << ") + sPDG=" << sPDG << "(M=" << sMass
03602 << ")=" << sum << G4endl;
03603 G4Exception("G4Quasmon::HadronizeQuasmon()", "HAD_CHPS_0053",
03604 FatalException, ed);
03605 }
03606 #ifdef debug
03607 G4cout<<"G4Q::HQ:Decay of Quasmon="<<q4Mom<<"->s="<<sPDG<<s4Mom<<"+R="<<resQ4Mom
03608 <<",f="<<ffin<<G4endl;
03609 #endif
03610 G4QHadron* candHadr = new G4QHadron(sPDG,s4Mom);
03611 if(ffin)
03612 {
03613
03614 theQHadrons.push_back(candHadr);
03615 G4QHadron* candHRes = new G4QHadron(rPDG,resQ4Mom);
03616 FillHadronVector(candHRes);
03617 ClearQuasmon();
03618 qEnv=theEnvironment;
03619 return theQHadrons;
03620 }
03621 else
03622 {
03623 G4QContent outQC = G4QPDGCode(sPDG).GetQuarkContent();
03624 G4int outChg = outQC.GetCharge();
03625 G4double outProb = 1.;
03626 if(theEnvironment.GetPDG()>NUCPDG)
03627 {
03628 G4int outBar = outQC.GetBaryonNumber();
03629 G4double outCB = theEnvironment.CoulombBarrier(outChg,outBar);
03630
03631 G4double outT = s4Mom.e()-s4Mom.m();
03632 outProb = theEnvironment.CoulBarPenProb(outCB,outT,outChg,outBar);
03633 }
03634 G4double rnd=G4UniformRand();
03635 #ifdef debug
03636 G4cout<<"G4Q::HQ: for "<<sPDG<<", rnd="<<rnd<<" < outP="<<outProb<<" ?"<<G4endl;
03637 #endif
03638 if(rnd<outProb)
03639 {
03640 FillHadronVector(candHadr);
03641 check+= s4Mom;
03642 ccheck+=outChg;
03643 q4Mom = resQ4Mom;
03644 valQ = curQ;
03645 status= 1;
03646 phot4M=zeroLV;
03647 piF=false;
03648 gaF=false;
03649 }
03650 else
03651 {
03652 status=3;
03653 delete candHadr;
03654 }
03655 }
03656 }
03657 else
03658 {
03659
03660 G4QContent outQC = G4QPDGCode(sPDG).GetQuarkContent();
03661 G4int outBar = outQC.GetBaryonNumber();
03662 G4int outChg = outQC.GetCharge();
03663 G4double outCB = theEnvironment.CoulombBarrier(outChg,outBar);
03664
03665 if(nucflag) rQ4Mom+=G4LorentzVector(dMass);
03666 G4QHadron tmpRQH(valQ+transQC,rQ4Mom);
03667
03668 G4double outT = fr4Mom.e()-fr4Mom.m();
03669 G4double outProb = theEnvironment.CoulBarPenProb(outCB,outT,outChg,outBar);
03670 if(G4UniformRand()<outProb)
03671 {
03672 theEnvironment.Reduce(pPDG);
03673 G4LorentzVector sumL=theEnvironment.Get4Momentum()+q4Mom;
03674 check += fr4Mom;
03675 ccheck+=G4QPDGCode(sPDG).GetCharge();
03676
03677
03678
03679
03680
03681
03682
03683
03684
03685
03686
03687
03688
03689
03690
03691
03692
03693
03694 q4Mom = rQ4Mom;
03695 if(sPDG>MINPDG) valQ += transQC;
03696 G4QHadron* candHadr = new G4QHadron(sPDG,fr4Mom);
03697 FillHadronVector(candHadr);
03698 #ifdef debug
03699 G4cout<<"G4Q::HQ:QuarkExchHadronizThroughCB Q="<<valQ<<",trQC="<<transQC<<G4endl;
03700 #endif
03701 sumL-=theEnvironment.Get4Momentum()+q4Mom+fr4Mom;
03702 #ifdef debug
03703 G4cout<<"G4Q::HQ:status=1, ------>> NuclearMatter SUBCHECK ---->>"<<sumL<<G4endl;
03704 #endif
03705 status=1;
03706 phot4M=zeroLV;
03707 piF=false;
03708 gaF=false;
03709 }
03710 else
03711 {
03712 #ifdef debug
03713 G4cout<<"G4Q::HQ:CBIsn'tPEN,P="<<outProb<<",T="<<outT<<",M="<<fr4Mom.m()<<G4endl;
03714 #endif
03715
03716
03717
03718 status=3;
03719 if(gaF)
03720 {
03721 phot4M=zeroLV;
03722 gaF=false;
03723 }
03724
03725 }
03726 }
03727 }
03728
03729 if(CheckGroundState())
03730
03731 {
03732 ClearQuasmon();
03733 qEnv=theEnvironment;
03734 return theQHadrons;
03735 }
03736 G4LorentzVector sumLor=theEnvironment.Get4Momentum()+q4Mom+check;
03737 #ifdef debug
03738 G4int eZ = theEnvironment.GetZ();
03739 G4int sumC = eZ+valQ.GetCharge()+ccheck;
03740 G4int curPDG=valQ.GetSPDGCode();
03741 G4cout<<"G4Q::HQ:Z="<<eZ<<valQ<<"***>FinalCHECK***>>4M="<<sumLor<<",Ch="<<sumC<<G4endl;
03742 if(!curPDG) G4cout<<"***G4Q::HQ: Quasmon-Tripolino QC="<<valQ<<G4endl;
03743 G4cout<<"G4Q::HQ:=------=> ResidualQ 4M="<<q4Mom<<", QC="<<valQ<<G4endl;
03744 #endif
03745 }
03746 #ifdef chdebug
03747 G4int ecSum=theEnvironment.GetZ()+valQ.GetCharge();
03748 G4int nHe=theQHadrons.size();
03749 if(nHe) for(int ih=0; ih<nHe; ih++) ecSum+=theQHadrons[ih]->GetCharge();
03750 if(ecSum!=cSum)
03751 {
03752 G4cerr<<"***G4Q::HQ:C"<<cSum<<",c="<<ecSum<<",E="<<theEnvironment<<",Q="<<valQ<<G4endl;
03753 G4cerr<<":G4Q::HQ:*END*,oE="<<oldEnv<<"oQ="<<oldCQC<<",oN="<<oldNH<<",N="<<nHe<<G4endl;
03754 if(nHe) for(G4int h=0; h<nHe; h++)
03755 {
03756 G4QHadron* cH = theQHadrons[h];
03757 G4cerr<<"::G4Q::HQ:#h"<<h<<",C="<<cH->GetCharge()<<",P="<<cH->GetPDGCode()<<G4endl;
03758 }
03759 }
03760 #endif
03761 #ifdef debug
03762 G4cout<<"G4Q::HQ: Q="<<q4Mom<<valQ<<",E="<<theEnvironment<<", status="<<status<<G4endl;
03763 #endif
03764 qEnv=theEnvironment;
03765 return theQHadrons;
03766 }
03767
03768
03769 G4QHadronVector* G4Quasmon::DecayQuasmon()
03770 {
03771 G4QHadron* thisQuasmon = new G4QHadron(valQ,q4Mom);
03772 FillHadronVector(thisQuasmon);
03773 G4QHadronVector* theFragments = new G4QHadronVector;
03774 G4int nHadrs=theQHadrons.size();
03775 #ifdef debug
03776 G4cout<<"G4Q::DecayQuasmon:After decay (FillHadronVector byItself) nH="<<nHadrs<<G4endl;
03777 #endif
03778 if(nHadrs) for (int hadron=0; hadron<nHadrs; hadron++)
03779 {
03780 G4QHadron* curHadr = new G4QHadron(theQHadrons[hadron]);
03781 theFragments->push_back(curHadr);
03782 }
03783 #ifdef pdebug
03784 else G4cerr<<"*******G4Quasmon::DecayQuasmon: *** Nothing is in the output ***"<<G4endl;
03785 #endif
03786 valQ=G4QContent(0,0,0,0,0,0);
03787 q4Mom=G4LorentzVector(0.,0.,0.,0.);
03788 return theFragments;
03789 }
03790
03791
03792 void G4Quasmon::FillHadronVector(G4QHadron* qH)
03793 {
03794
03795
03796
03797
03798
03799
03800
03801
03802
03803
03804
03805
03806 static const G4LorentzVector zeroLV(0.,0.,0.,0.);
03807 static const G4double mAlph = G4QPDGCode(2112).GetNuclMass(2,2,0);
03808 static const G4QContent neutQC(2,1,0,0,0,0);
03809 static const G4QContent protQC(1,2,0,0,0,0);
03810 static const G4QContent sigmQC(2,0,1,0,0,0);
03811 static const G4QContent lambQC(1,1,1,0,0,0);
03812 static const G4QContent sigpQC(0,2,1,0,0,0);
03813 static const G4QContent PiQC(0,1,0,1,0,0);
03814 static const G4QContent K0QC(1,0,0,0,0,1);
03815 static const G4QContent KpQC(0,1,0,0,0,1);
03816 static const G4double mNeut= G4QPDGCode(2112).GetMass();
03817 static const G4double mProt= G4QPDGCode(2212).GetMass();
03818 static const G4double mSigM= G4QPDGCode(3112).GetMass();
03819 static const G4double mLamb= G4QPDGCode(3122).GetMass();
03820 static const G4double mSigP= G4QPDGCode(3222).GetMass();
03821 static const G4double mPi = G4QPDGCode(211).GetMass();
03822 static const G4double mPi0 = G4QPDGCode(111).GetMass();
03823 static const G4double mK = G4QPDGCode(321).GetMass();
03824 static const G4double mK0 = G4QPDGCode(311).GetMass();
03825
03826 status=1;
03827 phot4M=zeroLV;
03828 G4int thePDG = qH->GetPDGCode();
03829 G4LorentzVector t = qH->Get4Momentum();
03830 #ifdef psdebug
03831 if(thePDG==113 && fabs(t.m()-770.)<.001)
03832 {
03833 G4cerr<<"G4Q::FillHadronVector: PDG="<<thePDG<<",M="<<t.m()<<G4endl;
03834
03835 G4Exception("G4Quasmon::FillHadronVector()", "HAD_CHPS_0000",
03836 FatalException, "Zero rho");
03837 }
03838 #endif
03839 #ifdef pdebug
03840 G4cout<<"G4Q::FillHadronVector:Hadron's PDG="<<thePDG<<",4Mom="<<t<<",m="<<t.m()<<G4endl;
03841 #endif
03842 if(thePDG>80000000 && (thePDG<90000000 || thePDG%1000>500 || thePDG%1000000>500000)
03843 && thePDG!=90002999 && thePDG!=89999003 && thePDG!=90003998 && thePDG!=89998004
03844 && thePDG!=90003999 && thePDG!=89999004 && thePDG!=90004998 && thePDG!=89998005)
03845 {
03846 if (thePDG==90999999) thePDG=-311;
03847 else if(thePDG==90999000) thePDG=-321;
03848 else if(thePDG==89000001) thePDG=311;
03849 else if(thePDG==89001000) thePDG=321;
03850 else if(thePDG==90000999) thePDG=211;
03851 else if(thePDG==89999001) thePDG=-211;
03852 else if(thePDG==89999999) thePDG=-2112;
03853 else if(thePDG==89999000) thePDG=-2212;
03854 else if(thePDG==89000000) thePDG=-3122;
03855 else if(thePDG==88999002) thePDG=-3222;
03856 else if(thePDG==89000999) thePDG=-3222;
03857 else if(thePDG==88000001) thePDG=-3322;
03858 else if(thePDG==88001000) thePDG=-3312;
03859 else if(thePDG==89999002) thePDG=1114;
03860 else if(thePDG==90001999) thePDG=2224;
03861 else if(thePDG==91000000) thePDG=3122;
03862 else if(thePDG==90999001) thePDG=3112;
03863 else if(thePDG==91000999) thePDG=3222;
03864 else if(thePDG==91999000) thePDG=3312;
03865 else if(thePDG==91999999) thePDG=3322;
03866 else if(thePDG==92998999) thePDG=3112;
03867 #ifdef pdebug
03868 else G4cerr<<"*G4Quasmon::FillQHV:PDG="<<thePDG<<",M="<<qH->Get4Momentum().m()<<G4endl;
03869 #endif
03870 qH->SetQPDG(G4QPDGCode(thePDG));
03871 }
03872 if (thePDG==10)
03873 {
03874 G4double rM = t.m();
03875 G4QContent chipQC = qH->GetQC();
03876 G4QContent h1QC = chipQC.SplitChipo(rM);
03877 G4QContent h2QC = chipQC - h1QC;
03878 G4int h1PDG = h1QC.GetSPDGCode();
03879 G4int h2PDG = h2QC.GetSPDGCode();
03880 if(!h1PDG || !h2PDG)
03881 {
03882 G4cerr<<"***FillHV:h1QC="<<h1QC<<"(PDG="<<h1PDG<<"),h2QC="<<h2QC<<"(PDG="<<h2PDG<<")"
03883 <<G4endl;
03884 G4Exception("G4Quasmon::FillHadronVector()", "HAD_CHPS_0001", FatalException,
03885 "Chipolino can't be defragmented");
03886 }
03887 G4QHadron* fHadr = new G4QHadron(h1PDG);
03888 G4QHadron* sHadr = new G4QHadron(h2PDG);
03889 G4LorentzVector f4Mom = fHadr->Get4Momentum();
03890 G4LorentzVector s4Mom = sHadr->Get4Momentum();
03891 if(!qH->DecayIn2(f4Mom,s4Mom))
03892 {
03893 delete fHadr;
03894 delete sHadr;
03895 G4cerr<<"***G4Q::FillHadrV:ChipQC"<<chipQC<<":PDG1="<<h1PDG<<",PDG2="<<h2PDG<<G4endl;
03896 theQHadrons.push_back(qH);
03897 }
03898 else
03899 {
03900 delete qH;
03901 fHadr->Set4Momentum(f4Mom);
03902 FillHadronVector(fHadr);
03903 sHadr->Set4Momentum(s4Mom);
03904 FillHadronVector(sHadr);
03905 }
03906 }
03907 else if(thePDG>80000000&&thePDG!=90000000)
03908 {
03909 G4double fragMas=t.m();
03910
03911 G4QNucleus qNuc(t,thePDG);
03912
03913
03914 G4double GSMass = G4QPDGCode(thePDG).GetMass();
03915 G4QContent totQC=qNuc.GetQCZNS();
03916 G4int nN =qNuc.GetN();
03917 G4int nZ =qNuc.GetZ();
03918 G4int nS =qNuc.GetS();
03919 G4int bA =qNuc.GetA();
03920 #ifdef pdebug
03921 G4cout<<"G4Quasm::FillHadrVect:Nucl="<<qNuc<<",nPDG="<<thePDG<<",GSM="<<GSMass<<G4endl;
03922 #endif
03923 if((nN<0 || nZ<0 || nS<0) && bA>0)
03924 {
03925 G4double m1=mPi;
03926 G4int PDG1=-211;
03927 G4QNucleus newNpm(totQC+PiQC);
03928 G4int newS=newNpm.GetStrangeness();
03929 if(newS>0) newNpm=G4QNucleus(totQC+PiQC+newS*K0QC);
03930 G4int PDG2=newNpm.GetPDG();
03931 G4double m2_value=newNpm.GetMZNS();
03932 if(nS<0)
03933 {
03934 m1 =mK;
03935 PDG1 =321;
03936 G4QNucleus newNp(totQC-KpQC);
03937 PDG2 =newNp.GetPDG();
03938 m2_value =newNp.GetMZNS();
03939 G4QNucleus newN0(totQC-K0QC);
03940 G4double m3_value=newN0.GetMZNS();
03941 if (m3_value+mK0<m2_value+mK)
03942 {
03943 m1 =mK0;
03944 PDG1=311;
03945 m2_value =m3_value;
03946 PDG2=newN0.GetPDG();
03947 }
03948 }
03949 else if(nS>0&&nZ+nN>0)
03950 {
03951 if(nN<0)
03952 {
03953 m1 =mSigP;
03954 PDG1 =3222;
03955 G4QNucleus newNp(totQC-sigpQC);
03956 PDG2 =newNp.GetPDG();
03957 m2_value =newNp.GetMZNS();
03958 }
03959 else
03960 {
03961 m1 =mSigM;
03962 PDG1 =3112;
03963 G4QNucleus newNp(totQC-sigmQC);
03964 PDG2 =newNp.GetPDG();
03965 m2_value =newNp.GetMZNS();
03966 }
03967 }
03968 else if(nN<0)
03969 {
03970 PDG1 =211;
03971 G4QNucleus newNpp(totQC-PiQC);
03972 PDG2 =newNpp.GetPDG();
03973 m2_value =newNpp.GetMZNS();
03974 }
03975 if(fragMas>m1+m2_value)
03976 {
03977 G4LorentzVector fq4M(0.,0.,0.,m1);
03978 G4LorentzVector qe4M(0.,0.,0.,m2_value);
03979 if(!qH->DecayIn2(fq4M,qe4M))
03980 {
03981 G4ExceptionDescription ed;
03982 ed << "Mes+ResA DecayIn2 did not succeed: QM=" << t.m() << "-> Mes=" << PDG1
03983 << "(M=" << m1 << ") + ResA=" << PDG2 << "(M=" << m2_value << ")" << G4endl;
03984 G4Exception("G4Quasmon::FillHadronVector()","HAD_CHPS_0002", FatalException, ed);
03985 }
03986 delete qH;
03987 G4QHadron* H1 = new G4QHadron(PDG1,fq4M);
03988 theQHadrons.push_back(H1);
03989 G4QHadron* H2 = new G4QHadron(PDG2,qe4M);
03990 FillHadronVector(H2);
03991 }
03992 else if(fabs(m1+m2_value-fragMas)<0.01)
03993 {
03994 G4double r1=m1/fragMas;
03995 G4double r2=1.-r1;
03996
03997
03998
03999 delete qH;
04000
04001 G4QHadron* H1 = new G4QHadron(PDG1,r1*t);
04002 theQHadrons.push_back(H1);
04003 G4QHadron* H2 = new G4QHadron(PDG2,r2*t);
04004 FillHadronVector(H2);
04005 }
04006 else
04007 {
04008 #ifdef debug
04009 G4cerr<<"-Warning-G4Q::FillHVec:PDG="<<thePDG<<"("<<t.m()<<","<<fragMas<<") < Mes="
04010 <<PDG1<<"("<<m1<<") + ResA="<<PDG2<<"("<<m2_value<<"), d="<<fragMas-m1-m2_value<<G4endl;
04011
04012 #endif
04013 theQHadrons.push_back(qH);
04014 }
04015 }
04016 else if(abs(fragMas-GSMass)<.1)
04017 {
04018 #ifdef pdebug
04019 G4cout<<"G4Quasm::FillHadrVect: Ground state"<<G4endl;
04020 #endif
04021 G4double nResM =1000000.;
04022 G4int nResPDG=0;
04023 if(nN>0&&bA>1)
04024 {
04025 G4QContent resQC=totQC-neutQC;
04026 G4QNucleus resN(resQC);
04027 nResPDG=resN.GetPDG();
04028 if (nResPDG==90000001) nResM=mNeut;
04029 else if(nResPDG==90001000) nResM=mProt;
04030 else if(nResPDG==91000000) nResM=mLamb;
04031 else nResM=resN.GetMZNS();
04032 }
04033 G4double pResM =1000000.;
04034 G4int pResPDG=0;
04035 if(nZ>0&&bA>1)
04036 {
04037 G4QContent resQC=totQC-protQC;
04038 G4QNucleus resN(resQC);
04039 pResPDG=resN.GetPDG();
04040 if (pResPDG==90000001) pResM=mNeut;
04041 else if(pResPDG==90001000) pResM=mProt;
04042 else if(pResPDG==91000000) pResM=mLamb;
04043 else pResM =resN.GetMZNS();
04044 }
04045 G4double lResM =1000000.;
04046 G4int lResPDG=0;
04047 if(nS>0&&bA>1)
04048 {
04049 G4QContent resQC=totQC-lambQC;
04050 G4QNucleus resN(resQC);
04051 lResPDG=resN.GetPDG();
04052 if (lResPDG==90000001) lResM=mNeut;
04053 else if(lResPDG==90001000) lResM=mProt;
04054 else if(lResPDG==91000000) lResM=mLamb;
04055 else lResM =resN.GetMZNS();
04056 }
04057 #ifdef debug
04058 G4cout<<"G4Quasm::FillHadrVec:rP="<<pResPDG<<",rN="<<nResPDG<<",rL="<<lResPDG<<",nN="
04059 <<nN<<",nZ="<<nZ<<",nL="<<nS<<",totM="<<fragMas<<",n="<<fragMas-nResM-mNeut
04060 <<",p="<<fragMas-pResM-mProt<<",l="<<fragMas-lResM-mLamb<<G4endl;
04061 #endif
04062 if ( thePDG == 90004004 ||
04063 (bA > 1 && ( (nN > 0 && fragMas > nResM+mNeut) ||
04064 (nZ > 0 && fragMas > pResM+mProt) ||
04065 (nS > 0 && fragMas > lResM+mLamb) ) ) )
04066 {
04067 G4int barPDG = 90002002;
04068 G4int resPDG = 90002002;
04069 G4double barM= mAlph;
04070 G4double resM= mAlph;
04071
04072 if (fragMas > nResM+mNeut) {
04073 barPDG = 90000001;
04074 resPDG = nResPDG;
04075 barM= mNeut;
04076 resM= nResM;
04077 }
04078 else if(fragMas>pResM+mProt)
04079 {
04080 barPDG=90001000;
04081 resPDG=pResPDG;
04082 barM =mProt;
04083 resM =pResM;
04084 }
04085 else if(fragMas>lResM+mLamb)
04086 {
04087 barPDG=91000000;
04088 resPDG=lResPDG;
04089 barM =mLamb;
04090 resM =lResM;
04091 }
04092 else if(thePDG!=90004004 && fragMas>GSMass)
04093 {
04094 barPDG=22;
04095 resPDG=thePDG;
04096 barM =0.;
04097 resM =pResM;
04098 }
04099 else if(thePDG!=90004004)
04100 {
04101 G4ExceptionDescription ed;
04102 ed << "Below GSM but cann't decay: PDG=" << thePDG << ",M=" << fragMas
04103 << "<GSM=" << GSMass << G4endl;
04104 G4Exception("G4Quasmon::FillHadronVector()","HAD_CHPS_0003", FatalException, ed);
04105 }
04106 G4LorentzVector a4Mom(0.,0.,0.,barM);
04107 G4LorentzVector b4Mom(0.,0.,0.,resM);
04108 if(!qH->DecayIn2(a4Mom,b4Mom))
04109 {
04110 theQHadrons.push_back(qH);
04111 G4cerr<<"---Warning---G4Q::FillHadronVector: Be8 decay did not succeed"<<G4endl;
04112 }
04113 else
04114 {
04115
04116
04117
04118 delete qH;
04119
04120 G4QHadron* HadrB = new G4QHadron(barPDG,a4Mom);
04121 FillHadronVector(HadrB);
04122 G4QHadron* HadrR = new G4QHadron(resPDG,b4Mom);
04123 FillHadronVector(HadrR);
04124 }
04125 }
04126 else
04127 {
04128 #ifdef debug
04129 G4cout<<"G4Quasm::FillHadrVect: Leave as it is"<<G4endl;
04130 #endif
04131 theQHadrons.push_back(qH);
04132 }
04133 }
04134 else if (fragMas < GSMass)
04135 {
04136 G4cerr<<"***G4Quasmon::FillHV:M="<<fragMas<<">GSM="<<GSMass<<"(PDG="<<thePDG<<"),d="
04137 <<fragMas-GSMass<<", NZS="<<nN<<","<<nZ<<","<<nS<<G4endl;
04138
04139 G4cout<<"****>>G4Quasm::FillHadrVect: Leave as it is Instead of Exception"<<G4endl;
04140 theQHadrons.push_back(qH);
04141 }
04142 else if (bA==1 && fragMas>GSMass)
04143 {
04144 G4int gamPDG=22;
04145 G4double gamM=0.;
04146 if(fragMas>mPi0+GSMass)
04147 {
04148 gamPDG=111;
04149 gamM=mPi0;
04150 }
04151 G4LorentzVector a4Mom(0.,0.,0.,gamM);
04152 G4LorentzVector b4Mom(0.,0.,0.,GSMass);
04153 if(!qH->DecayIn2(a4Mom,b4Mom))
04154 {
04155 theQHadrons.push_back(qH);
04156 G4cerr<<"---Warning---G4Q::FillHadrVect:N*->gamma/pi0+N decay error"<<G4endl;
04157 }
04158 else
04159 {
04160 G4QHadron* HadrB = new G4QHadron(gamPDG,a4Mom);
04161 FillHadronVector(HadrB);
04162 qH->Set4Momentum(b4Mom);
04163 theQHadrons.push_back(qH);
04164 }
04165 }
04166 else
04167 {
04168 #ifdef ppdebug
04169 G4cout<<"G4Quasm::FillHadrVect:Evaporate "<<thePDG<<",tM="<<fragMas<<" > GS="<<GSMass
04170 <<qNuc.Get4Momentum()<<", m="<<qNuc.Get4Momentum().m()<<G4endl;
04171 #endif
04172 G4QHadron* bHadron = new G4QHadron;
04173 G4QHadron* rHadron = new G4QHadron;
04174 if(!qNuc.EvaporateBaryon(bHadron,rHadron))
04175 {
04176 G4cerr<<"---Warning---G4Q::FillHV:Evaporate PDG="<<thePDG<<",M="<<fragMas<<G4endl;
04177 delete bHadron;
04178 delete rHadron;
04179 theQHadrons.push_back(qH);
04180 }
04181 else
04182 {
04183 #ifdef debug
04184 G4cout<<"G4Q::FlHV:Done,b="<<bHadron->GetQPDG()<<",r="<<rHadron->GetQPDG()<<G4endl;
04185 #endif
04186 delete qH;
04187 if(bHadron->GetPDGCode() < -1111)
04188 {
04189 G4QHadronVector* tmpQHadVec=DecayQHadron(bHadron);
04190 G4int tmpS=tmpQHadVec->size();
04191 theQHadrons.resize(tmpS+theQHadrons.size());
04192 copy( tmpQHadVec->begin(), tmpQHadVec->end(), theQHadrons.end()-tmpS);
04193 tmpQHadVec->clear();
04194 delete tmpQHadVec;
04195 }
04196 else FillHadronVector(bHadron);
04197 if(rHadron->GetPDGCode() < -1111)
04198 {
04199 G4QHadronVector* tmpQHadVec=DecayQHadron(rHadron);
04200 G4int tmpS=tmpQHadVec->size();
04201 theQHadrons.resize(tmpS+theQHadrons.size());
04202 copy( tmpQHadVec->begin(), tmpQHadVec->end(), theQHadrons.end()-tmpS);
04203 tmpQHadVec->clear();
04204 delete tmpQHadVec;
04205 }
04206 else FillHadronVector(rHadron);
04207 }
04208 }
04209 }
04210 else
04211 {
04212 #ifdef pdebug
04213 G4cout<<"G4Q::FillHV: ---DECAY--- QH="<<qH->GetPDGCode()<<qH->Get4Momentum()<<G4endl;
04214 #endif
04215 G4QHadronVector* tmpQHadVec=DecayQHadron(qH);
04216 #ifdef pdebug
04217 G4cout<<"G4Q::FillHV: ---DECAY IS DONE--- with nH="<<tmpQHadVec->size()<<G4endl;
04218 #endif
04219 G4int tmpS=tmpQHadVec->size();
04220 theQHadrons.resize(tmpS+theQHadrons.size());
04221 copy( tmpQHadVec->begin(), tmpQHadVec->end(), theQHadrons.end()-tmpS);
04222 #ifdef pdebug
04223 G4cout<<"G4Q::FillHV: -->Products are added to QHV, nQHV="<<theQHadrons.size()<<G4endl;
04224 #endif
04225 tmpQHadVec->clear();
04226 delete tmpQHadVec;
04227 #ifdef pdebug
04228 G4cout<<"G4Q::FillHV: TemporaryQHV of DecayProducts is deleted"<<G4endl;
04229 #endif
04230 }
04231 }
04232
04233
04234
04235
04236 G4double G4Quasmon::GetQPartonMomentum(G4double kMax, G4double mC2)
04237 {
04238
04239 #ifdef debug
04240
04241 G4cout<<"G4Quas::GetQPartonMomentum:***called*** kMax="<<kMax<<",mC="<<sqrt(mC2)<<G4endl;
04242 #endif
04243 G4double qMass = q4Mom.m();
04244 G4double kLim = qMass/2.;
04245 G4double twM = qMass+qMass;
04246 G4double kMin = mC2/twM;
04247
04248
04249
04250
04251
04252
04253
04254
04255
04256
04257
04258
04259 if (kLim<kMax) kMax = kLim;
04260 if (kMin<0 || kMax<0 || qMass<=0. || nOfQ<2)
04261 {
04262 G4ExceptionDescription ed;
04263 ed << "Can not generate quark-parton: kMax=" << kMax << ", kMin=" << kMin
04264 << ", kLim=" << kLim << ", MQ=" << qMass << ", n=" << nOfQ << G4endl;
04265 G4Exception("G4Quasmon::GetQPartonMomentum()", "HAD_CHPS_0000", FatalException, ed);
04266 }
04267 #ifdef debug
04268 G4cout<<"G4Q::GetQPM: kLim="<<kLim<<",kMin="<<kMin<<",kMax="<<kMax<<",nQ="<<nOfQ<<G4endl;
04269 #endif
04270 if(kMin>kMax||nOfQ==2) return kMax;
04271 G4int n=nOfQ-2;
04272 G4double fn=n;
04273 G4int dn=n;
04274 G4double vRndm = G4UniformRand();
04275
04276
04277 if (kMin>0.)
04278 {
04279 G4double xMin=kMin/kLim;
04280 if (kMax>=kLim) vRndm = vRndm*pow((1.-xMin),n)*(1.+n*xMin);
04281 else
04282 {
04283 G4double xMax=kMax/kLim;
04284 G4double vRmin = pow((1.-xMin),n)*(1.+n*xMin);
04285 G4double vRmax = pow((1.-xMax),n)*(1.+n*xMax);
04286 vRndm = vRmax + vRndm*(vRmin-vRmax);
04287 }
04288 }
04289 else if (kMax<kLim)
04290 {
04291 G4double xMax=kMax/kLim;
04292 G4double vRmax = pow((1.-xMax),n)*(1.+n*xMax);
04293 vRndm = vRmax + vRndm*(1.-vRmax);
04294 }
04295 if (vRndm<=0. || vRndm>1.)
04296 {
04297
04298
04299 if(vRndm<=0.) vRndm=1.e-9;
04300 else if(vRndm>1.) vRndm=1.;
04301 }
04302 if (n==1) return kLim*sqrt(1.-vRndm);
04303 else
04304 {
04305 G4double x = 1.-pow(vRndm*(1+n*vRndm)/(fn+1.),1./fn);
04306 G4double ox = x;
04307 G4int it = 0;
04308 G4double d = 1.;
04309 G4double df = 1./static_cast<double>(nOfQ);
04310 G4double f = df*(static_cast<int>(nOfQ*nOfQ*n*x/5.)+(nOfQ/2));
04311 G4double xMin=.0001;
04312 G4double xMax=.9999;
04313 if(kLim>0)
04314 {
04315 xMin=kMin/kLim;
04316 xMax=kMax/kLim;
04317 }
04318 if(f>27.)
04319 {
04320 #ifdef debug
04321 G4cout<<"G4Q::GetQPMom: f="<<f<<" is changed to 99"<<G4endl;
04322 #endif
04323 f = 27.;
04324 }
04325 if(x<1.e-27) x=1.e-27;
04326 else if(x>.999999999) x=.999999999;
04327 G4double r = 1.-x;
04328 G4double p = r;
04329 if (n>2) p = pow(r,n-1);
04330 G4double nx = n*x;
04331 G4double c = p*r*(1.+nx);
04332 G4double od = c - vRndm;
04333 #ifdef debug
04334 G4cout<<"G4Q::GetQPMom:--->> First x="<<x<<", n="<<n<<", f="<<f<<", d/R(first)="
04335 <<od/vRndm<<G4endl;
04336 #endif
04337 G4int nitMax=dn+dn;
04338 if(nitMax>100)nitMax=100;
04339 while( abs(d/vRndm) > 0.001 && it <= nitMax)
04340 {
04341 x = x + f*od/(r*nx*(fn+1.));
04342 if(x<1.e-27) x=1.e-27;
04343 else if(x>.999999999) x=.999999999;
04344 r = 1.-x;
04345 if (n>2) p = pow(r,n-1);
04346 else p = r;
04347 nx = n*x;
04348 c = p*r*(1.+nx);
04349 d = c - vRndm;
04350 if ((od>0&&d<0)||(od<0&&d>0))
04351 {
04352 if (f>1.0001) f=1.+(f-1.)/2.;
04353 if (f<0.9999) f=1.+(1.-f)*2;
04354 x = (x + ox)/2.;
04355 if(x<1.e-27) x=1.e-27;
04356 else if(x>.999999999) x=.999999999;
04357 r = 1.-x;
04358 if (n>2) p = pow(r,n-1);
04359 else p = r;
04360 nx = n*x;
04361 c = p*r*(1.+nx);
04362 d = c - vRndm;
04363 }
04364 else
04365 {
04366 if (f>1.0001&&f<27.) f=1.+(f-1.)*2;
04367 if (f<0.99999999999) f=1.+(1.-f)/2.;
04368 if (f>=27.) f=27.;
04369 }
04370 #ifdef debug
04371 G4cout<<"G4Q::GetQPMom: Iter#"<<it<<": (c="<<c<<" - R="<<vRndm<<")/R ="<<d/vRndm
04372 <<", x="<<x<<", f="<<f<<G4endl;
04373 #endif
04374 if(x>xMax) x=xMax;
04375 if(x<xMin) x=xMin;
04376 if(fabs(d)>fabs(od) && n>99 && x!=xMin && x!=xMax)
04377 {
04378 x=ox;
04379 break;
04380 }
04381 od = d;
04382 ox = x;
04383 it++;
04384 }
04385 #ifdef debug
04386 if(it>nitMax) G4cout<<"G4Q::GetQPMom: a#of iterations > nitMax="<<nitMax<<G4endl;
04387 #endif
04388 if(x>xMax) x=xMax;
04389 if(x<xMin) x=xMin;
04390 G4double kCand=kLim*x;
04391 if(kCand>=kMax)kCand=kMax-.001;
04392 if(kCand<=kMin)kCand=kMin+.001;
04393 return kCand;
04394 }
04395
04396
04397
04398
04399
04400
04401
04402
04403
04404
04405
04406
04407
04408
04409
04410
04411
04412
04413
04414
04415
04417
04418
04419
04420 }
04421
04422
04423 G4int G4Quasmon::CalculateNumberOfQPartons(G4double qMass)
04424 {
04425 static const G4double mK0 = G4QPDGCode(311).GetMass();
04426
04427
04428
04429
04430
04431
04432 G4double qMOverT = qMass/Temperature;
04433 G4int valc = valQ.GetTot();
04434
04435
04441
04442
04443
04444
04446
04447
04448
04449
04450
04451
04452
04453
04454
04455
04456
04457
04458
04459 nOfQ = RandomPoisson((1.+sqrt(1.+qMOverT*qMOverT))/2.);
04460 G4int ev = valc%2;
04461 if (!ev && nOfQ<2) nOfQ=2;
04462 else if ( ev && nOfQ<3) nOfQ=3;
04463
04464 #ifdef debug
04465 G4cout<<"G4Q::Calc#ofQP:QM="<<q4Mom<<qMass<<",T="<<Temperature<<",QC="<<valQ<<",n="<<nOfQ
04466 <<G4endl;
04467 #endif
04468 G4int absb = abs(valQ.GetBaryonNumber());
04469 G4int tabn = 0;
04470 if(absb)tabn=3*absb;
04471 else tabn=4;
04472 if (nOfQ<tabn) nOfQ=tabn;
04473 G4int nSeaPairs = (nOfQ-valc)/2;
04474 G4int stran = abs(valQ.GetS());
04475 G4int astra = abs(valQ.GetAS());
04476 if(astra>stran) stran=astra;
04477 G4int nMaxStrangeSea=static_cast<int>((qMass-stran*mK0)/(mK0+mK0));
04478 if (absb) nMaxStrangeSea=static_cast<int>((qMass-absb)/672.);
04479 #ifdef debug
04480 G4cout<<"G4Q::Calc#ofQP:"<<valQ<<",INtot="<<valc<<",nOfQ="<<nOfQ<<",SeaPairs="<<nSeaPairs
04481 <<G4endl;
04482 #endif
04483 if (nSeaPairs)
04484 {
04485 #ifdef debug
04486 G4int morDec=0;
04487 #endif
04488 if(nSeaPairs>0)valQ.IncQAQ(nSeaPairs,SSin2Gluons);
04489 #ifdef debug
04490 else morDec=valQ.DecQAQ(-nSeaPairs);
04491 if(morDec) G4cout<<"G4Q::Calc#ofQP: "<<morDec<<" pairs can be reduced more"<<G4endl;
04492 #endif
04493 G4int sSea=valQ.GetS();
04494 G4int asSea=valQ.GetAS();
04495 if(asSea<sSea) sSea=asSea;
04496 if(sSea>nMaxStrangeSea)
04497 {
04498 #ifdef debug
04499 G4cout<<"G4Q::Calc#ofQP:**Reduce** S="<<sSea<<",aS="<<asSea<<",maxS="<<nMaxStrangeSea
04500 <<G4endl;
04501 #endif
04502 sSea-=nMaxStrangeSea;
04503 valQ.DecS(sSea);
04504 valQ.DecAS(sSea);
04505 valQ.IncQAQ(sSea,0.);
04506 }
04507 }
04508
04509
04510
04511
04512
04513 #ifdef debug
04514 G4cout<<"G4Quasmon::Calc#ofQP: *** RESULT IN*** nQ="<<nOfQ<<", FinalQC="<<valQ<<G4endl;
04515 #endif
04516 return nOfQ;
04517 }
04518
04519
04520 void G4Quasmon::ModifyInMatterCandidates()
04521 {
04523 G4double envM = theEnvironment.GetMass();
04524 G4QContent envQC=theEnvironment.GetQCZNS();
04525 G4int eP = theEnvironment.GetZ();
04526 G4int eN = theEnvironment.GetN();
04527 G4int eL = theEnvironment.GetS();
04528 G4QContent totQC=theEnvironment.GetQC()+valQ;
04529 G4int tP = totQC.GetP();
04530 G4int tN = totQC.GetN();
04531 G4int tL = totQC.GetL();
04532 G4double totM=G4QNucleus(totQC).GetMZNS();
04533 for (unsigned ind=0; ind<theQCandidates.size(); ind++)
04534 {
04535 G4QCandidate* curCand=theQCandidates[ind];
04536 G4int cPDG = curCand->GetPDGCode();
04537 G4bool poss = curCand->GetPossibility();
04538 G4QContent tmpTQ=totQC-curCand->GetQC();
04539 G4QNucleus tmpT(tmpTQ);
04540 G4double tmpTM=tmpT.GetMZNS();
04541 G4QPDGCode cQPDG(cPDG);
04542 G4double frM=cQPDG.GetMass();
04543 if(cPDG>80000000&&cPDG!=90000000)
04544 {
04545 if(totMass<tmpTM+frM)
04546 {
04547 #ifdef sdebug
04548 G4cout<<"G4Q::ModInMatCand:C="<<cPDG<<tmpT<<tmpTM<<"+"<<frM<<"="<<tmpTM+frM<<">tM="
04549 <<totMass<<G4endl;
04550 #endif
04551 curCand->SetPossibility(false);
04552 }
04553 G4QNucleus cNuc(cPDG);
04554 G4int cP = cNuc.GetZ();
04555 G4int cN = cNuc.GetN();
04556 G4int cL = cNuc.GetS();
04557
04558 #ifdef debug
04559 if(cPDG==90001000) G4cout<<"G4Q::MIM:->>cPDG=90001000<<-,possibility="<<poss<<G4endl;
04560 #endif
04561 if(eP>=cP&&eN>=cN&&eL>=cL&&poss)
04562 {
04563 G4double clME = 0.;
04564 G4double clMN = 0.;
04565 G4double renvM = 0;
04566 if(cP==eP&&cN==eN&&cL==eL)clME=cQPDG.GetMass();
04567 else
04568 {
04569 renvM = cQPDG.GetNuclMass(eP-cP,eN-cN,eL-cL);
04570 clME = envM-renvM;
04571 }
04572 if(cP==tP&&cN==tN&&cL==tL)clMN=cQPDG.GetMass();
04573 else
04574 {
04575 renvM = cQPDG.GetNuclMass(tP-cP,tN-cN,tL-cL);
04576 clMN = totM-renvM;
04577 }
04578 curCand->SetParPossibility(true);
04579 curCand->SetEBMass(clME);
04580 curCand->SetNBMass(clMN);
04581 #ifdef sdebug
04582 G4int envPDGC = theEnvironment.GetPDGCode();
04583 G4cout<<"G4Q:ModInMatCand:C="<<cPDG<<cNuc<<clME<<","<<clMN<<",E="<<envPDGC<<",M="
04584 <<renvM<<G4endl;
04585 #endif
04586 }
04587 else curCand->SetParPossibility(false);
04588 }
04589 }
04590 }
04591
04592
04593 void G4Quasmon::CalculateHadronizationProbabilities
04594 (G4double E, G4double kVal, G4LorentzVector k4M,G4bool piF, G4bool gaF, G4bool )
04595
04596 {
04597 static const G4double mPi0 = G4QPDGCode(111).GetMass();
04599 G4double kLS=E;
04600 kLS=k4M.e();
04601 G4int vap = nOfQ-3;
04602
04603
04604 G4double mQ2 = q4Mom.m2();
04605 G4double eQ = q4Mom.e();
04606 G4double mQ = sqrt(mQ2);
04607 G4double dk = kVal + kVal;
04608 G4double rQ2 = mQ2-dk*mQ;
04610 G4double mQk = mQ-dk;
04611 G4double var = theEnvironment.GetProbability();
04612 G4double vaf = 0;
04613 if(vap>0)vaf = var*mQk/kVal/vap;
04614
04615 G4double accumulatedProbability = 0.;
04616 G4double secondAccumProbability = 0.;
04617 G4int qBar =valQ.GetBaryonNumber();
04618 G4int nofU = valQ.GetU()- valQ.GetAU();
04619 G4int dofU = nofU+nofU;
04620 G4int nofD = valQ.GetD()- valQ.GetAD();
04621 G4int dofD = nofD+nofD;
04622 G4int qChg = valQ.GetCharge();
04623 G4int qIso = qBar-qChg-qChg;
04626 G4int maxC = theQCandidates.size();
04627 G4double totZ = theEnvironment.GetZ() + valQ.GetCharge();
04629 G4double envM = theEnvironment.GetMass();
04630 G4int envPDGC = theEnvironment.GetPDGCode();
04631 G4int envN = theEnvironment.GetN();
04632 G4int envZ = theEnvironment.GetZ();
04634 G4int envA = theEnvironment.GetA();
04635 G4QContent envQC=theEnvironment.GetQCZNS();
04637 #ifdef debug
04638 G4int absb = abs(qBar);
04639 G4int maxB = theEnvironment.GetMaxClust();
04640 G4cout<<"G4Q::CalcHadronizationProbab:Q="<<mQ<<valQ<<",v="<<vaf<<",r="<<var<<",mC="<<maxB
04641 <<",vap="<<vap<<",k="<<kVal<<G4endl;
04642 #endif
04643
04644 unsigned nHC=theQCandidates.size();
04645 #ifdef debug
04646 G4cout<<"G4Q::CHP: *** nHC="<<nHC<<G4endl;
04647 #endif
04648 if(nHC) for (unsigned index=0; index<nHC; index++)
04649 {
04650 G4QCandidate* curCand=theQCandidates[index];
04651 G4int cPDG = curCand->GetPDGCode();
04652 G4int aPDG = abs(cPDG);
04653 curCand->ClearParClustVector();
04654 G4double probability = 0.;
04655 G4double secondProbab = 0.;
04656 if ( (aPDG > 80000000 && envA > 0) || aPDG < 80000000)
04657 {
04658 G4int resPDG=0;
04659 G4double comb=0.;
04660 G4QContent candQC = curCand->GetQC();
04661 G4QContent tmpTQ=envQC+valQ-candQC;
04662 G4QNucleus tmpT(tmpTQ);
04663 G4double tmpTM=tmpT.GetMZNS();
04664 G4QPDGCode cQPDG(cPDG);
04665 G4double frM=cQPDG.GetMass();
04666 G4int cU=candQC.GetU()-candQC.GetAU();
04668 G4int cD=candQC.GetD()-candQC.GetAD();
04670 G4int dUD=abs(cU-cD);
04672 G4bool pos=curCand->GetPossibility()&&totMass>tmpTM+frM;
04673
04674 #ifdef pdebug
04675 G4bool pPrint= (abs(cPDG)%10 <3 && cPDG <80000000) || (cPDG >80000000 && frM <5000.);
04676
04677
04678
04679
04680
04681 if(pPrint) G4cout<<"G4Q::CHP:==****==>> c="<<cPDG<<",dUD="<<dUD<<",pos="<<pos<<",eA="
04682 <<envA<<",tM="<<totMass<<" > tmpTM+frM="<<tmpTM+frM<<G4endl;
04683 #endif
04684
04685 if(pos&&(cPDG<80000000||(cPDG>80000000&&cPDG!=90000000&&dUD<2)))
04686
04687
04688
04689 {
04690 G4QContent curQ = valQ;
04691 G4int baryn= candQC.GetBaryonNumber();
04692 G4int cC = candQC.GetCharge();
04693 G4double CB=0.;
04694 if(envA) CB=theEnvironment.CoulombBarrier(cC,baryn);
04696
04697 G4int cNQ= candQC.GetTot()-2;
04698
04699
04700 G4double resM=0.;
04701 #ifdef debug
04702 if(pPrint)G4cout<<"G4Q::CHP:B="<<baryn<<",C="<<cC<<",CB="<<CB<<",#q="<<cNQ<<G4endl;
04703 #endif
04704 if(cPDG>80000000&&cPDG!=90000000&&baryn<=envA)
04705 {
04706 G4int parentCounter=0;
04707 G4double pcomb=0.;
04708 G4double frM2=frM*frM;
04709 G4double qMax=frM+CB-kLS;
04710
04712 G4int iQmin=0;
04713 G4int iQmax=3;
04714 G4int oQmin=0;
04715 G4int oQmax=2;
04716 if (dofU<=nofD) iQmax=1;
04717 else if(dofD<=nofU) iQmin=1;
04718
04719 if(piF)
04720 {
04721 iQmin=0;
04722 iQmax=1;
04723 }
04724 #ifdef debug
04725 if(pPrint)G4cout<<"G4Q::CHP:***!!!***>>F="<<cPDG<<",mF="<<frM<<",iq:"<<iQmin<<","
04726 <<iQmax<<",kLS="<<kLS<<",kQCM="<<kVal<<",eA="<<envA<<G4endl;
04727 #endif
04728 if(iQmax>iQmin) for(int iq=iQmin; iq<iQmax; iq++)
04729 {
04730 G4double qFact=1.;
04731
04732 if(iq==1&&gaF)
04733 {
04734 qFact=4.;
04735 #ifdef debug
04736 if(pPrint) G4cout<<"G4Q::CHP:photon cap(gaF) is enhanced for Uquark"<<G4endl;
04737 #endif
04738 }
04739 G4double nqInQ=0.;
04740 if (!iq) nqInQ=valQ.GetD();
04741 else if(iq==1) nqInQ=valQ.GetU();
04742 else if(iq==2) nqInQ=valQ.GetS();
04743 comb=0.;
04744 #ifdef sdebug
04745 G4cout<<"G4Q::CHP:i="<<iq<<",cU="<<cU<<",cD="<<cD<<",omi="<<oQmin<<",oma="
04746 <<oQmax<<G4endl;
04747 #endif
04748 if(oQmax>oQmin) for(int oq=oQmin; oq<oQmax; oq++)
04749 {
04750 G4int shift= cQPDG.GetRelCrossIndex(iq, oq);
04751 G4QContent ioQC=cQPDG.GetExQContent(iq, oq);
04752 G4QContent resQC=valQ+ioQC;
04753 #ifdef sdebug
04754 G4cout<<"G4Q::CHP:iq="<<iq<<",oq="<<oq<<",QC="<<ioQC<<",rQC="<<resQC<<G4endl;
04755 #endif
04756 G4QPDGCode resQPDG(resQC);
04757 resPDG=resQPDG.GetPDGCode();
04758 G4int resQ=resQPDG.GetQCode();
04759 #ifdef pdebug
04760 if(pPrint) G4cout<<"G4Q::CHP:i="<<iq<<",o="<<oq<<ioQC<<",s="<<shift
04761 <<",cQPDG="<<cQPDG<<", residQC="<<resQC<<resQPDG<<G4endl;
04762 #endif
04763 G4int resD=resQC.GetD()-resQC.GetAD();
04764 G4int resU=resQC.GetU()-resQC.GetAU();
04765 G4int resS=resQC.GetS()-resQC.GetAS();
04766 G4int resA=resQC.GetBaryonNumber();
04767 G4bool rI=resA>0 && resU>=0 && resD>=0 &&
04768 (resU+resS>resD+resD||resD+resS>resU+resU);
04769
04770
04771
04772
04773
04774
04775 if (resQ > -2 && resPDG && resPDG != 10 && !rI &&
04776 (!piF || ( piF && (cPDG != 90001000 || G4UniformRand() < .333333) &&
04777 cPDG != 90002001 && cPDG != 90002002
04778 )
04779 )
04780 )
04781
04782
04783
04784
04785
04786
04787
04788 {
04789 G4int is=index+shift;
04790 if(shift!=7&&is<maxC)
04791 {
04792 G4QCandidate* parCand=theQCandidates[is];
04793 G4QContent parQC = parCand->GetQC();
04794 G4int barot = parQC.GetBaryonNumber();
04795 G4int charge= parQC.GetCharge();
04796 G4int possib=parCand->GetParPossibility();
04797 #ifdef pdebug
04798 if(pPrint) G4cout<<"G4Q::CHP:parentPossibility="<<possib<<",pZ="<<charge
04799 <<" <= envZ="<<envZ<<", pN="<<barot-charge<<" <= envN="
04800 <<envN<<", cPDG="<<cPDG<<G4endl;
04801 #endif
04802 if(possib && charge<=envZ && barot-charge<=envN)
04803 {
04804
04807 G4int isos = barot-charge-charge;
04808
04809
04810 G4double pUD= 1.;
04811 if(barot>2) pUD= pow(2.,abs(isos)-1);
04812
04813
04814 if(barot!=baryn) G4cerr<<"--Warning--G4Q::CHP:c="<<candQC<<",p="<<parQC
04815 <<",s="<<shift<<",i="<<index<<",s="<<is<<G4endl;
04816 G4int dI=qIso-isos;
04817 G4int dC=cC-charge;
04818 G4int dS=dI+dC;
04819 #ifdef pdebug
04820 if(pPrint)G4cout<<"G4Q::CHP: dS="<<dS<<", dI="<<dI<<", dC="<<dC<<", I="
04821 <<qIso<<",i="<<isos<<", C="<<cC<<",c="<<charge<<G4endl;
04822 #endif
04823
04824
04825
04826
04827
04828
04829
04830
04831
04832
04833
04834
04835
04836
04837
04838
04839
04840
04841
04842
04843
04844
04845
04846
04847
04848
04849
04850
04851
04852
04853
04854
04855
04856
04857
04858
04859
04860
04861
04862
04863
04864
04865
04866
04869 if ( (!piF && baryn < 5 ) ||
04870 ( piF && abs(dS) < 3) )
04871
04872
04873
04874
04875 {
04876 G4double pPP=parCand->GetPreProbability();
04877
04878 G4int parPDG=parCand->GetPDGCode();
04879 G4double boundM=parCand->GetEBMass();
04880 G4double nucBM =parCand->GetNBMass();
04881 #ifdef debug
04882 if(pPrint) G4cout<<"G4Q::CHP:c="<<cPDG<<",p="<<parPDG<<",bM="<<boundM
04883 <<",i="<<is<<",adr="<<parCand<<",pPP="<<pPP<<G4endl;
04884 #endif
04885
04886 G4double minM =0.;
04887 if (resPDG==10)minM=G4QChipolino(resQC).GetMass();
04888 else if(resPDG)minM=G4QPDGCode(resPDG).GetMass();
04889 G4double bNM2=nucBM*nucBM;
04890 G4double nDelta=0.;
04891 if(nucBM)nDelta=(frM2-bNM2)/(nucBM+nucBM);
04892 #ifdef pdebug
04893 G4int iniPDG =valQ.GetSPDGCode();
04894 G4double iniQM = G4QPDGCode(iniPDG).GetMass();
04895 G4double freeE = (mQ-iniQM)*iniQM;
04896 G4double kCut=boundM/2.+freeE/(iniQM+boundM);
04897 if(pPrint)G4cout<<"G4Q::CHP:r="<<resPDG<<",M="<<minM<<",k="<<kLS<<"<"
04898 <<kCut<<",E="<<E<<">"<<nDelta<<",p="<<parPDG<<G4endl;
04899 #endif
04900 if(resPDG && minM>0.)
04901 {
04902 #ifdef debug
04903 if(pPrint) G4cout<<"G4Q::CHP:fM="<<frM<<",bM="<<boundM<<",rM="
04904 <<tmpTM<<",tM="<<totMass<<G4endl;
04905 #endif
04906 G4double pmk=rMo*boundM/kLS;
04907
04908 G4double bM2=boundM*boundM;
04909 G4double eDelta=(frM2-bM2)/(boundM+boundM);
04910 G4double ked =kLS-eDelta;
04911 G4double dked=ked+ked;
04913 G4double kd =kLS-nDelta;
04914 G4double dkd=kd+kd;
04915
04916 G4double dkLS=kLS+kLS;
04917
04918
04919 G4double Em=(E-nDelta-CB)*(1.-frM/totMass);
04920
04921 G4double ne=1.-dked/(boundM+dkLS);
04922 G4double kf=1.;
04923 if(ne>0.&&ne<1.)kf=pow(ne,cNQ);
04924 #ifdef debug
04925 if(pPrint)G4cout<<"G4Q::CHP:<qi_DEF>="<<ne<<",k="<<kf<<",dk="<<dked
04926 <<",dkLS="<<dkLS<<",M="<<boundM<<",C="<<CB<<G4endl;
04927 #endif
04928
04929
04930 G4QContent rtQC=valQ+ioQC;
04931 if(envA-barot>bEn) rtQC+=bEnQC;
04932 else rtQC+=envQC-parQC;
04933 G4QNucleus rtN(rtQC);
04934 G4double rtM=rtN.GetGSMass();
04935 G4double rtEP=rEP;
04936
04937 if(envA-barot>bEn) rtEP+=mbEn;
04938 else rtEP+=envM-boundM;
04939 G4double rtE=rtEP-rMo;
04941 #ifdef debug
04942 G4QContent tmpEQ=envQC-parQC;
04943 if(pPrint) G4cout<<"G4Q::CHP:RN="<<tmpEQ<<"="<<envM-boundM<<"=eM="
04944 <<envM<<"-bM="<<boundM<<",E="<<rtE<<",eQC="<<envQC
04945 <<",pQC="<<parQC<<G4endl;
04946 #endif
04947 G4double mintM2=rtM*rtM+.1;
04948 G4double rtQ2=rtE*rtE-rMo*rMo;
04949
04950 G4double minBM=minM;
04951
04952 if ( (envA-barot <= bEn && envM > boundM) || envA-barot > bEn)
04953
04954 {
04955 minBM=rtM;
04956
04957 if(envA-barot > bEn) minBM-=mbEn;
04958 else minBM-=envM-boundM;
04959 }
04960 G4double minBM2=minBM*minBM+.1;
04961 G4double minM2=minM*minM+.1;
04962 #ifdef debug
04963 G4double ph=kf;
04964 if(pPrint) G4cout<<"G4Q::CHP:M2="<<minM2<<",R="<<rQ2<<",m="<<mintM2
04965 <<",RN2="<<rtQ2<<",q="<<(minM2-rQ2)/rEP/2<<",qN="
04966 <<(mintM2-rtQ2)/rtEP/2<<G4endl;
04967 #endif
04968 G4double newh=1.;
04969
04970
04971 G4double nc=1.-(dkLS-E-E+CB+CB)/boundM;
04972 G4double newl=0.;
04973 #ifdef debug
04974 if(pPrint) G4cout<<"G4Q::CHP:qi_k-E="<<nc<<",k="<<kLS<<",E="<<E
04975 <<",M="<<boundM<<G4endl;
04976 #endif
04977 if(nc > 0. && nc < 1. && nc < ne)
04978 {
04979 ne=nc;
04980 newh=pow(nc,cNQ);
04981 if(newh < kf) kf=newh;
04982 }
04983 else if(nc <= 0.) kf=0.;
04984
04985 G4double nk=1.-(dkd-Em-Em)/boundM;
04986 #ifdef debug
04987 if(pPrint) G4cout<<"G4Q::CHP:qi_R="<<nk<<",kd="<<kd<<",E="<<Em
04988 <<",M="<<totMass<<G4endl;
04989 #endif
04990 if(nk > 0. && nk < 1. && nk < ne)
04991 {
04992 ne=nk;
04993 newh=pow(nk,cNQ);
04994 if(newh<kf) kf=newh;
04995 }
04996 else if(nk <= 0.) kf=0.;
04997
04998
04999
05000
05001 #ifdef debug
05002
05003
05004 #endif
05005
05006
05007
05008
05009
05010
05011
05012
05013
05014
05015 G4double mix=boundM+E-CB;
05017 G4double st=0.;
05018 if(mix > frM) st=sqrt(mix*mix-frM2);
05019 G4double nq=1.-(dkLS-st-st)/boundM;
05020 #ifdef debug
05021 if(pPrint) G4cout<<"G4Q::CHP:qi_k-st="<<nq<<",st="<<st<<",m="
05022 <<mix<<",M="<<frM<<G4endl;
05023 #endif
05024 if(nq > 0. && nq < 1. && nq < ne)
05025
05026 {
05027 ne=nq;
05028 newh=pow(nq,cNQ);
05029 if(newh < kf) kf=newh;
05030 }
05031 else if(nq<=0.)kf=0.;
05032
05033 G4LorentzVector rq4M=q4Mom-k4M;
05034 G4ThreeVector k3V=k4M.vect().unit();
05035 G4ThreeVector rq3V=rq4M.vect().unit();
05036 G4bool atrest=(eQ-mQ)/mQ<.001||k3V.dot(rq3V)<-.999;
05037
05038
05039 if(mintM2>rtQ2)
05040
05041 {
05042 G4double nz=0.;
05043 if(atrest) nz=1.-(mintM2-rtQ2+pmk*dked)/(boundM*(rtEP+pmk));
05044 else nz=1.-(mintM2-rtQ2)/(boundM*rtEP);
05045
05046
05047 #ifdef debug
05048 if(pPrint) G4cout<<"G4Q::CHP:q="<<nz<<",a="<<atrest<<",M2="
05049 <<mintM2<<">"<<rtQ2<<G4endl;
05050 #endif
05051 if(nz > 0. && nz < 1. && nz < ne)
05052 {
05053 ne=nz;
05054 newh=pow(nz,cNQ);
05055 if(newh < kf) kf=newh;
05056 }
05057 else if(nz <= 0.) kf=0.;
05058 }
05059
05060
05061
05062
05063
05064
05065
05066
05067
05068
05069
05070
05071
05072
05073
05074
05075
05076
05077
05078
05079
05080 if (minBM2 > rQ2 &&
05081 (!piF ||
05082 (piF && cPDG!=90000001 && cPDG!=90001001 && cPDG!=90001002)
05083 )
05084 )
05085
05086
05087 {
05088 G4double nz=0.;
05089 if(atrest) nz=1.-(minBM2-rQ2+pmk*dked)/(boundM*(rEP+pmk));
05090 else nz=1.-(minBM2-rQ2)/(boundM*rEP);
05091
05092
05093 #ifdef debug
05094 if(pPrint) G4cout<<"G4Q::CHP:q="<<nz<<",a="<<atrest<<",QM2="
05095 <<minM2<<">"<<rQ2<<G4endl;
05096 #endif
05097 if(nz>0.&&nz<1.&&nz<ne)
05098 {
05099 ne=nz;
05100 newh=pow(nz,cNQ);
05101 if(newh<kf) kf=newh;
05102 }
05103 else if(nz<=0.)kf=0.;
05104 }
05105
05106
05107
05108
05109
05110
05111
05112
05113
05114
05115
05116
05117
05118
05119
05120
05121 if ( minM2 > rQ2 &&
05122 ( (!piF && baryn > 4) ||
05123
05124 (piF && cPDG != 90000001
05125 && cPDG != 90001001)
05126 )
05127 )
05128
05129
05130
05131 {
05132 G4double nz=0.;
05133 if(atrest) nz=1.-(minM2-rQ2+pmk*dked)/(boundM*(rEP+pmk));
05134 else nz=1.-(minM2-rQ2)/(boundM*rEP);
05135
05136
05137 #ifdef debug
05138 if(pPrint) G4cout<<"G4Q::CHP:q="<<nz<<",a="<<atrest<<",QM2="
05139 <<minM2<<">"<<rQ2<<G4endl;
05140 #endif
05141 if(nz>0.&&nz<1.&&nz<ne)
05142 {
05143 ne=nz;
05144 newh=pow(nz,cNQ);
05145 if(newh<kf) kf=newh;
05146 }
05147 else if(nz<=0.)kf=0.;
05148 }
05149 if(kf<0.)kf=0.;
05150 if(kf>1.)kf=1.;
05151 G4double high = kf;
05152 #ifdef debug
05153 if(pPrint) G4cout<<"G4Q::CHP:"<<kf<<",minM2="<<minM2<<",rQ2="<<rQ2
05154 <<G4endl;
05155 #endif
05156 G4double lz=1.-dked/boundM;
05157
05158
05159
05160
05161 G4double low=0.;
05162 if(lz>0.&&lz<1.)low=pow(lz,cNQ);
05163 else if(lz>=1.)low=1.;
05164 #ifdef debug
05165 G4double pl=low;
05166 if(pPrint) G4cout<<"G4Q::CHP:<qa_DEF>="<<lz<<", eDel="<<eDelta
05167 <<",nDel="<<nDelta<<G4endl;
05168 #endif
05169
05170
05171 G4double tms=kLS+eDelta+Em;
05172 G4double le=1.-(tms+tms)/boundM;
05173 #ifdef debug
05174 if(pPrint) G4cout<<"G4Q::CHP:qa_t="<<le<<",k="<<kLS<<",E="<<Em
05175 <<",bM="<<boundM<<G4endl;
05176 #endif
05177 if(le>0.&&le<1.&&le>lz)
05178 {
05179 lz=le;
05180 newl=pow(le,cNQ);
05181 if(newl>low) low=newl;
05182 }
05183 else if(le>=1.)low=1.;
05184
05185
05186
05187 G4double lk=1.-(dkLS+E+E-CB-CB)/boundM;
05188 #ifdef debug
05189 if(pPrint) G4cout<<"G4Q::CHP:qa_k+E="<<lk<<",k="<<kLS<<",E="<<E
05190 <<",M="<<boundM<<G4endl;
05191 #endif
05192 if(lk>0.&&lk<1.&&lk>lz)
05193 {
05194 lz=lk;
05195 newl=pow(lk,cNQ);
05196 if(newl>low) low=newl;
05197 }
05198 else if(lk>=1.)low=1.;
05199
05200
05201
05202
05203 #ifdef debug
05204
05205
05206 #endif
05207
05208
05209
05210
05211
05212
05213
05214
05215
05216 G4double lp=1.-(dkLS+sr+sr)/boundM;
05217 #ifdef debug
05218 if(pPrint) G4cout<<"G4Q::CHP:qa_k+sr="<<lp<<",sr="<<sr
05219 <<",M="<<frM<<G4endl;
05220 #endif
05221 if(lp>0.&&lp<1.&&lp>lz)
05222 {
05223 lz=lp;
05224 newl=pow(lp,cNQ);
05225 if(newl>low) low=newl;
05226 }
05227 else if(lp>=1.)low=1.;
05228
05229
05230 if(totZ > cC)
05231
05232 {
05233 G4double qmaCB=boundM-qMax;
05234
05235 G4double nz=1.-(qmaCB+qmaCB)/boundM;
05236 #ifdef debug
05237 if(pPrint) G4cout<<"G4Q::CHP:<qa_CB>="<<nz<<",m="<<qmaCB<<",CB="
05238 <<CB<<G4endl;
05239 #endif
05240 if(nz>0.&&nz>lz)
05241 {
05242 newl=pow(nz,cNQ);
05243 if(newl>low) low=newl;
05244 }
05245 else if(nz>1.) low=10.;
05246 }
05247
05248 kf-=low;
05249 #ifdef debug
05250 if(pPrint) G4cout<<"G4Q::CHP:>>"<<cPDG<<",l="<<low<<",h="<<high
05251 <<",ol="<<pl<<",oh="<<ph<<",nl="<<newl<<",nh="
05252 <<newh<<",kf="<<kf<<",d="<<eDelta<<G4endl;
05253 #endif
05254 G4double probab=0.;
05255 if(kf>0)
05256 {
05257 kf*=boundM/kLS/cNQ;
05258 G4int noc=cQPDG.GetNumOfComb(iq, oq);
05259
05261 probab=qFact*kf*nqInQ*pPP*noc/pUD;
05262
05263
05264
05265
05266
05267
05268
05269 G4QContent rQQC = valQ+ioQC;
05270 G4int BarRQC=rQQC.GetBaryonNumber();
05271 G4int StrRQC=rQQC.GetStrangeness();
05272 if(BarRQC==2 && !StrRQC)
05273 {
05274 G4int ChgRQC=rQQC.GetCharge();
05275 if(ChgRQC==1) probab/=2;
05276 else probab*=2;
05277 }
05278 #ifdef debug
05279 if(pPrint)G4cout<<"G4Q::CHP:prob="<<probab<<",qF="<<qFact<<",iq="
05280 <<iq<<",oq="<<oq<<",Pho4M="<<phot4M<<",pUD="<<pUD
05281 <<",pPP="<<pPP<<G4endl;
05282 #endif
05283 if(probab<0.) probab=0.;
05284 }
05285 pcomb += probab;
05286 G4QParentCluster* curParC = new G4QParentCluster(parPDG,pcomb);
05287 curParC->SetTransQC(ioQC);
05288 curParC->SetLow(low);
05289 curParC->SetHigh(high);
05290 curParC->SetEBMass(boundM);
05291 curParC->SetNBMass(nucBM);
05292 curParC->SetEBind(eDelta);
05293 curParC->SetNBind(nDelta);
05294 curParC->SetNQPart2(cNQ);
05295 #ifdef sdebug
05296 G4cout<<"G4Q::CalcHP: FillParentClaster="<<*curParC<<G4endl;
05297 #endif
05298 curCand->FillPClustVec(curParC);
05299 comb += probab;
05300 #ifdef pdebug
05301 if(pPrint) G4cout<<"G4Q::CHP:in="<<index<<",cPDG="<<cPDG<<",pc"<<parentCounter
05302 <<parQC<<",Env="<<theEnvironment<<",comb="<<comb
05303 <<",posib="<<parCand->GetParPossibility()<<G4endl;
05304 #endif
05305 parentCounter++;
05306 }
05307 }
05308 #ifdef sdebug
05309 else G4cout<<"***G4Q::CHP:dI="<<dI<<",cC="<<cC<<G4endl;
05310 #endif
05311 }
05312 }
05313 }
05314 probability+=comb;
05315 #ifdef pdebug
05316 if(pPrint) G4cout<<"G4Q::CHPr: probab="<<probability<<"("<<comb<<"),iq="<<iq
05317 <<",oq="<<oq<<G4endl;
05318 #endif
05319 }
05320 }
05321 }
05322 else if(cPDG<80000000)
05323 {
05324
05325 G4int curnh=theQHadrons.size();
05326 G4int npip=0;
05327 G4int npin=0;
05328 G4int npiz=0;
05329 for (G4int ind=0; ind<curnh; ind++)
05330 {
05331 G4int curhPDG=theQHadrons[ind]->GetPDGCode();
05332 if (curhPDG== 111) npiz++;
05333 if (curhPDG== 211) npip++;
05334 if (curhPDG==-211) npin++;
05335 }
05336
05337 comb = valQ.NOfCombinations(candQC);
05338 if(!comb)
05339 {
05340 if ( (aPDG==111)|(aPDG==211) ) comb=1.;
05341 else if ( (aPDG==311)|(aPDG==321) ) comb=SSin2Gluons;
05342 }
05343 if(cPDG== 211&&npip>0) comb*=(npip+1);
05344 if(cPDG==-211&&npip>0) comb*=(npin+1);
05345 if(cPDG==111||cPDG==221||cPDG==331||cPDG==113||cPDG==223||cPDG==333||cPDG==115||
05346 cPDG==225||cPDG==335||cPDG==117||cPDG==227||cPDG==337||cPDG==110||cPDG==220||
05347 cPDG==330)
05348 {
05349 G4QContent tQCd(1,0,0,1,0,0);
05350 G4QContent tQCu(0,1,0,0,1,0);
05351 G4QContent tQCs(0,0,1,0,0,1);
05352 G4double cmd=valQ.NOfCombinations(tQCd);
05353 G4double cmu=valQ.NOfCombinations(tQCu);
05354 G4double cms=valQ.NOfCombinations(tQCs);
05355 if(cPDG!=333&&cPDG!=335&&cPDG!=337) comb=(cmd+cmu)/2.;
05356
05357 if(cPDG==331||cPDG==221) comb =(comb + cms)/4.;
05358 if(cPDG==113) comb*=4.;
05359 if(cPDG==223) comb*=2.;
05360 if(cPDG==111&&npiz>0) comb*=(npiz+1);
05361 #ifdef debug
05362 if(abs(cPDG)<3) G4cout<<"G4Q::CHP:comb="<<comb<<",cmd="<<cmd<<",cmuu="<<cmu
05363 <<",cms="<<cms<<G4endl;
05364 #endif
05365 }
05366 curQ -= candQC;
05367 resPDG = curQ.GetSPDGCode();
05368 G4QContent resTQC = curQ+envQC;
05369 G4double resTM=G4QPDGCode(resTQC.GetSPDGCode()).GetMass();
05370 #ifdef debug
05371 G4bool priCon = aPDG < 10000 && aPDG%10 < 3;
05372 if(priCon) G4cout<<"G4Q::CHP:***>>cPDG="<<cPDG<<",cQC="<<candQC<<",comb="<<comb
05373 <<",curQC="<<curQ<<",mQ="<<mQ<<",ab="<<absb<<G4endl;
05374 #endif
05375 if(resPDG==221 || resPDG==331)
05376 {
05377 resPDG=111;
05378 resTM=mPi0;
05379 }
05380 #ifdef debug
05381 if(priCon) G4cout<<"G4Q::CHP:cPDG="<<cPDG<<",c="<<comb<<",rPDG/QC="<<resPDG<<curQ
05382 <<",tM="<<totMass<<">"<<frM-CB+resTM<<"=fM="<<frM<<"+rM="<<resTM
05383 <<"-CB="<<CB<<G4endl;
05384 #endif
05385 if (comb && resPDG && totMass > frM-CB+resTM &&
05386 ((resPDG > 80000000 && resPDG != 90000000) || resPDG<10000) )
05387 {
05388 #ifdef debug
05389 if(priCon) G4cout<<"G4Q::CHP:ind="<<index<<",qQC="<<valQ<<mQ<<",cPDG="<<cPDG
05390 <<",rPDG="<<resPDG<<curQ<<G4endl;
05391 #endif
05392 if(resPDG!=10)resM=G4QPDGCode(resPDG).GetMass();
05393 else resM=G4QChipolino(curQ).GetMass();
05394 G4int resQCode=G4QPDGCode(curQ).GetQCode();
05395 #ifdef debug
05396 if(priCon) G4cout<<"G4Q::CHP:rM/QC="<<resM<<curQ<<",E="<<envPDGC<<",rQC="
05397 <<resQCode<<G4endl;
05398 #endif
05399
05400 if(envPDGC>80000000 && envPDGC!=90000000 && resM>0. &&
05401 resPDG!=10 && resPDG!=1114 && resPDG!=2224)
05402 {
05403 G4QContent rtQC=curQ+envQC;
05404 G4QNucleus rtN(rtQC);
05405 G4double rtM =rtN.GetMZNS();
05406 G4double bnRQ=rtM-envM;
05407 #ifdef debug
05408 if(priCon) G4cout<<"G4Q::CHP: **Rec**,RQMass="<<bnRQ<<",envM="<<envM<<",rtM="
05409 <<rtM<<G4endl;
05410 #endif
05411
05412 if(bnRQ<resM) resM=bnRQ;
05413 }
05414 #ifdef debug
05415 if(aPDG<10000&&aPDG%10<3)
05416
05417 G4cout<<"G4Q::CHP: resM="<<resM<<", resQCode="<<resQCode<<G4endl;
05418 #endif
05419 if(resM>0. && resQCode>-2)
05420 {
05421 G4double limM=mQ-resM;
05422 G4double rndM=GetRandomMass(cPDG,limM);
05423 #ifdef debug
05424 G4double cMass=G4QPDGCode(cPDG).GetMass();
05425 if(aPDG<10000&&aPDG%10<3)
05426
05427 G4cout<<"G4Q::CHP:rndM="<<rndM<<",limM="<<limM<<" > cM="<<cMass<<" ,rM+fM="
05428 <<resM+rndM<<" < mQ="<<mQ<<G4endl;
05429 #endif
05430
05431 if(rndM>0. && resM+rndM<mQ)
05432 {
05433 curCand->SetEBMass(rndM);
05434 curCand->SetNBMass(rndM);
05435 G4double mH2 = rndM*rndM;
05436 G4double rHk = mH2/dk;
05437 G4double zMax = 1.-rHk/mQ;
05438 G4double mR2 = resM*resM;
05439 G4double zMin=0.;
05440
05441
05442 if(qBar) zMin= mR2/mQ/(mQ-dk);
05443 G4double possibility=zMax-zMin;
05444 #ifdef debug
05445 if(priCon) G4cout<<"G4Q::CHP:M="<<rndM<<",ps="<<possibility<<",zMax="<<zMax
05446 <<",rHk="<<rHk<<",mQ="<<mQ<<",dk="<<dk<<",zMin="<<zMin
05447 <<",mR2="<<mR2<<",rM="<<resM<<"; "<<mQ*(mQ-dk)<<G4endl;
05448 #endif
05449 if (resPDG==10)
05450 {
05451 G4double rM2 = mQk*(mQ-rHk);
05452 if(rM2<resM*resM) possibility = 0.;
05453 }
05454 if (possibility>0. && vap>0 && zMax>zMin)
05455 {
05456 probability = vaf*(pow(zMax, vap)-pow(zMin, vap));
05457 #ifdef debug
05458 if(priCon) G4cout<<"G4Q::CHP:#"<<index<<",mH2="<<mH2<<",nQ="<<nOfQ<<",p="
05459 <<probability<<",vf="<<vaf<<",vp="<<vap<<",zMax="<<zMax
05460 <<",zMin="<<zMin<<G4endl;
05461 #endif
05462
05463
05464
05465
05466
05467
05468
05469
05470
05471
05472 if(cPDG==110||cPDG==220||cPDG==330) probability*=comb;
05473 else probability*=comb*(abs(cPDG)%10);
05474 G4int BarRQC=curQ.GetBaryonNumber();
05475 G4int StrRQC=curQ.GetStrangeness();
05476 if(BarRQC==2 && !StrRQC)
05477 {
05478 G4int ChgRQC=curQ.GetCharge();
05479 if(ChgRQC==1) probability/=2;
05480 else probability*=2;
05481 }
05482
05483 }
05484 }
05485 else
05486 {
05487 #ifdef debug
05488 if(priCon) G4cout<<"G4Q::CHP:cM=0[cPDG"<<cPDG<<"],mQ/QC="<<mQ<<valQ<<",rM="
05489 <<resM<<curQ<<G4endl;
05490 #endif
05491 }
05492 }
05493 else
05494 {
05495 #ifdef debug
05496 if(priCon) G4cout<<"***G4Q::CHP: M=0, #"<<index<<valQ<<",cPDH="<<cPDG<<"+rP="
05497 <<resPDG<<curQ<<G4endl;
05498 #endif
05499 }
05500 }
05501 else
05502 {
05503 probability=0.;
05504 #ifdef debug
05505 if(priCon) G4cout<<"G4Q::CHP:"<<index<<valQ<<",PDG="<<cPDG<<"+r="<<resPDG<<curQ
05506 <<":c=0(!) || tM="<<totMass<<"<"<<frM-CB+resTM<<" = fM="<<frM
05507 <<"+rTM="<<resTM<<"-CB="<<CB<< G4endl;
05508 #endif
05509 }
05510
05511 }
05512 else probability=0.;
05513 #ifdef debug
05514 G4int aPDG = abs(cPDG);
05515 if((cPDG>90000000 && baryn<5) || (aPDG<10000 && aPDG%10<3))
05516 G4cout<<"G4Q::CHP:^^cPDG="<<cPDG<<",p="<<pos<<",rPDG="<<resPDG<<curQ<<resM<<",p="
05517 <<probability<<",a="<<accumulatedProbability<<",sp="<<secondProbab<<G4endl;
05518 #endif
05519 }
05520 }
05521 curCand->SetRelProbability(probability);
05522 accumulatedProbability += probability;
05523 curCand->SetIntegProbability(accumulatedProbability);
05524 curCand->SetSecondRelProb(secondProbab);
05525 secondAccumProbability += secondProbab;
05526 curCand->SetSecondIntProb(secondAccumProbability);
05527 }
05528 }
05529
05530
05531 G4bool G4Quasmon::CheckGroundState(G4bool corFlag)
05532 {
05533 static const G4QContent neutQC(2,1,0,0,0,0);
05534 static const G4QContent protQC(1,2,0,0,0,0);
05535 static const G4QContent lambQC(1,1,1,0,0,0);
05536 static const G4QContent deutQC(3,3,0,0,0,0);
05537 static const G4QContent alphQC(6,6,0,0,0,0);
05538 static const G4double mNeut= G4QPDGCode(2112).GetMass();
05539 static const G4double mProt= G4QPDGCode(2212).GetMass();
05540 static const G4double mLamb= G4QPDGCode(3122).GetMass();
05541 static const G4double mDeut= G4QPDGCode(2112).GetNuclMass(1,1,0);
05542 static const G4double mAlph= G4QPDGCode(2112).GetNuclMass(2,2,0);
05543 static const G4QNucleus vacuum= G4QNucleus(G4LorentzVector(0.,0.,0.,0.),90000000);
05547 G4int resQPDG=valQ.GetSPDGCode();
05548 G4double resQMa=G4QPDGCode(resQPDG).GetMass();
05549 G4double resEMa=0.;
05550 G4int bsCond=0;
05551 G4LorentzVector enva4M=G4LorentzVector(0.,0.,0.,0.);
05552 G4LorentzVector reTLV=q4Mom;
05553 G4double resSMa=resQMa;
05554 #ifdef debug
05555 G4cout<<"G4Q::CheckGS: EnvPDG="<<theEnvironment.GetPDG()<<",Quasmon="<<resQPDG<<G4endl;
05556 #endif
05557 if(theEnvironment.GetPDG()!=90000000)
05558 {
05559 resEMa=theEnvironment.GetMZNS();
05560 #ifdef debug
05561 G4cout<<"G4Q::CheckGS: Environment Mass="<<resEMa<<G4endl;
05562 #endif
05563 enva4M=G4LorentzVector(0.,0.,0.,resEMa);
05564 reTLV+=enva4M;
05565 resSMa+=resEMa;
05566 }
05567 else
05568 {
05569 G4QNucleus tmpQN(valQ,reTLV);
05570
05571 bsCond = tmpQN.SplitBaryon();
05572 #ifdef debug
05573 G4cout<<"G4Q::CheckGS: No environment, theOnlyQ="<<tmpQN<<",bsCond="<<bsCond<<G4endl;
05574 #endif
05575 if(bsCond)
05576 {
05577 G4QContent fragmQC=protQC;
05578 G4double fragmM=mProt;
05579 if(bsCond==2112)
05580 {
05581 fragmQC=neutQC;
05582 fragmM=mNeut;
05583 }
05584 else if(bsCond==3122)
05585 {
05586 fragmQC=lambQC;
05587 fragmM=mLamb;
05588 }
05589 else if(bsCond==90001001)
05590 {
05591 fragmQC=deutQC;
05592 fragmM=mDeut;
05593 }
05594 else if(bsCond==90002002)
05595 {
05596 fragmQC=alphQC;
05597 fragmM=mAlph;
05598 }
05599 G4QContent rsQC=valQ-fragmQC;
05600 G4QNucleus rsQN(rsQC);
05601 G4double rsMass=rsQN.GetGSMass();
05602 G4LorentzVector quas4M = G4LorentzVector(0.,0.,0.,rsMass);
05603 G4QHadron* quasH = new G4QHadron(rsQC, quas4M);
05604 G4LorentzVector frag4M = G4LorentzVector(0.,0.,0.,fragmM);
05605 G4QHadron* fragH = new G4QHadron(fragmQC, frag4M);
05606 if(G4QHadron(reTLV).DecayIn2(frag4M,quas4M))
05607 {
05608 quasH->Set4Momentum(quas4M);
05609 FillHadronVector(quasH);
05610 fragH->Set4Momentum(frag4M);
05611 FillHadronVector(fragH);
05612 return true;
05613 }
05614 else
05615 {
05616 delete quasH;
05617 delete fragH;
05618 }
05619 }
05620 }
05621 G4QContent envaQC = theEnvironment.GetQCZNS();
05622 G4double resTMa=reTLV.m();
05623
05624 G4int nOfOUT = theQHadrons.size();
05625 #ifdef debug
05626 G4cout<<"G4Q::CheckGS: (totM="<<resTMa<<" < rQM+rEM="<<resSMa<<" || rEM="<<resEMa
05627 <<"=0 && "<<bsCond<<"=0) && n="<<nOfOUT<<" >0"<<G4endl;
05628 #endif
05629 if ( (resTMa < resSMa || (!resEMa && !bsCond) ) && nOfOUT > 0 && corFlag)
05630 {
05631
05632 G4QHadron* theLast = theQHadrons[nOfOUT-1];
05633 if(!(theLast->GetNFragments()) && theLast->GetPDGCode()!=22)
05634 {
05635 G4LorentzVector hadr4M=theLast->Get4Momentum();
05636 G4double hadrMa=hadr4M.m();
05637 G4LorentzVector tmpTLV=reTLV+hadr4M;
05638 #ifdef debug
05639 G4cout<<"G4Q::CheckGS:YES,T="<<tmpTLV<<tmpTLV.m()<<">rM+hM="<<resSMa+hadrMa<<G4endl;
05640 #endif
05641 if(tmpTLV.m()>resSMa+hadrMa)
05642 {
05643 if(resEMa)
05644 {
05645 G4LorentzVector quas4M = G4LorentzVector(0.,0.,0.,resQMa);
05646 if(!G4QHadron(tmpTLV).DecayIn3(hadr4M,quas4M,enva4M))
05647 {
05648 G4cerr<<"---Warning---G4Q::CheckGS:DecIn Fragm+ResQ+ResEnv Error"<<G4endl;
05649 return false;
05650 }
05651 else
05652 {
05653
05654
05655 G4QHadron* envH = new G4QHadron(envaQC,enva4M);
05656 FillHadronVector(envH);
05657 theEnvironment = vacuum;
05658 G4QHadron* quasH = new G4QHadron(valQ, quas4M);
05659
05660 FillHadronVector(quasH);
05661 theLast->Set4Momentum(hadr4M);
05662 }
05663 }
05664 else
05665 {
05666 G4LorentzVector quas4M = G4LorentzVector(0.,0.,0.,resQMa);
05667 G4QHadron* quasH = new G4QHadron(valQ, quas4M);
05668 if(!G4QHadron(tmpTLV).DecayIn2(hadr4M,quas4M))
05669 {
05670 delete quasH;
05671 G4cerr<<"---Warning---G4Q::CheckGS: Decay in Fragm+ResQ Error"<<G4endl;
05672 return false;
05673 }
05674 else
05675 {
05676
05677 theLast->Set4Momentum(hadr4M);
05678 quasH->Set4Momentum(quas4M);
05679 FillHadronVector(quasH);
05680 }
05681 }
05682 }
05683 else
05684 {
05685 if(nOfOUT>1 && corFlag)
05686 {
05687 G4QHadron* thePrev = theQHadrons[nOfOUT-2];
05688 if(thePrev->GetNFragments()||thePrev->GetPDGCode()==22)return false;
05689 G4LorentzVector prev4M=thePrev->Get4Momentum();
05690 G4double prevMa=prev4M.m();
05691 tmpTLV+=prev4M;
05692 G4QContent totQC=valQ+envaQC;
05693 G4int totPDG=totQC.GetSPDGCode();
05694 G4double totQMa=G4QPDGCode(totPDG).GetMass();
05695 G4double totNMa=tmpTLV.m();
05696 #ifdef debug
05697 G4cout<<"G4Q::CheckGS:NO, M="<<tmpTLV<<totNMa<<">"<<totQMa+hadrMa+prevMa<<G4endl;
05698 #endif
05699 if(totNMa>totQMa+hadrMa+prevMa)
05700 {
05701 G4LorentzVector nuc4M = G4LorentzVector(0.,0.,0.,totQMa);
05702 if(!G4QHadron(tmpTLV).DecayIn3(hadr4M,prev4M,nuc4M))
05703 {
05704 G4cerr<<"---Warning---G4Q::CheckGS:DecIn3 ResN+Last+Prev Error"<<G4endl;
05705 return false;
05706 }
05707 else
05708 {
05709
05710 G4QHadron* envH = new G4QHadron(totQC,nuc4M);
05711 FillHadronVector(envH);
05712 theEnvironment = vacuum;
05713 theLast->Set4Momentum(hadr4M);
05714 thePrev->Set4Momentum(prev4M);
05715 #ifdef debug
05716 G4cout<<"G4Q::CheckGS: Yes, Check D4M="<<tmpTLV-hadr4M-prev4M-nuc4M<<G4endl;
05717 #endif
05718 }
05719 }
05720 else
05721 {
05722 G4QContent tmpLNQ=totQC+thePrev->GetQC();
05723 G4int resLPDG =tmpLNQ.GetSPDGCode();
05724 G4double resLastM=G4QPDGCode(resLPDG).GetMass();
05725 G4QContent tmpPNQ=totQC+theLast->GetQC();
05726 G4int resPPDG =tmpPNQ.GetSPDGCode();
05727 G4double resPrevM=G4QPDGCode(resPPDG).GetMass();
05729 #ifdef debug
05730 G4cout<<"G4Quasm::CheckGS: NO, tM="<<totNMa<<" > rp+l="<<resLastM+hadrMa
05731 <<" || > rl+p="<<resPrevM+prevMa<<G4endl;
05732 #endif
05733 if (totNMa>resLastM+hadrMa)
05734 {
05735 theQHadrons.pop_back();
05736 theQHadrons.pop_back();
05737 theQHadrons.push_back(theLast);
05738 delete thePrev;
05739 thePrev=theLast;
05740 resPPDG=resLPDG;
05741 resPrevM=resLastM;
05742 prev4M = hadr4M;
05743 }
05744 else if (totNMa>resPrevM+prevMa)
05745 {
05746 theQHadrons.pop_back();
05747 delete theLast;
05748 }
05749 else return false;
05750 G4LorentzVector nuc4M = G4LorentzVector(0.,0.,0.,resPrevM);
05751 if(!G4QHadron(tmpTLV).DecayIn2(prev4M,nuc4M))
05752 {
05753 G4cerr<<"---Warning---G4Q::CheckGS:DecIn2 (ResN+Last)+Prev Error"<<G4endl;
05754 return false;
05755 }
05756 else
05757 {
05758
05759 G4QHadron* envH = new G4QHadron(resPPDG,nuc4M);
05760 FillHadronVector(envH);
05761 theEnvironment = vacuum;
05762 thePrev->Set4Momentum(prev4M);
05763 #ifdef debug
05764 G4cout<<"G4Q::CheckGS:Yes, Check D4M="<<tmpTLV-prev4M-nuc4M<<G4endl;
05765 #endif
05766 }
05767 }
05768 }
05769 else return false;
05770 }
05771 }
05772 else return false;
05773 }
05774 else return false;
05775 return true;
05776 }
05777
05778
05779 G4QHadronVector* G4Quasmon::DecayQHadron(G4QHadron* qH)
05780 {
05781 G4QHadronVector* theFragments = new G4QHadronVector;
05782 G4QPDGCode theQPDG = qH->GetQPDG();
05783 G4int thePDG = theQPDG.GetPDGCode();
05784 G4int pap = 0;
05785 if(thePDG<0) pap = 1;
05786 G4LorentzVector t = qH->Get4Momentum();
05787 G4double m_value = t.m();
05788
05789 G4QDecayChanVector decV = theWorld->GetQParticle(theQPDG)->GetDecayVector();
05790 G4int nChan = decV.size();
05791 #ifdef debug
05792 G4cout<<"G4Quasm::DecQHadron: PDG="<<thePDG<<",m="<<m_value<<",("<<nChan<<" channels)"<<G4endl;
05793 #endif
05794 if(nChan)
05795 {
05796 G4int i=0;
05797 if(nChan>1)
05798 {
05799 G4double rnd = G4UniformRand();
05800 for(i=0; i<nChan; i++)
05801 {
05802 G4QDecayChan* dC = decV[i];
05803 #ifdef debug
05804 G4cout<<"G4Quasmon::DecaQHadr:i="<<i<<",r="<<rnd<<"<dl="<<dC->GetDecayChanLimit()
05805 <<", mm="<<dC->GetMinMass()<<G4endl;
05806 #endif
05807 if(rnd<dC->GetDecayChanLimit() && m_value>dC->GetMinMass()) break;
05808 }
05809 if(i>nChan-1) i=nChan-1;
05810 }
05811 G4QPDGCodeVector cV=decV[i]->GetVecOfSecHadrons();
05812 G4int nPart=cV.size();
05813 #ifdef debug
05814 G4cout<<"G4Quasmon::DecayQHadron: resi="<<i<<",nP="<<nPart<<":"<<cV[0]->GetPDGCode()
05815 <<","<<cV[1]->GetPDGCode();
05816 if(nPart>2) G4cout<<","<<cV[2]->GetPDGCode();
05817 G4cout<<G4endl;
05818 #endif
05819 if(nPart<2||nPart>3)
05820 {
05821 G4cerr<<"---Warning---G4Q::DecayQHadr:n="<<nPart<<",ch#"<<i<<",PDG="<<thePDG<<G4endl;
05822 theFragments->push_back(qH);
05823 return theFragments;
05824 }
05825 #ifdef debug
05826 G4cout<<"G4Q::DecQH:Decay("<<ElMaDecays<<") PDG="<<thePDG<<t<<m_value<<",nP="<<nPart<<G4endl;
05827 #endif
05828 if(nPart==2)
05829 {
05830 G4QHadron* fHadr;
05831 G4QHadron* sHadr;
05832 G4int fPDG=cV[0]->GetPDGCode();
05833 G4int sPDG=cV[1]->GetPDGCode();
05834
05835 if ( (fPDG != 22 && sPDG != 22) || ElMaDecays) {
05836 #ifdef debug
05837 G4cout<<"G4Q::DecQH:Yes2,fPDG="<<fPDG<<",sPDG="<<sPDG<<",EMF="<<ElMaDecays<<G4endl;
05838 #endif
05839 if(cV[0]->GetWidth()==0.)
05840 {
05841 fHadr = new G4QHadron(cV[0]->GetPDGCode());
05842 if(cV[1]->GetWidth()==0.)sHadr = new G4QHadron(sPDG);
05843 else
05844 {
05845 G4QParticle* sPart=theWorld->GetQParticle(cV[1]);
05846 G4double sdm = m_value - fHadr->GetMass();
05847 sHadr = new G4QHadron(sPart,sdm);
05848 if(sPDG<0) sHadr->MakeAntiHadron();
05849 }
05850 }
05851 else
05852 {
05853 G4QParticle* sPart=theWorld->GetQParticle(cV[1]);
05854 G4double mim = sPart->MinMassOfFragm();
05855 G4double fdm = m_value - mim;
05856 G4QParticle* fPart=theWorld->GetQParticle(cV[0]);
05857 fHadr = new G4QHadron(fPart,fdm);
05858 if(fPDG<0) fHadr->MakeAntiHadron();
05859 G4double fm=fHadr->GetMass();
05860 G4double sdm = m_value - fm;
05861 sHadr = new G4QHadron(sPart,sdm);
05862 if(sPDG<0) sHadr->MakeAntiHadron();
05863 #ifdef debug
05864 G4cout<<"G4Q::DQH:M="<<m_value<<",mi="<<mim<<",fd="<<fdm<<",fm="<<fm<<",sd="<<sdm
05865 <<",sm="<<sHadr->GetMass()<<G4endl;
05866 #endif
05867 }
05868 #ifdef debug
05869 G4cout<<"G4Q::DQH:(DecayIn2)1="<<fHadr->GetMass()<<",2="<<sHadr->GetMass()<<G4endl;
05870 #endif
05871 if(pap)
05872 {
05873 fHadr->MakeAntiHadron();
05874 sHadr->MakeAntiHadron();
05875 }
05876 G4LorentzVector f4Mom = fHadr->Get4Momentum();
05877 G4LorentzVector s4Mom = sHadr->Get4Momentum();
05878 if(!qH->DecayIn2(f4Mom,s4Mom))
05879 {
05880 delete fHadr;
05881 delete sHadr;
05882 #ifdef debug
05883 G4cerr<<"---Warning---G4Q::DecayQHadron:in2,PDGC="<<thePDG<<", ch#"<<i<<": 4M="
05884 <<qH->Get4Momentum()<<"("<<qH->GetMass()<<")->"<<f4Mom<<"+"<<s4Mom<<G4endl;
05885
05886 #endif
05887 theFragments->push_back(qH);
05888 return theFragments;
05889 }
05890 else
05891 {
05892
05893
05894
05895 delete qH;
05896
05897 fHadr->Set4Momentum(f4Mom);
05898 G4QHadronVector* theTmpQHV=DecayQHadron(fHadr);
05899 G4int nProd=theTmpQHV->size();
05900 #ifdef debug
05901 G4cout<<"G4Q::DecayQHadr:(DecayIn2) nOfProdForQH1="<<nProd<<G4endl;
05902 #endif
05903 if(nProd==1) theFragments->push_back((*theTmpQHV)[0]);
05904 else for(G4int ip1=0; ip1<nProd; ip1++)
05905 {
05906 G4QHadronVector* intTmpQHV = DecayQHadron((*theTmpQHV)[ip1]);
05907 G4int tmpS=intTmpQHV->size();
05908 if(tmpS==1)theFragments->push_back((*intTmpQHV)[0]);
05909 else
05910 {
05911 theFragments->resize(tmpS+theFragments->size());
05912 copy(intTmpQHV->begin(), intTmpQHV->end(), theFragments->end()-tmpS);
05913 }
05914 #ifdef debug
05915 G4cout<<"G4Q::DecayQHadr:(DecayIn2) Copy Sec11 nProd="<<tmpS<<G4endl;
05916 #endif
05917 intTmpQHV->clear();
05918 delete intTmpQHV;
05919 }
05920 theTmpQHV->clear();
05921 delete theTmpQHV;
05922 sHadr->Set4Momentum(s4Mom);
05923 theTmpQHV=DecayQHadron(sHadr);
05924 nProd=theTmpQHV->size();
05925 #ifdef debug
05926 G4cout<<"G4Q::DecayQHadr:(DecayIn2) nOfProdForQH2="<<nProd<<G4endl;
05927 #endif
05928 if(nProd==1) theFragments->push_back((*theTmpQHV)[0]);
05929 else for(G4int ip1=0; ip1<nProd; ip1++)
05930 {
05931 G4QHadronVector* intTmpQHV = DecayQHadron((*theTmpQHV)[ip1]);
05932 G4int tmpS=intTmpQHV->size();
05933 if(tmpS==1)theFragments->push_back((*intTmpQHV)[0]);
05934 else
05935 {
05936 theFragments->resize(tmpS+theFragments->size());
05937 copy(intTmpQHV->begin(), intTmpQHV->end(), theFragments->end()-tmpS);
05938 }
05939 #ifdef debug
05940 G4cout<<"G4Q::DecayQHadr:(DecayIn2) Copy Sec12 nProd="<<tmpS<<G4endl;
05941 #endif
05942 intTmpQHV->clear();
05943 delete intTmpQHV;
05944 }
05945 theTmpQHV->clear();
05946 delete theTmpQHV;
05947 }
05948 #ifdef debug
05949 G4cout<<"G4Q::DecQHadr: DecayIn2 is made with nH="<<theFragments->size()<<G4endl;
05950 #endif
05951 }
05952 else
05953 {
05954 #ifdef debug
05955 if(thePDG==89999003||thePDG==90002999)G4cerr<<"*G4Q::DQH:8999003/90002999"<<G4endl;
05956 #endif
05957 theFragments->push_back(qH);
05958 }
05959 }
05960 else if(nPart==3)
05961 {
05962 G4QHadron* fHadr;
05963 G4QHadron* sHadr;
05964 G4QHadron* tHadr;
05965 G4int fPDG=cV[0]->GetPDGCode();
05966 G4int sPDG=cV[1]->GetPDGCode();
05967 G4int tPDG=cV[2]->GetPDGCode();
05968
05969 if ( (fPDG != 22 && sPDG != 22 && tPDG != 22) || ElMaDecays)
05970 {
05971 #ifdef debug
05972 G4cout<<"G4Q::DQH:Y,f="<<fPDG<<",s="<<sPDG<<",t="<<tPDG<<",F="<<ElMaDecays<<G4endl;
05973 #endif
05974 if(cV[0]->GetWidth()==0.)
05975 {
05976 fHadr = new G4QHadron(fPDG);
05977 if(cV[1]->GetWidth()==0.)
05978 {
05979 sHadr = new G4QHadron(sPDG);
05980 if(cV[2]->GetWidth()==0.)tHadr = new G4QHadron(tPDG);
05981 else
05982 {
05983 G4QParticle* tPart=theWorld->GetQParticle(cV[2]);
05984 G4double tdm = m_value-fHadr->GetMass()-sHadr->GetMass();
05985 tHadr = new G4QHadron(tPart,tdm);
05986 if(tPDG<0) tHadr->MakeAntiHadron();
05987 }
05988 }
05989 else
05990 {
05991 m_value-=fHadr->GetMass();
05992 G4QParticle* tPart=theWorld->GetQParticle(cV[2]);
05993 G4double mim = tPart->MinMassOfFragm();
05994 G4double sdm = m_value - mim;
05995 G4QParticle* sPart=theWorld->GetQParticle(cV[1]);
05996 sHadr = new G4QHadron(sPart,sdm);
05997 if(sPDG<0) sHadr->MakeAntiHadron();
05998 G4double tdm = m_value - sHadr->GetMass();
05999 tHadr = new G4QHadron(tPart,tdm);
06000 if(tPDG<0) tHadr->MakeAntiHadron();
06001 }
06002 }
06003 else
06004 {
06005 G4QParticle* sPart=theWorld->GetQParticle(cV[1]);
06006 G4double smim = sPart->MinMassOfFragm();
06007 G4QParticle* tPart=theWorld->GetQParticle(cV[2]);
06008 G4double tmim = tPart->MinMassOfFragm();
06009 G4double fdm = m_value - smim - tmim;
06010 G4QParticle* fPart=theWorld->GetQParticle(cV[0]);
06011 fHadr = new G4QHadron(fPart,fdm);
06012 if(fPDG<0) fHadr->MakeAntiHadron();
06013 m_value-=fHadr->GetMass();
06014 G4double sdm = m_value - tmim;
06015 sHadr = new G4QHadron(sPart,sdm);
06016 if(sPDG<0) sHadr->MakeAntiHadron();
06017 G4double tdm = m_value - sHadr->GetMass();
06018 tHadr = new G4QHadron(tPart,tdm);
06019 if(tPDG<0) tHadr->MakeAntiHadron();
06020 }
06021 #ifdef debug
06022 G4cout<<"G4Quasmon::DecayQHadron:3Dec. m1="<<fHadr->GetMass()
06023 <<",m2="<<sHadr->GetMass()<<",m3="<<tHadr->GetMass()<<G4endl;
06024 #endif
06025 if(pap)
06026 {
06027 fHadr->MakeAntiHadron();
06028 sHadr->MakeAntiHadron();
06029 tHadr->MakeAntiHadron();
06030 }
06031 G4LorentzVector f4Mom = fHadr->Get4Momentum();
06032 G4LorentzVector s4Mom = sHadr->Get4Momentum();
06033 G4LorentzVector t4Mom = tHadr->Get4Momentum();
06034 if(!qH->DecayIn3(f4Mom,s4Mom,t4Mom))
06035 {
06036 delete fHadr;
06037 delete sHadr;
06038 delete tHadr;
06039 G4cerr<<"---Warning---G4Q::DecayQHadron:in3,PDGC="<<thePDG<<", ch#"<<i<<G4endl;
06040 theFragments->push_back(qH);
06041 return theFragments;
06042 }
06043 else
06044 {
06045
06046
06047
06048 delete qH;
06049
06050 fHadr->Set4Momentum(f4Mom);
06051 G4QHadronVector* theTmpQHV=DecayQHadron(fHadr);
06052 G4int nProd=theTmpQHV->size();
06053 #ifdef debug
06054 G4cout<<"G4Q::DecayQHadr:(DecayIn3) nOfProdForQH1="<<nProd<<G4endl;
06055 #endif
06056 if(nProd==1) theFragments->push_back((*theTmpQHV)[0]);
06057 else for(G4int ip1=0; ip1<nProd; ip1++)
06058 {
06059 G4QHadronVector* intTmpQHV = DecayQHadron((*theTmpQHV)[ip1]);
06060 G4int tmpS=intTmpQHV->size();
06061 if(tmpS==1)theFragments->push_back((*intTmpQHV)[0]);
06062 else
06063 {
06064 theFragments->resize(tmpS+theFragments->size());
06065 copy(intTmpQHV->begin(), intTmpQHV->end(), theFragments->end()-tmpS);
06066 }
06067 #ifdef debug
06068 G4cout<<"G4Q::DecayQHadr:(DecayIn3) Copy Sec11 nProd="<<tmpS<<G4endl;
06069 #endif
06070 intTmpQHV->clear();
06071 delete intTmpQHV;
06072 }
06073 theTmpQHV->clear();
06074 delete theTmpQHV;
06075
06076 sHadr->Set4Momentum(s4Mom);
06077 theTmpQHV=DecayQHadron(sHadr);
06078 nProd=theTmpQHV->size();
06079 #ifdef debug
06080 G4cout<<"G4Q::DecayQHadr:(DecayIn3) nOfProdForQH2="<<nProd<<G4endl;
06081 #endif
06082 if(nProd==1) theFragments->push_back((*theTmpQHV)[0]);
06083 else for(G4int ip1=0; ip1<nProd; ip1++)
06084 {
06085 G4QHadronVector* intTmpQHV = DecayQHadron((*theTmpQHV)[ip1]);
06086 G4int tmpS=intTmpQHV->size();
06087 if(tmpS==1)theFragments->push_back((*intTmpQHV)[0]);
06088 else
06089 {
06090 theFragments->resize(tmpS+theFragments->size());
06091 copy(intTmpQHV->begin(), intTmpQHV->end(), theFragments->end()-tmpS);
06092 }
06093 #ifdef debug
06094 G4cout<<"G4Q::DecayQHadr:(DecayIn3) Copy Sec12 nProd="<<tmpS<<G4endl;
06095 #endif
06096 intTmpQHV->clear();
06097 delete intTmpQHV;
06098 }
06099 theTmpQHV->clear();
06100 delete theTmpQHV;
06101
06102 tHadr->Set4Momentum(t4Mom);
06103 theTmpQHV=DecayQHadron(tHadr);
06104 nProd=theTmpQHV->size();
06105 #ifdef debug
06106 G4cout<<"G4Q::DecayQHadr:(DecayIn3) nOfProdForQH3="<<nProd<<G4endl;
06107 #endif
06108 if(nProd==1) theFragments->push_back((*theTmpQHV)[0]);
06109 else for(G4int ip1=0; ip1<nProd; ip1++)
06110 {
06111 G4QHadronVector* intTmpQHV = DecayQHadron((*theTmpQHV)[ip1]);
06112 G4int tmpS=intTmpQHV->size();
06113 if(tmpS==1)theFragments->push_back((*intTmpQHV)[0]);
06114 else
06115 {
06116 theFragments->resize(tmpS+theFragments->size());
06117 copy(intTmpQHV->begin(), intTmpQHV->end(), theFragments->end()-tmpS);
06118 }
06119 #ifdef debug
06120 G4cout<<"G4Q::DecayQHadr:(DecayIn3) Copy Sec13 nProd="<<tmpS<<G4endl;
06121 #endif
06122 intTmpQHV->clear();
06123 delete intTmpQHV;
06124 }
06125 theTmpQHV->clear();
06126 delete theTmpQHV;
06127
06128 }
06129 #ifdef debug
06130 G4cout<<"G4Q::DecQHadr: DecayIn3 is made with nH="<<theFragments->size()<<G4endl;
06131 #endif
06132 }
06133 else theFragments->push_back(qH);
06134 }
06135 }
06136 else
06137 {
06138 #ifdef debug
06139 G4cout<<"G4Quas::DecQHadr:Fill PDG= "<<thePDG<<t<<m_value<<" as it is ***0***>>"<<G4endl;
06140 #endif
06141 if(thePDG==89999003||thePDG==90002999)G4cerr<<"-War-G4Q::DQH:8999003/90002999"<<G4endl;
06142 theFragments->push_back(qH);
06143 }
06144 #ifdef debug
06145 G4cout<<"G4Q::DecQHadr:=-= HADRON IS DECAYED =-= with nH="<<theFragments->size()<<G4endl;
06146 #endif
06147 return theFragments;
06148 }
06149
06150
06151 G4int G4Quasmon::RandomPoisson(G4double meanValue)
06152 {
06153 if (meanValue<=0.)
06154 {
06155 G4cerr<<"---Warning---G4Q::RandomPoisson:Negative(zero) MeanValue="<<meanValue<<G4endl;
06156
06157 return -1;
06158 }
06159 G4double r=G4UniformRand();
06160 G4double t=exp(-meanValue);
06161 G4double s_value=t;
06162 if (r<s_value) return 0;
06163 t*=meanValue;
06164 s_value+=t;
06165 if (r<s_value) return 1;
06166 G4int i=1;
06167 while ( s_value<r && i<100 )
06168 {
06169 i++;
06170 t*=meanValue/i;
06171 s_value+=t;
06172 }
06173 return i;
06174 }
06175
06176
06177
06178 G4QHadronVector* G4Quasmon::Fragment(G4QNucleus& nucEnviron, G4int nQ)
06179 {
06180 #ifdef debug
06181 G4cout<<"G4Quasmon::Fragment called E="<<nucEnviron<<nucEnviron.GetProbability()<<G4endl;
06182 #endif
06183 G4int nQs=nQ;
06184 HadronizeQuasmon(nucEnviron,nQs);
06185 G4int nHadrs=theQHadrons.size();
06186 #ifdef debug
06187 G4cout<<"G4Quasm::Fragm:after HadronizeQuasmon nH="<<nHadrs<<",Env="<<nucEnviron<<G4endl;
06188 #endif
06189 G4QHadronVector* theFragments = new G4QHadronVector;
06190 if(nHadrs) for (int hadron=0; hadron<nHadrs; hadron++)
06191 {
06192 G4QHadron* curHadr = new G4QHadron(theQHadrons[hadron]);
06193 theFragments->push_back(curHadr);
06194 }
06195 #ifdef pdebug
06196 else G4cerr<<"*******G4Quasmon::Fragment *** Nothing is in the output ***"<<G4endl;
06197 #endif
06198 return theFragments;
06199 }
06200
06201
06202 void G4Quasmon::Boost(const G4LorentzVector& boost4M)
06203 {
06204
06205 G4double bm=boost4M.mag();
06206 G4double factor = (q4Mom.vect()*boost4M.vect()/(boost4M.e()+bm) - q4Mom.e())/bm;
06207 q4Mom.setE(q4Mom.dot(boost4M)/bm);
06208 q4Mom.setVect(factor*boost4M.vect() + q4Mom.vect());
06209 }