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 #include <cmath>
00046
00047 #include "G4QHadron.hh"
00048 #include "G4PhysicalConstants.hh"
00049 #include "G4SystemOfUnits.hh"
00050
00051 using namespace std;
00052
00053 G4double G4QHadron::StrangeSuppress = 0.48;
00054 G4double G4QHadron::sigmaPt = 1.7*GeV;
00055 G4double G4QHadron::widthOfPtSquare = 0.01*GeV*GeV;
00056
00057 G4QHadron::G4QHadron(): theMomentum(0.,0.,0.,0.), theQPDG(0), valQ(0,0,0,0,0,0), nFragm(0),
00058 thePosition(0.,0.,0.), theCollisionCount(0), isSplit(false), Direction(true),
00059 Color(), AntiColor(), bindE(0.), formTime(0.) {}
00060
00061 G4QHadron::G4QHadron(G4LorentzVector p): theMomentum(p), theQPDG(0), valQ(0,0,0,0,0,0),
00062 nFragm(0), thePosition(0.,0.,0.), theCollisionCount(0), isSplit(false), Direction(true),
00063 Color(), AntiColor(), bindE(0.), formTime(0.) {}
00064
00065
00066 G4QHadron::G4QHadron(G4int PDGCode, G4LorentzVector p): theMomentum(p), theQPDG(PDGCode),
00067 nFragm(0),thePosition(0.,0.,0.),theCollisionCount(0),isSplit(false),Direction(true),
00068 Color(), AntiColor(), bindE(0.), formTime(0.)
00069 {
00070 #ifdef debug
00071 G4cout<<"G4QHadron must be created with PDG="<<PDGCode<<", 4M="<<p<<G4endl;
00072 #endif
00073 if(GetQCode()>-1)
00074 {
00075 if(theMomentum.e()==0.) theMomentum.setE(theQPDG.GetMass());
00076 valQ=theQPDG.GetQuarkContent();
00077 }
00078 else if(PDGCode>80000000) DefineQC(PDGCode);
00079 else G4cerr<<"***G4QHadron:(P) PDG="<<PDGCode<<", use other constructor"<<G4endl;
00080 #ifdef debug
00081 G4cout<<"G4QHadron is created with QCode="<<GetQCode()<<", QC="<<valQ<<G4endl;
00082 #endif
00083 }
00084
00085
00086 G4QHadron::G4QHadron(G4QPDGCode QPDG, G4LorentzVector p): theMomentum(p), theQPDG(QPDG),
00087 nFragm(0), thePosition(0.,0.,0.), theCollisionCount(0), isSplit(false), Direction(true),
00088 Color(), AntiColor(), bindE(0.), formTime(0.)
00089 {
00090 if(theQPDG.GetQCode()>-1)
00091 {
00092 if(theMomentum.e()==0.) theMomentum.setE(theQPDG.GetMass());
00093 valQ=theQPDG.GetQuarkContent();
00094 }
00095 else
00096 {
00097 G4int cPDG=theQPDG.GetPDGCode();
00098 if(cPDG>80000000) DefineQC(cPDG);
00099 else G4cerr<<"***G4QHadr:(QP) PDG="<<cPDG<<" use other constructor"<<G4endl;
00100 }
00101 }
00102
00103
00104 G4QHadron::G4QHadron(G4QContent QC, G4LorentzVector p): theMomentum(p),theQPDG(0),valQ(QC),
00105 nFragm(0), thePosition(0.,0.,0.), theCollisionCount(0), isSplit(false), Direction(true),
00106 Color(), AntiColor(), bindE(0.), formTime(0.)
00107 {
00108 G4int curPDG=valQ.GetSPDGCode();
00109 if(curPDG==10&&valQ.GetBaryonNumber()>0) curPDG=valQ.GetZNSPDGCode();
00110 if(curPDG&&curPDG!=10) theQPDG.SetPDGCode(curPDG);
00111 else theQPDG.InitByQCont(QC);
00112 }
00113
00114 G4QHadron::G4QHadron(G4int PDGCode, G4double aMass, G4QContent QC) :
00115 theMomentum(0.,0.,0.,aMass), theQPDG(PDGCode), valQ(QC), nFragm(0),thePosition(0.,0.,0.),
00116 theCollisionCount(0), isSplit(false), Direction(true), Color(), AntiColor(), bindE(0.),
00117 formTime(0.)
00118 {}
00119
00120 G4QHadron::G4QHadron(G4QPDGCode QPDG, G4double aMass, G4QContent QC) :
00121 theMomentum(0.,0.,0.,aMass), theQPDG(QPDG), valQ(QC), nFragm(0), thePosition(0.,0.,0.),
00122 theCollisionCount(0), isSplit(false), Direction(true), Color(), AntiColor(), bindE(0.),
00123 formTime(0.)
00124 {}
00125
00126 G4QHadron::G4QHadron(G4int PDGCode, G4LorentzVector p, G4QContent QC) : theMomentum(p),
00127 theQPDG(PDGCode), valQ(QC), nFragm(0), thePosition(0.,0.,0.), theCollisionCount(0),
00128 isSplit(false), Direction(true), Color(), AntiColor(), bindE(0.), formTime(0.)
00129 {}
00130
00131 G4QHadron::G4QHadron(G4QPDGCode QPDG, G4LorentzVector p, G4QContent QC) : theMomentum(p),
00132 theQPDG(QPDG), valQ(QC), nFragm(0), thePosition(0.,0.,0.), theCollisionCount(0),
00133 isSplit(false), Direction(true), Color(), AntiColor(), bindE(0.), formTime(0.)
00134 {}
00135
00136 G4QHadron::G4QHadron(G4QParticle* pPart, G4double maxM) : theMomentum(0.,0.,0.,0.),
00137 theQPDG(pPart->GetQPDG()), nFragm(0), thePosition(0.,0.,0.), theCollisionCount(0),
00138 isSplit(false), Direction(true), Color(), AntiColor(), bindE(0.), formTime(0.)
00139 {
00140 #ifdef debug
00141 G4cout<<"G4QHadron is created & randomized with maxM="<<maxM<<G4endl;
00142 #endif
00143 G4int PDGCode = theQPDG.GetPDGCode();
00144 if(PDGCode<2)G4cerr<<"***G4QHadron:(M) PDGC="<<PDGCode<<" use other constructor"<<G4endl;
00145 valQ=theQPDG.GetQuarkContent();
00146 theMomentum.setE(RandomizeMass(pPart, maxM));
00147 }
00148
00149 G4QHadron::G4QHadron(const G4QHadron& right)
00150 {
00151 theMomentum = right.theMomentum;
00152 theQPDG = right.theQPDG;
00153 valQ = right.valQ;
00154 nFragm = right.nFragm;
00155 thePosition = right.thePosition;
00156 theCollisionCount = 0;
00157 isSplit = false;
00158 Direction = right.Direction;
00159 bindE = right.bindE;
00160 formTime = right.formTime;
00161 }
00162
00163 G4QHadron::G4QHadron(const G4QHadron* right)
00164 {
00165 theMomentum = right->theMomentum;
00166 theQPDG = right->theQPDG;
00167 valQ = right->valQ;
00168 nFragm = right->nFragm;
00169 thePosition = right->thePosition;
00170 theCollisionCount = 0;
00171 isSplit = false;
00172 Direction = right->Direction;
00173 bindE = right->bindE;
00174 formTime = right->formTime;
00175 }
00176
00177 G4QHadron::G4QHadron(const G4QHadron* right, G4int C, G4ThreeVector P, G4LorentzVector M)
00178 {
00179 theMomentum = M;
00180 theQPDG = right->theQPDG;
00181 valQ = right->valQ;
00182 nFragm = right->nFragm;
00183 thePosition = P;
00184 theCollisionCount = C;
00185 isSplit = false;
00186 Direction = right->Direction;
00187 bindE = right->bindE;
00188 formTime = right->formTime;
00189 }
00190
00191 const G4QHadron& G4QHadron::operator=(const G4QHadron &right)
00192 {
00193 if(this != &right)
00194 {
00195 theMomentum = right.theMomentum;
00196 theQPDG = right.theQPDG;
00197 valQ = right.valQ;
00198 nFragm = right.nFragm;
00199 thePosition = right.thePosition;
00200 theCollisionCount = 0;
00201 isSplit = false;
00202 Direction = right.Direction;
00203 bindE = right.bindE;
00204 }
00205 return *this;
00206 }
00207
00208 G4QHadron::~G4QHadron()
00209 {
00210 std::list<G4QParton*>::iterator ipos = Color.begin();
00211 std::list<G4QParton*>::iterator epos = Color.end();
00212 for( ; ipos != epos; ipos++) {delete [] *ipos;}
00213 Color.clear();
00214
00215 ipos = AntiColor.begin();
00216 epos = AntiColor.end();
00217 for( ; ipos != epos; ipos++) {delete [] *ipos;}
00218 AntiColor.clear();
00219 }
00220
00221
00222 void G4QHadron::DefineQC(G4int PDGCode)
00223 {
00224
00225 G4int szn=PDGCode-90000000;
00226 G4int ds=0;
00227 G4int dz=0;
00228 G4int dn=0;
00229 if(szn<-100000)
00230 {
00231 G4int ns_value=(-szn)/1000000+1;
00232 szn+=ns_value*1000000;
00233 ds+=ns_value;
00234 }
00235 else if(szn<-100)
00236 {
00237 G4int nz=(-szn)/1000+1;
00238 szn+=nz*1000;
00239 dz+=nz;
00240 }
00241 else if(szn<0)
00242 {
00243 G4int nn=-szn;
00244 szn=0;
00245 dn+=nn;
00246 }
00247 G4int sz =szn/1000;
00248 G4int n =szn%1000;
00249 if(n>700)
00250 {
00251 n-=1000;
00252 dz--;
00253 }
00254 G4int z =sz%1000-dz;
00255 if(z>700)
00256 {
00257 z-=1000;
00258 ds--;
00259 }
00260 G4int Sq =sz/1000-ds;
00261 G4int zns=z+n+Sq;
00262 G4int Dq=n+zns;
00263 G4int Uq=z+zns;
00264 if (Dq<0&&Uq<0&&Sq<0)valQ=G4QContent(0 ,0 ,0 ,-Dq,-Uq,-Sq);
00265 else if (Uq<0&&Sq<0) valQ=G4QContent(Dq,0 ,0 ,0 ,-Uq,-Sq);
00266 else if (Dq<0&&Sq<0) valQ=G4QContent(0 ,Uq,0 ,-Dq,0 ,-Sq);
00267 else if (Dq<0&&Uq<0) valQ=G4QContent(0 ,0 ,Sq,-Dq,-Uq,0 );
00268 else if (Uq<0) valQ=G4QContent(Dq,0 ,Sq,0 ,-Uq,0 );
00269 else if (Sq<0) valQ=G4QContent(Dq,Uq,0 ,0 ,0 ,-Sq);
00270 else if (Dq<0) valQ=G4QContent(0 ,Uq,Sq,-Dq,0 ,0 );
00271 else valQ=G4QContent(Dq,Uq,Sq,0 ,0 ,0 );
00272 }
00273
00274
00275 void G4QHadron::SetQPDG(const G4QPDGCode& newQPDG)
00276 {
00277 theQPDG = newQPDG;
00278 G4int PDG= newQPDG.GetPDGCode();
00279 G4int Q = newQPDG.GetQCode();
00280 #ifdef debug
00281 G4cout<<"G4QHadron::SetQPDG is called with PDGCode="<<PDG<<", QCode="<<Q<<G4endl;
00282 #endif
00283 if (Q>-1) valQ=theQPDG.GetQuarkContent();
00284 else if(PDG>80000000) DefineQC(PDG);
00285 else
00286 {
00287
00288
00289 G4ExceptionDescription ed;
00290 ed << "Impossible QPDG Probably a Chipolino: QPDG=" << newQPDG << G4endl;
00291 G4Exception("G4QHadron::SetQPDG()", "HAD_CHPS_0000", FatalException, ed);
00292 }
00293 }
00294
00295
00296 G4bool G4QHadron::RelDecayIn2(G4LorentzVector& f4Mom, G4LorentzVector& s4Mom,
00297 G4LorentzVector& dir, G4double maxCost, G4double minCost)
00298 {
00299 G4double fM2 = f4Mom.m2();
00300 G4double fM = sqrt(fM2);
00301 G4double sM2 = s4Mom.m2();
00302 G4double sM = sqrt(sM2);
00303 G4double iM2 = theMomentum.m2();
00304 G4double iM = sqrt(iM2);
00305 G4double vP = theMomentum.rho();
00306 G4double dE = theMomentum.e();
00307 if(dE<vP)
00308 {
00309 G4cerr<<"***G4QHad::RelDecIn2: Tachionic 4-mom="<<theMomentum<<", E-p="<<dE-vP<<G4endl;
00310 G4double accuracy=.000001*vP;
00311 G4double emodif=std::fabs(dE-vP);
00312
00313
00314 G4cerr<<"G4QHadron::RelDecIn2: *Boost* E-p shift is corrected to "<<emodif<<G4endl;
00315 theMomentum.setE(vP+emodif+.01*accuracy);
00316
00317 }
00318 G4ThreeVector ltb = theMomentum.boostVector();
00319 G4ThreeVector ltf = -ltb;
00320 G4LorentzVector cdir = dir;
00321 #ifdef ppdebug
00322 if(cdir.e()+.001<cdir.rho()) G4cerr<<"*G4QH::RDIn2:*Boost* cd4M="<<cdir<<",e-p="
00323 <<cdir.e()-cdir.rho()<<G4endl;
00324 #endif
00325 cdir.boost(ltf);
00326 G4ThreeVector vdir = cdir.vect();
00327 #ifdef ppdebug
00328 G4cout<<"G4QHad::RelDI2:dir="<<dir<<",ltf="<<ltf<<",cdir="<<cdir<<",vdir="<<vdir<<G4endl;
00329 #endif
00330 G4ThreeVector vx(0.,0.,1.);
00331 G4ThreeVector vy(0.,1.,0.);
00332 G4ThreeVector vz(1.,0.,0.);
00333 if(vdir.mag2() > 0.)
00334 {
00335 vx = vdir.unit();
00336 #ifdef ppdebug
00337 G4cout<<"G4QH::RelDecIn2:Vx="<<vx<<",M="<<theMomentum<<",d="<<dir<<",c="<<cdir<<G4endl;
00338 #endif
00339 G4ThreeVector vv= vx.orthogonal();
00340 vy = vv.unit();
00341 vz = vx.cross(vy);
00342 }
00343 #ifdef ppdebug
00344 G4cout<<"G4QHad::RelDecIn2:iM="<<iM<<"=>fM="<<fM<<"+sM="<<sM<<",ob="<<vx<<vy<<vz<<G4endl;
00345 #endif
00346 if(maxCost> 1.) maxCost= 1.;
00347 if(minCost<-1.) minCost=-1.;
00348 if(maxCost<-1.) maxCost=-1.;
00349 if(minCost> 1.) minCost= 1.;
00350 if(minCost> maxCost) minCost=maxCost;
00351 if(fabs(iM-fM-sM)<.00000001)
00352 {
00353 G4double fR=fM/iM;
00354 G4double sR=sM/iM;
00355 f4Mom=fR*theMomentum;
00356 s4Mom=sR*theMomentum;
00357 return true;
00358 }
00359 else if (iM+.001<fM+sM || iM==0.)
00360 {
00361 G4cerr<<"***G4QH::RelDecIn2: fM="<<fM<<"+sM="<<sM<<">iM="<<iM<<",d="<<iM-fM-sM<<G4endl;
00362 return false;
00363 }
00364 G4double d2 = iM2-fM2-sM2;
00365 G4double p2 = (d2*d2/4.-fM2*sM2)/iM2;
00366 if(p2<0.)
00367 {
00368 #ifdef ppdebug
00369 G4cout<<"**G4QH:RDIn2:p2="<<p2<<"<0,d2^2="<<d2*d2/4.<<"<4*fM2*sM2="<<4*fM2*sM2<<G4endl;
00370 #endif
00371 p2=0.;
00372 }
00373 G4double p = sqrt(p2);
00374 G4double ct = maxCost;
00375 if(maxCost>minCost)
00376 {
00377 G4double dcost=maxCost-minCost;
00378 ct = minCost+dcost*G4UniformRand();
00379 }
00380 G4double phi= twopi*G4UniformRand();
00381 G4double ps=0.;
00382 if(fabs(ct)<1.) ps = p * sqrt(1.-ct*ct);
00383 else
00384 {
00385 #ifdef ppdebug
00386 G4cout<<"**G4QH::RDIn2:ct="<<ct<<",mac="<<maxCost<<",mic="<<minCost<<G4endl;
00387
00388 #endif
00389 if(ct>1.) ct=1.;
00390 if(ct<-1.) ct=-1.;
00391 }
00392 G4ThreeVector pVect=(ps*sin(phi))*vz+(ps*cos(phi))*vy+p*ct*vx;
00393 #ifdef ppdebug
00394 G4cout<<"G4QH::RelDIn2:ct="<<ct<<",p="<<p<<",ps="<<ps<<",ph="<<phi<<",v="<<pVect<<G4endl;
00395 #endif
00396
00397 f4Mom.setVect(pVect);
00398 f4Mom.setE(sqrt(fM2+p2));
00399 s4Mom.setVect((-1)*pVect);
00400 s4Mom.setE(sqrt(sM2+p2));
00401
00402 #ifdef ppdebug
00403 G4cout<<"G4QHadr::RelDecIn2:p2="<<p2<<",v="<<ltb<<",f4M="<<f4Mom<<" + s4M="<<s4Mom<<" = "
00404 <<f4Mom+s4Mom<<", M="<<iM<<G4endl;
00405 #endif
00406 if(f4Mom.e()+.001<f4Mom.rho())G4cerr<<"*G4QH::RDIn2:*Boost* f4M="<<f4Mom<<",e-p="
00407 <<f4Mom.e()-f4Mom.rho()<<G4endl;
00408 f4Mom.boost(ltb);
00409 if(s4Mom.e()+.001<s4Mom.rho())G4cerr<<"*G4QH::RDIn2:*Boost* s4M="<<s4Mom<<",e-p="
00410 <<s4Mom.e()-s4Mom.rho()<<G4endl;
00411 s4Mom.boost(ltb);
00412 #ifdef ppdebug
00413 G4cout<<"G4QHadron::RelDecayIn2:Output, f4Mom="<<f4Mom<<" + s4Mom="<<s4Mom<<" = "
00414 <<f4Mom+s4Mom<<", d4M="<<theMomentum-f4Mom-s4Mom<<G4endl;
00415 #endif
00416 return true;
00417 }
00418
00419
00420 G4bool G4QHadron::CopDecayIn2(G4LorentzVector& f4Mom, G4LorentzVector& s4Mom,
00421 G4LorentzVector& dir, G4double cosp)
00422 {
00423 G4double fM2 = f4Mom.m2();
00424 G4double fM = sqrt(fM2);
00425 G4double sM2 = s4Mom.m2();
00426 G4double sM = sqrt(sM2);
00427 G4double iM2 = theMomentum.m2();
00428 G4double iM = sqrt(iM2);
00429 G4double vP = theMomentum.rho();
00430 G4double dE = theMomentum.e();
00431 G4bool neg=false;
00432 if(cosp<0)
00433 {
00434 cosp=-cosp;
00435 neg=true;
00436 }
00437 if(dE<vP)
00438 {
00439 G4cerr<<"***G4QHad::CopDecIn2: Tachionic 4-mom="<<theMomentum<<", E-p="<<dE-vP<<G4endl;
00440 G4double accuracy=.000001*vP;
00441 G4double emodif=std::fabs(dE-vP);
00442
00443
00444 G4cerr<<"G4QHadron::CopDecIn2: *Boost* E-p shift is corrected to "<<emodif<<G4endl;
00445 theMomentum.setE(vP+emodif+.01*accuracy);
00446
00447 }
00448 G4ThreeVector ltb = theMomentum.boostVector();
00449 G4ThreeVector ltf = -ltb;
00450 G4LorentzVector cdir = dir;
00451 #ifdef ppdebug
00452 if(cdir.e()+.001<cdir.rho()) G4cerr<<"*G4QH::RDIn2:*Boost* cd4M="<<cdir<<",e-p="
00453 <<cdir.e()-cdir.rho()<<G4endl;
00454 #endif
00455 cdir.boost(ltf);
00456 G4ThreeVector vdir = cdir.vect();
00457 #ifdef ppdebug
00458 G4cout<<"G4QHad::CopDI2:dir="<<dir<<",ltf="<<ltf<<",cdir="<<cdir<<",vdir="<<vdir<<G4endl;
00459 #endif
00460 G4ThreeVector vx(0.,0.,1.);
00461 G4ThreeVector vy(0.,1.,0.);
00462 G4ThreeVector vz(1.,0.,0.);
00463 if(vdir.mag2() > 0.)
00464 {
00465 vx = vdir.unit();
00466 #ifdef ppdebug
00467 G4cout<<"G4QH::CopDecIn2:Vx="<<vx<<",M="<<theMomentum<<",d="<<dir<<",c="<<cdir<<G4endl;
00468 #endif
00469 G4ThreeVector vv= vx.orthogonal();
00470 vy = vv.unit();
00471 vz = vx.cross(vy);
00472 }
00473 #ifdef ppdebug
00474 G4cout<<"G4QHad::CopDecIn2:iM="<<iM<<"=>fM="<<fM<<"+sM="<<sM<<",ob="<<vx<<vy<<vz<<G4endl;
00475 #endif
00476 if(fabs(iM-fM-sM)<.00000001)
00477 {
00478 G4double fR=fM/iM;
00479 G4double sR=sM/iM;
00480 f4Mom=fR*theMomentum;
00481 s4Mom=sR*theMomentum;
00482 return true;
00483 }
00484 else if (iM+.001<fM+sM || iM==0.)
00485 {
00486 G4cerr<<"***G4QH::CopDecIn2: fM="<<fM<<"+sM="<<sM<<">iM="<<iM<<",d="<<iM-fM-sM<<G4endl;
00487 return false;
00488 }
00489 G4double d2 = iM2-fM2-sM2;
00490 G4double p2 = (d2*d2/4.-fM2*sM2)/iM2;
00491 if(p2<0.)
00492 {
00493 #ifdef ppdebug
00494 G4cout<<"*G4QH:CopDI2:p2="<<p2<<"<0,d4/4="<<d2*d2/4.<<"<4*fM2*sM2="<<4*fM2*sM2<<G4endl;
00495 #endif
00496 p2=0.;
00497 }
00498 G4double p = sqrt(p2);
00499 G4double ct = 0;
00500 G4double rn = pow(G4UniformRand(),cosp+1.);
00501 if(neg) ct = rn+rn-1.;
00502 else ct = 1.-rn-rn;
00503
00504 G4double phi= twopi*G4UniformRand();
00505 G4double ps=0.;
00506 if(fabs(ct)<1.) ps = p * sqrt(1.-ct*ct);
00507 else
00508 {
00509 #ifdef ppdebug
00510 G4cout<<"**G4QH::CopDecayIn2:ct="<<ct<<",mac="<<maxCost<<",mic="<<minCost<<G4endl;
00511
00512 #endif
00513 if(ct>1.) ct=1.;
00514 if(ct<-1.) ct=-1.;
00515 }
00516 G4ThreeVector pVect=(ps*sin(phi))*vz+(ps*cos(phi))*vy+p*ct*vx;
00517 #ifdef ppdebug
00518 G4cout<<"G4QH::CopDIn2:ct="<<ct<<",p="<<p<<",ps="<<ps<<",ph="<<phi<<",v="<<pVect<<G4endl;
00519 #endif
00520
00521 f4Mom.setVect(pVect);
00522 f4Mom.setE(sqrt(fM2+p2));
00523 s4Mom.setVect((-1)*pVect);
00524 s4Mom.setE(sqrt(sM2+p2));
00525
00526 #ifdef ppdebug
00527 G4cout<<"G4QHadr::CopDecIn2:p2="<<p2<<",v="<<ltb<<",f4M="<<f4Mom<<" + s4M="<<s4Mom<<" = "
00528 <<f4Mom+s4Mom<<", M="<<iM<<G4endl;
00529 #endif
00530 if(f4Mom.e()+.001<f4Mom.rho())G4cerr<<"*G4QH::RDIn2:*Boost* f4M="<<f4Mom<<",e-p="
00531 <<f4Mom.e()-f4Mom.rho()<<G4endl;
00532 f4Mom.boost(ltb);
00533 if(s4Mom.e()+.001<s4Mom.rho())G4cerr<<"*G4QH::RDIn2:*Boost* s4M="<<s4Mom<<",e-p="
00534 <<s4Mom.e()-s4Mom.rho()<<G4endl;
00535 s4Mom.boost(ltb);
00536 #ifdef ppdebug
00537 G4cout<<"G4QHadron::CopDecayIn2:Output, f4Mom="<<f4Mom<<" + s4Mom="<<s4Mom<<" = "
00538 <<f4Mom+s4Mom<<", d4M="<<theMomentum-f4Mom-s4Mom<<G4endl;
00539 #endif
00540 return true;
00541 }
00542
00543
00544 G4bool G4QHadron::DecayIn2(G4LorentzVector& f4Mom, G4LorentzVector& s4Mom)
00545 {
00546 G4double fM2 = f4Mom.m2();
00547 if(fM2<0.) fM2=0.;
00548 G4double fM = sqrt(fM2);
00549 G4double sM2 = s4Mom.m2();
00550 if(sM2<0.) sM2=0.;
00551 G4double sM = sqrt(sM2);
00552 G4double iM2 = theMomentum.m2();
00553 if(iM2<0.) iM2=0.;
00554 G4double iM = sqrt(iM2);
00555 #ifdef debug
00556 G4cout<<"G4QHadron::DecIn2: iM="<<iM<<" => fM="<<fM<<" + sM="<<sM<<" = "<<fM+sM<<G4endl;
00557 #endif
00558
00559 if (fabs(iM-fM-sM)<.0000001)
00560 {
00561 G4double fR=fM/iM;
00562 G4double sR=sM/iM;
00563 f4Mom=fR*theMomentum;
00564 s4Mom=sR*theMomentum;
00565 return true;
00566 }
00567 else if (iM+.001<fM+sM || iM==0.)
00568 {
00569 #ifdef debug
00570 G4cerr<<"***G4QHadron::DecayIn2*** fM="<<fM<<" + sM="<<sM<<"="<<fM+sM<<" > iM="<<iM
00571 <<", d="<<iM-fM-sM<<G4endl;
00572 #endif
00573 return false;
00574 }
00575
00576 G4double d2 = iM2-fM2-sM2;
00577 G4double p2 = (d2*d2/4.-fM2*sM2)/iM2;
00578 if (p2<0.)
00579 {
00580 #ifdef debug
00581 G4cerr<<"***G4QH::DI2:p2="<<p2<<"<0,d2^2="<<d2*d2/4.<<"<4*fM2*sM2="<<4*fM2*sM2<<G4endl;
00582 #endif
00583 p2=0.;
00584 }
00585 G4double p = sqrt(p2);
00586 G4double ct = 1.-2*G4UniformRand();
00587 #ifdef debug
00588 G4cout<<"G4QHadron::DecayIn2: ct="<<ct<<", p="<<p<<G4endl;
00589 #endif
00590 G4double phi= twopi*G4UniformRand();
00591 G4double ps = p * sqrt(1.-ct*ct);
00592 G4ThreeVector pVect(ps*sin(phi),ps*cos(phi),p*ct);
00593
00594 f4Mom.setVect(pVect);
00595 f4Mom.setE(sqrt(fM2+p2));
00596 s4Mom.setVect((-1)*pVect);
00597 s4Mom.setE(sqrt(sM2+p2));
00598
00599 if(theMomentum.e()<theMomentum.rho())
00600 {
00601 G4cerr<<"*G4QH::DecIn2:*Boost* 4M="<<theMomentum<<",e-p="
00602 <<theMomentum.e()-theMomentum.rho()<<G4endl;
00603
00604
00605 return false;
00606 }
00607 G4double vP = theMomentum.rho();
00608 G4double dE = theMomentum.e();
00609 if(dE<vP)
00610 {
00611 G4cerr<<"***G4QHad::RelDecIn2: Tachionic 4-mom="<<theMomentum<<", E-p="<<dE-vP<<G4endl;
00612 G4double accuracy=.000001*vP;
00613 G4double emodif=std::fabs(dE-vP);
00614 if(emodif<accuracy)
00615 {
00616 G4cerr<<"G4QHadron::DecayIn2: *Boost* E-p shift is corrected to "<<emodif<<G4endl;
00617 theMomentum.setE(vP+emodif+.01*accuracy);
00618 }
00619 }
00620 G4ThreeVector ltb = theMomentum.boostVector();
00621 #ifdef debug
00622 G4cout<<"G4QHadron::DecIn2:LorTrans v="<<ltb<<",f4Mom="<<f4Mom<<",s4Mom="<<s4Mom<<G4endl;
00623 #endif
00624 if(f4Mom.e()+.001<f4Mom.rho())G4cerr<<"*G4QH::DecIn2:*Boost* f4M="<<f4Mom<<G4endl;
00625 f4Mom.boost(ltb);
00626 if(s4Mom.e()+.001<s4Mom.rho())G4cerr<<"*G4QH::DecIn2:*Boost* s4M="<<s4Mom<<G4endl;
00627 s4Mom.boost(ltb);
00628 #ifdef debug
00629 G4cout<<"G4QHadron::DecayIn2: ROOT OUTPUT f4Mom="<<f4Mom<<", s4Mom="<<s4Mom<<G4endl;
00630 #endif
00631 return true;
00632 }
00633
00634
00635 G4bool G4QHadron::CorMDecayIn2(G4double corM, G4LorentzVector& fr4Mom)
00636 {
00637 G4double fM = fr4Mom.m();
00638 G4LorentzVector comp=theMomentum+fr4Mom;
00639 G4double iM = comp.m();
00640 #ifdef debug
00641 G4cout<<"G4QH::CMDIn2: iM="<<iM<<comp<<"=>fM="<<fM<<"+corM="<<corM<<"="<<fM+corM<<G4endl;
00642 #endif
00643 G4double dE=iM-fM-corM;
00644
00645 if (fabs(dE)<.001)
00646 {
00647 G4double fR=fM/iM;
00648 G4double cR=corM/iM;
00649 fr4Mom=fR*comp;
00650 theMomentum=cR*comp;
00651 return true;
00652 }
00653 else if (dE<-.001 || iM==0.)
00654 {
00655 G4cerr<<"***G4QH::CorMDIn2***fM="<<fM<<" + cM="<<corM<<" > iM="<<iM<<",d="<<dE<<G4endl;
00656 return false;
00657 }
00658 G4double corM2= corM*corM;
00659 G4double fM2 = fM*fM;
00660 G4double iM2 = iM*iM;
00661 G4double d2 = iM2-fM2-corM2;
00662 G4double p2 = (d2*d2/4.-fM2*corM2)/iM2;
00663 if (p2<0.)
00664 {
00665 #ifdef debug
00666 G4cerr<<"**G4QH::CMDI2:p2="<<p2<<"<0,d="<<d2*d2/4.<<"<4*fM2*hM2="<<4*fM2*corM2<<G4endl;
00667 #endif
00668 p2=0.;
00669 }
00670 G4double p = sqrt(p2);
00671 if(comp.e()<comp.rho())G4cerr<<"*G4QH::CorMDecayIn2:*Boost* comp4M="<<comp<<",e-p="
00672 <<comp.e()-comp.rho()<<G4endl;
00673 G4ThreeVector ltb = comp.boostVector();
00674 G4ThreeVector ltf = -ltb;
00675 G4LorentzVector cm4Mom=fr4Mom;
00676 if(cm4Mom.e()<cm4Mom.rho())
00677 {
00678 G4cerr<<"*G4QH::CorMDecIn2:*Boost* c4M="<<cm4Mom<<G4endl;
00679
00680 return false;
00681 }
00682 cm4Mom.boost(ltf);
00683 G4double pfx= cm4Mom.px();
00684 G4double pfy= cm4Mom.py();
00685 G4double pfz= cm4Mom.pz();
00686 G4double pt2= pfx*pfx+pfy*pfy;
00687 G4double tx=0.;
00688 G4double ty=0.;
00689 if(pt2<=0.)
00690 {
00691 G4double phi= 360.*deg*G4UniformRand();
00692 tx=sin(phi);
00693 ty=cos(phi);
00694 }
00695 else
00696 {
00697 G4double pt=sqrt(pt2);
00698 tx=pfx/pt;
00699 ty=pfy/pt;
00700 }
00701 G4double pc2=pt2+pfz*pfz;
00702 G4double ct=0.;
00703 if(pc2<=0.)
00704 {
00705 G4double rnd= G4UniformRand();
00706 ct=1.-rnd-rnd;
00707 }
00708 else
00709 {
00710 G4double pc_value=sqrt(pc2);
00711 ct=pfz/pc_value;
00712 }
00713 #ifdef debug
00714 G4cout<<"G4QHadron::CorMDecayIn2: ct="<<ct<<", p="<<p<<G4endl;
00715 #endif
00716 G4double ps = p * sqrt(1.-ct*ct);
00717 G4ThreeVector pVect(ps*tx,ps*ty,p*ct);
00718 fr4Mom.setVect(pVect);
00719 fr4Mom.setE(sqrt(fM2+p2));
00720 theMomentum.setVect((-1)*pVect);
00721 theMomentum.setE(sqrt(corM2+p2));
00722 #ifdef debug
00723 G4LorentzVector dif2=comp-fr4Mom-theMomentum;
00724 G4cout<<"G4QH::CorMDIn2:c="<<comp<<"-f="<<fr4Mom<<"-4M="<<theMomentum<<"="<<dif2<<G4endl;
00725 #endif
00726 if(fr4Mom.e()+.001<fr4Mom.rho())G4cerr<<"*G4QH::CorMDecIn2:*Boost*fr4M="<<fr4Mom<<G4endl;
00727 fr4Mom.boost(ltb);
00728 if(theMomentum.e()<theMomentum.rho())
00729 {
00730 G4cerr<<"*G4QH::CMDI2:4="<<theMomentum<<G4endl;
00731 theMomentum.setE(1.0000001*theMomentum.rho());
00732 }
00733 theMomentum.boost(ltb);
00734 #ifdef debug
00735 G4LorentzVector dif3=comp-fr4Mom-theMomentum;
00736 G4cout<<"G4QH::CorMDecIn2:OUTPUT:f4M="<<fr4Mom<<",h4M="<<theMomentum<<"d="<<dif3<<G4endl;
00737 #endif
00738 return true;
00739 }
00740
00741
00742
00743 G4bool G4QHadron::CorEDecayIn2(G4double corE, G4LorentzVector& fr4Mom)
00744 {
00745 G4double fE = fr4Mom.m();
00746 #ifdef debug
00747 G4cout<<"G4QH::CorEDecIn2:fE="<<fE<<fr4Mom<<">corE="<<corE<<",h4M="<<theMomentum<<G4endl;
00748 #endif
00749 if (fE+.001<=corE)
00750 {
00751 #ifdef debug
00752 G4cerr<<"***G4QHadron::CorEDecIn2*** fE="<<fE<<"<corE="<<corE<<", d="<<corE-fE<<G4endl;
00753 #endif
00754 return false;
00755 }
00756 G4double fM2=fr4Mom.m2();
00757 if(fM2<0.) fM2=0.;
00758 G4double iPx=fr4Mom.px();
00759 G4double iPy=fr4Mom.py();
00760 G4double iPz=fr4Mom.pz();
00761 G4double fP2=iPx*iPx+iPy*iPy+iPz*iPz;
00762 G4double finE = fE - corE;
00763 G4double rP = sqrt((finE*finE-fM2)/fP2);
00764 G4double fPx=iPx*rP;
00765 G4double fPy=iPy*rP;
00766 G4double fPz=iPz*rP;
00767 fr4Mom= G4LorentzVector(fPx,fPy,fPz,finE);
00768 G4double Px=theMomentum.px()+iPx-fPx;
00769 G4double Py=theMomentum.py()+iPy-fPy;
00770 G4double Pz=theMomentum.pz()+iPz-fPz;
00771 G4double mE=theMomentum.e();
00773 theMomentum= G4LorentzVector(Px,Py,Pz,mE+corE);
00774 #ifdef debug
00775 G4double difF=fr4Mom.m2()-fM2;
00776 G4cout<<"G4QH::CorEDecIn2: dF="<<difF<<",out:"<<theMomentum<<fr4Mom<<G4endl;
00777 #endif
00778 return true;
00779 }
00780
00781
00782 G4bool G4QHadron::DecayIn3
00783 (G4LorentzVector& f4Mom, G4LorentzVector& s4Mom, G4LorentzVector& t4Mom)
00784 {
00785 #ifdef debug
00786 G4cout<<"G4QH::DIn3:"<<theMomentum<<"=>pf="<<f4Mom<<"+ps="<<s4Mom<<"+pt="<<t4Mom<<G4endl;
00787 #endif
00788 G4double iM = theMomentum.m();
00789 G4double fM = f4Mom.m();
00790 G4double sM = s4Mom.m();
00791 G4double tM = t4Mom.m();
00792 G4double eps = 0.001;
00793 if (fabs(iM-fM-sM-tM)<=eps)
00794 {
00795 G4double fR=fM/iM;
00796 G4double sR=sM/iM;
00797 G4double tR=tM/iM;
00798 f4Mom=fR*theMomentum;
00799 s4Mom=sR*theMomentum;
00800 t4Mom=tR*theMomentum;
00801 return true;
00802 }
00803 if (iM+eps<fM+sM+tM)
00804 {
00805 G4cout<<"***G4QHadron::DecayIn3:fM="<<fM<<" + sM="<<sM<<" + tM="<<tM<<" > iM="<<iM
00806 <<",d="<<iM-fM-sM-tM<<G4endl;
00807 return false;
00808 }
00809 G4double fM2 = fM*fM;
00810 G4double sM2 = sM*sM;
00811 G4double tM2 = tM*tM;
00812 G4double iM2 = iM*iM;
00813 G4double m13sBase=(iM-sM)*(iM-sM)-(fM+tM)*(fM+tM);
00814 G4double m12sMin =(fM+sM)*(fM+sM);
00815 G4double m12sBase=(iM-tM)*(iM-tM)-m12sMin;
00816 G4double rR = 0.;
00817 G4double rnd= 1.;
00818 #ifdef debug
00819 G4int tr = 0;
00820 #endif
00821 G4double m12s = 0.;
00822 while (rnd > rR)
00823 {
00824 m12s = m12sMin + m12sBase*G4UniformRand();
00825 G4double e1=m12s+fM2-sM2;
00826 G4double e2=iM2-m12s-tM2;
00827 G4double four12=4*m12s;
00828 G4double m13sRange=0.;
00829 G4double dif=(e1*e1-four12*fM2)*(e2*e2-four12*tM2);
00830 if(dif<0.)
00831 {
00832 #ifdef debug
00833 if(dif<-.01) G4cerr<<"*G4QHadron::DecayIn3:iM="<<iM<<",tM="<<tM<<",sM="<<sM<<",fM="
00834 <<fM<<",m12(s+f)="<<sqrt(m12s)<<", d="<<iM-fM-sM-tM<<G4endl;
00835 #endif
00836 }
00837 else m13sRange=sqrt(dif)/m12s;
00838 rR = m13sRange/m13sBase;
00839 rnd= G4UniformRand();
00840 #ifdef debug
00841 G4cout<<"G4QHadron::DecayIn3: try to decay #"<<++tr<<", rR="<<rR<<",rnd="<<rnd<<G4endl;
00842 #endif
00843 }
00844 G4double m12 = sqrt(m12s);
00845 G4LorentzVector dh4Mom(0.,0.,0.,m12);
00846
00847 if(!DecayIn2(t4Mom,dh4Mom))
00848 {
00849 G4cerr<<"***G4QHadron::DecayIn3: Exception1"<<G4endl;
00850
00851 return false;
00852 }
00853 #ifdef debug
00854 G4cout<<"G4QHadron::DecayIn3: Now the last decay of m12="<<dh4Mom.m()<<G4endl;
00855 #endif
00856 if(!G4QHadron(dh4Mom).DecayIn2(f4Mom,s4Mom))
00857 {
00858 G4cerr<<"***G4QHadron::DecayIn3: Error in DecayIn2 -> Exception2"<<G4endl;
00859
00860 return false;
00861 }
00862 return true;
00863 }
00864
00865
00866 G4bool G4QHadron::RelDecayIn3(G4LorentzVector& f4Mom, G4LorentzVector& s4Mom,
00867 G4LorentzVector& t4Mom, G4LorentzVector& dir,
00868 G4double maxCost, G4double minCost)
00869 {
00870 #ifdef debug
00871 G4cout<<"G4QH::RelDIn3:"<<theMomentum<<"=>f="<<f4Mom<<"+s="<<s4Mom<<"+t="<<t4Mom<<G4endl;
00872 #endif
00873 G4double iM = theMomentum.m();
00874 G4double fM = f4Mom.m();
00875 G4double sM = s4Mom.m();
00876 G4double tM = t4Mom.m();
00877 G4double eps = 0.001;
00878 if (fabs(iM-fM-sM-tM)<=eps)
00879 {
00880 G4double fR=fM/iM;
00881 G4double sR=sM/iM;
00882 G4double tR=tM/iM;
00883 f4Mom=fR*theMomentum;
00884 s4Mom=sR*theMomentum;
00885 t4Mom=tR*theMomentum;
00886 return true;
00887 }
00888 if (iM+eps<fM+sM+tM)
00889 {
00890 G4cout<<"***G4QHadron::RelDecayIn3:fM="<<fM<<" + sM="<<sM<<" + tM="<<tM<<" > iM="<<iM
00891 <<",d="<<iM-fM-sM-tM<<G4endl;
00892 return false;
00893 }
00894 G4double fM2 = fM*fM;
00895 G4double sM2 = sM*sM;
00896 G4double tM2 = tM*tM;
00897 G4double iM2 = iM*iM;
00898 G4double m13sBase=(iM-sM)*(iM-sM)-(fM+tM)*(fM+tM);
00899 G4double m12sMin =(fM+sM)*(fM+sM);
00900 G4double m12sBase=(iM-tM)*(iM-tM)-m12sMin;
00901 G4double rR = 0.;
00902 G4double rnd= 1.;
00903 #ifdef debug
00904 G4int tr = 0;
00905 #endif
00906 G4double m12s = 0.;
00907 while (rnd > rR)
00908 {
00909 m12s = m12sMin + m12sBase*G4UniformRand();
00910 G4double e1=m12s+fM2-sM2;
00911 G4double e2=iM2-m12s-tM2;
00912 G4double four12=4*m12s;
00913 G4double m13sRange=0.;
00914 G4double dif=(e1*e1-four12*fM2)*(e2*e2-four12*tM2);
00915 if(dif<0.)
00916 {
00917 #ifdef debug
00918 if(dif<-.01) G4cerr<<"G4QHadron::RelDecayIn3:iM="<<iM<<",tM="<<tM<<",sM="<<sM<<",fM="
00919 <<fM<<",m12(s+f)="<<sqrt(m12s)<<", d="<<iM-fM-sM-tM<<G4endl;
00920 #endif
00921 }
00922 else m13sRange=sqrt(dif)/m12s;
00923 rR = m13sRange/m13sBase;
00924 rnd= G4UniformRand();
00925 #ifdef debug
00926 G4cout<<"G4QHadron::RelDecayIn3: try decay #"<<++tr<<", rR="<<rR<<",rnd="<<rnd<<G4endl;
00927 #endif
00928 }
00929 G4double m12 = sqrt(m12s);
00930 G4LorentzVector dh4Mom(0.,0.,0.,m12);
00931
00932 if(!RelDecayIn2(t4Mom,dh4Mom,dir,maxCost,minCost))
00933 {
00934 G4cerr<<"***G4QHadron::RelDecayIn3: Exception1"<<G4endl;
00935
00936 return false;
00937 }
00938 #ifdef debug
00939 G4cout<<"G4QHadron::RelDecayIn3: Now the last decay of m12="<<dh4Mom.m()<<G4endl;
00940 #endif
00941 if(!G4QHadron(dh4Mom).DecayIn2(f4Mom,s4Mom))
00942 {
00943 G4cerr<<"***G4QHadron::RelDecayIn3: Error in DecayIn2 -> Exception2"<<G4endl;
00944
00945 return false;
00946 }
00947 return true;
00948 }
00949
00950
00951 G4bool G4QHadron::CopDecayIn3(G4LorentzVector& f4Mom, G4LorentzVector& s4Mom,
00952 G4LorentzVector& t4Mom, G4LorentzVector& dir, G4double cosp)
00953 {
00954 #ifdef debug
00955 G4cout<<"G4QH::CopDIn3:"<<theMomentum<<"=>f="<<f4Mom<<"+s="<<s4Mom<<"+t="<<t4Mom<<G4endl;
00956 #endif
00957 G4double iM = theMomentum.m();
00958 G4double fM = f4Mom.m();
00959 G4double sM = s4Mom.m();
00960 G4double tM = t4Mom.m();
00961 G4double eps = 0.001;
00962 if (fabs(iM-fM-sM-tM)<=eps)
00963 {
00964 G4double fR=fM/iM;
00965 G4double sR=sM/iM;
00966 G4double tR=tM/iM;
00967 f4Mom=fR*theMomentum;
00968 s4Mom=sR*theMomentum;
00969 t4Mom=tR*theMomentum;
00970 return true;
00971 }
00972 if (iM+eps<fM+sM+tM)
00973 {
00974 G4cout<<"***G4QHadron::CopDecayIn3:fM="<<fM<<" + sM="<<sM<<" + tM="<<tM<<" > iM="<<iM
00975 <<",d="<<iM-fM-sM-tM<<G4endl;
00976 return false;
00977 }
00978 G4double fM2 = fM*fM;
00979 G4double sM2 = sM*sM;
00980 G4double tM2 = tM*tM;
00981 G4double iM2 = iM*iM;
00982 G4double m13sBase=(iM-sM)*(iM-sM)-(fM+tM)*(fM+tM);
00983 G4double m12sMin =(fM+sM)*(fM+sM);
00984 G4double m12sBase=(iM-tM)*(iM-tM)-m12sMin;
00985 G4double rR = 0.;
00986 G4double rnd= 1.;
00987 #ifdef debug
00988 G4int tr = 0;
00989 #endif
00990 G4double m12s = 0.;
00991 while (rnd > rR)
00992 {
00993 m12s = m12sMin + m12sBase*G4UniformRand();
00994 G4double e1=m12s+fM2-sM2;
00995 G4double e2=iM2-m12s-tM2;
00996 G4double four12=4*m12s;
00997 G4double m13sRange=0.;
00998 G4double dif=(e1*e1-four12*fM2)*(e2*e2-four12*tM2);
00999 if(dif<0.)
01000 {
01001 #ifdef debug
01002 if(dif<-.01) G4cerr<<"G4QHadron::CopDecayIn3:iM="<<iM<<",tM="<<tM<<",sM="<<sM<<",fM="
01003 <<fM<<",m12(s+f)="<<sqrt(m12s)<<", d="<<iM-fM-sM-tM<<G4endl;
01004 #endif
01005 }
01006 else m13sRange=sqrt(dif)/m12s;
01007 rR = m13sRange/m13sBase;
01008 rnd= G4UniformRand();
01009 #ifdef debug
01010 G4cout<<"G4QHadron::CopDecayIn3: try decay #"<<++tr<<", rR="<<rR<<",rnd="<<rnd<<G4endl;
01011 #endif
01012 }
01013 G4double m12 = sqrt(m12s);
01014 G4LorentzVector dh4Mom(0.,0.,0.,m12);
01015
01016 if(!CopDecayIn2(t4Mom,dh4Mom,dir,cosp))
01017 {
01018 G4cerr<<"***G4QHadron::CopDecayIn3: Exception1"<<G4endl;
01019
01020 return false;
01021 }
01022 #ifdef debug
01023 G4cout<<"G4QHadron::DecayIn3: Now the last decay of m12="<<dh4Mom.m()<<G4endl;
01024 #endif
01025 if(!G4QHadron(dh4Mom).DecayIn2(f4Mom,s4Mom))
01026 {
01027 G4cerr<<"***G4QHadron::CopDecayIn3: Error in DecayIn2 -> Exception2"<<G4endl;
01028
01029 return false;
01030 }
01031 return true;
01032 }
01033
01034
01035 G4double G4QHadron::RandomizeMass(G4QParticle* pPart, G4double maxM)
01036 {
01037 G4double meanM = theQPDG.GetMass();
01038 G4double width = theQPDG.GetWidth()/2.;
01039 #ifdef debug
01040 G4cout<<"G4QHadron::RandomizeMass: meanM="<<meanM<<", halfWidth="<<width<<G4endl;
01041 #endif
01042 if(maxM<meanM-3*width)
01043 {
01044 #ifdef debug
01045 G4cout<<"***G4QH::RandM:m=0 maxM="<<maxM<<"<meanM="<<meanM<<"-3*halfW="<<width<<G4endl;
01046 #endif
01047 return 0.;
01048 }
01050 if(width==0.)
01051 {
01052 #ifdef debug
01053 if(meanM>maxM) G4cerr<<"***G4QHadron::RandM:Stable m="<<meanM<<">maxM="<<maxM<<G4endl;
01054 #endif
01055 return meanM;
01056
01057 }
01058 else if(width<0.)
01059 {
01060
01061
01062 G4ExceptionDescription ed;
01063 ed << "width of the Hadron < 0. : width=" << width << "<0,PDGC="
01064 << theQPDG.GetPDGCode() << G4endl;
01065 G4Exception("G4QHadron::RandomizeMass()", "HAD_CHPS_0000", FatalException, ed);
01066 }
01067 G4double minM = pPart->MinMassOfFragm();
01068 if(minM>maxM)
01069 {
01070 #ifdef debug
01071 G4cout<<"***G4QHadron::RandomizeMass:for PDG="<<theQPDG.GetPDGCode()<<" minM="<<minM
01072 <<" > maxM="<<maxM<<G4endl;
01073 #endif
01074 return 0.;
01075 }
01076
01077 G4double v1=atan((minM-meanM)/width);
01078 G4double v2=atan((maxM-meanM)/width);
01079 G4double dv=v2-v1;
01080 #ifdef debug
01081 G4cout<<"G4QHadr::RandM:Mi="<<minM<<",i="<<v1<<",Ma="<<maxM<<",a="<<v2<<","<<dv<<G4endl;
01082 #endif
01083 return meanM+width*tan(v1+dv*G4UniformRand());
01084 }
01085
01086
01087 void G4QHadron::SplitUp()
01088 {
01089 if (isSplit) return;
01090 #ifdef pdebug
01091 G4cout<<"G4QHadron::SplitUp ***IsCalled***, before Splitting nC="<<Color.size()
01092 <<", SoftColCount="<<theCollisionCount<<G4endl;
01093 #endif
01094 isSplit=true;
01095 if (!Color.empty()) return;
01096 if (!theCollisionCount)
01097 {
01098 G4QParton* Left = 0;
01099 G4QParton* Right = 0;
01100 GetValenceQuarkFlavors(Left, Right);
01101 G4ThreeVector Pos=GetPosition();
01102 Left->SetPosition(Pos);
01103 Right->SetPosition(Pos);
01104
01105 G4double theMomPlus = theMomentum.plus();
01106 G4double theMomMinus = theMomentum.minus();
01107 #ifdef pdebug
01108 G4cout<<"G4QHadron::SplitUp: *Dif* possition="<<Pos<<", 4M="<<theMomentum<<G4endl;
01109 #endif
01110
01111
01112 G4double pt2 = theMomentum.perp2();
01113 G4double transverseMass2 = theMomPlus*theMomMinus;
01114 if(transverseMass2<0.) transverseMass2=0.;
01115 G4double maxAvailMomentum2 = sqr(std::sqrt(transverseMass2) - std::sqrt(pt2));
01116 G4ThreeVector pt(0., 0., 0.);
01117 if(maxAvailMomentum2/widthOfPtSquare > 0.01)
01118 pt=GaussianPt(widthOfPtSquare, maxAvailMomentum2);
01119 #ifdef pdebug
01120 G4cout<<"G4QHadron::SplitUp: *Dif* maxMom2="<<maxAvailMomentum2<<", pt="<<pt<<G4endl;
01121 #endif
01122 G4LorentzVector LeftMom(pt, 0.);
01123 G4LorentzVector RightMom;
01124 RightMom.setPx(theMomentum.px() - pt.x());
01125 RightMom.setPy(theMomentum.py() - pt.y());
01126 #ifdef pdebug
01127 G4cout<<"G4QHadron::SplitUp: *Dif* right4m="<<RightMom<<", left4M="<<LeftMom<<G4endl;
01128 #endif
01129 G4double Local1 = theMomMinus + (RightMom.perp2() - LeftMom.perp2()) / theMomPlus;
01130 G4double Local2 = std::sqrt(std::max(0., Local1*Local1 -
01131 4*RightMom.perp2()*theMomMinus / theMomPlus));
01132 #ifdef pdebug
01133 G4cout<<"G4QHadron::SplitUp:Dif,L1="<<Local1<<",L2="<<Local2<<",D="<<Direction<<G4endl;
01134 #endif
01135 if (Direction) Local2 = -Local2;
01136 G4double RightMinus = 0.5*(Local1 + Local2);
01137 G4double LeftMinus = theMomentum.minus() - RightMinus;
01138 #ifdef pdebug
01139 G4cout<<"G4QHadron::SplitUp: *Dif* Rminus="<<RightMinus<<",Lminus="<<LeftMinus<<",hmm="
01140 <<theMomentum.minus()<<G4endl;
01141 #endif
01142 G4double LeftPlus = 0.;
01143 if(LeftMinus) LeftPlus = LeftMom.perp2()/LeftMinus;
01144 G4double RightPlus = theMomentum.plus() - LeftPlus;
01145 #ifdef pdebug
01146 G4cout<<"G4QHadron::SplitUp: *Dif* Rplus="<<RightPlus<<", Lplus="<<LeftPlus<<G4endl;
01147 #endif
01148 LeftMom.setPz(0.5*(LeftPlus - LeftMinus));
01149 LeftMom.setE (0.5*(LeftPlus + LeftMinus));
01150 RightMom.setPz(0.5*(RightPlus - RightMinus));
01151 RightMom.setE (0.5*(RightPlus + RightMinus));
01152
01153 #ifdef pdebug
01154 G4cout<<"G4QHadron::SplitUp: *Dif* -final- R4m="<<RightMom<<", L4M="<<LeftMom<<", L+R="
01155 <<RightMom+LeftMom<<", D4M="<<theMomentum-RightMom+LeftMom<<G4endl;
01156 #endif
01157 Left->Set4Momentum(LeftMom);
01158 Right->Set4Momentum(RightMom);
01159 Color.push_back(Left);
01160 AntiColor.push_back(Right);
01161 }
01162 else
01163 {
01164
01165
01166 G4ThreeVector SumP(0.,0.,0.);
01167 G4ThreeVector Pos = GetPosition();
01168 G4int nSeaPair = theCollisionCount-1;
01169 #ifdef pdebug
01170 G4cout<<"G4QHadron::SplitUp:*Soft* Pos="<<Pos<<", nSeaPair="<<nSeaPair<<G4endl;
01171 #endif
01172
01173
01174 for (G4int aSeaPair = 0; aSeaPair < nSeaPair; aSeaPair++)
01175 {
01176
01177 G4int aPDGCode = 1 + (G4int)(G4UniformRand()/StrangeSuppress);
01178
01179
01180
01181 G4QParton* aParton = BuildSeaQuark(false, aPDGCode);
01182
01183
01184 G4int firstPartonColour = aParton->GetColour();
01185 G4double firstPartonSpinZ = aParton->GetSpinZ();
01186 #ifdef pdebug
01187 G4cout<<"G4QHadron::SplitUp:*Soft* Part1 PDG="<<aPDGCode<<", Col="<<firstPartonColour
01188 <<", SpinZ="<<firstPartonSpinZ<<", 4M="<<aParton->Get4Momentum()<<G4endl;
01189 #endif
01190 SumP+=aParton->Get4Momentum();
01191 Color.push_back(aParton);
01192
01193
01194 aParton = BuildSeaQuark(true, aPDGCode);
01195 aParton->SetSpinZ(-firstPartonSpinZ);
01196 aParton->SetColour(-firstPartonColour);
01197 #ifdef pdebug
01198 G4cout<<"G4QHadron::SplUp:Sft,P2="<<aParton->Get4Momentum()<<",i="<<aSeaPair<<G4endl;
01199 #endif
01200
01201 SumP+=aParton->Get4Momentum();
01202 AntiColor.push_back(aParton);
01203 #ifdef pdebug
01204 G4cout<<"G4QHadron::SplUp:*Sft* Antiquark is filled, i="<<aSeaPair<<G4endl;
01205 #endif
01206 }
01207
01208 G4QParton* pColorParton = 0;
01209 G4QParton* pAntiColorParton = 0;
01210 GetValenceQuarkFlavors(pColorParton, pAntiColorParton);
01211 G4int ColorEncoding = pColorParton->GetPDGCode();
01212 #ifdef pdebug
01213 G4int AntiColorEncoding = pAntiColorParton->GetPDGCode();
01214 G4cout<<"G4QHadron::SplUp:*Sft*,C="<<ColorEncoding<<", AC="<<AntiColorEncoding<<G4endl;
01215 #endif
01216 G4ThreeVector ptr = GaussianPt(sigmaPt, DBL_MAX);
01217 SumP += ptr;
01218 #ifdef pdebug
01219 G4cout<<"G4QHadron::SplitUp: *Sft*, ptr="<<ptr<<G4endl;
01220 #endif
01221
01222 if (ColorEncoding < 0)
01223 {
01224 G4LorentzVector ColorMom(-SumP, 0);
01225 pColorParton->Set4Momentum(ColorMom);
01226 G4LorentzVector AntiColorMom(ptr, 0.);
01227 pAntiColorParton->Set4Momentum(AntiColorMom);
01228 }
01229 else
01230 {
01231 G4LorentzVector ColorMom(ptr, 0);
01232 pColorParton->Set4Momentum(ColorMom);
01233 G4LorentzVector AntiColorMom(-SumP, 0);
01234 pAntiColorParton->Set4Momentum(AntiColorMom);
01235 }
01236 Color.push_back(pColorParton);
01237 AntiColor.push_back(pAntiColorParton);
01238 #ifdef pdebug
01239 G4cout<<"G4QHadron::SplitUp: *Soft* Col&Anticol are filled PDG="<<GetPDGCode()<<G4endl;
01240 #endif
01241
01242 G4int nColor=Color.size();
01243 G4int nAntiColor=AntiColor.size();
01244 if(nColor!=nAntiColor || nColor != nSeaPair+1)
01245 {
01246 G4cerr<<"***G4QHadron::SplitUp: nA="<<nAntiColor<<",nAC="<<nColor<<",nSea="<<nSeaPair
01247 <<G4endl;
01248 G4Exception("G4QHadron::SplitUp:","72",FatalException,"Colours&AntiColours notSinc");
01249 }
01250 #ifdef pdebug
01251 G4cout<<"G4QHad::SpUp:,nPartons="<<nColor+nColor<<G4endl;
01252 #endif
01253 G4int dnCol=nColor+nColor;
01254
01255 G4double* xs=RandomX(dnCol);
01256
01257
01258 #ifdef pdebug
01259 G4cout<<"G4QHadron::SplitUp:*Sft* Loop ColorX="<<ColorX<<G4endl;
01260 #endif
01261 std::list<G4QParton*>::iterator icolor = Color.begin();
01262 std::list<G4QParton*>::iterator ecolor = Color.end();
01263 std::list<G4QParton*>::iterator ianticolor = AntiColor.begin();
01264 std::list<G4QParton*>::iterator eanticolor = AntiColor.end();
01265 G4int xi=-1;
01266
01267 for ( ; icolor != ecolor && ianticolor != eanticolor; ++icolor, ++ianticolor)
01268 {
01269 (*icolor)->SetX(xs[++xi]);
01270
01271
01272
01273
01274
01275 (*icolor)->DefineEPz(theMomentum);
01276 (*ianticolor)->SetX(xs[++xi]);
01277
01278
01279
01280
01281
01282 (*ianticolor)->DefineEPz(theMomentum);
01283 }
01284 delete[] xs;
01285 #ifdef pdebug
01286 G4cout<<"G4QHadron::SplitUp: *Soft* ===> End, ColSize="<<Color.size()<<G4endl;
01287 #endif
01288 return;
01289 }
01290 }
01291
01292
01293 void G4QHadron::Boost(const G4LorentzVector& boost4M)
01294 {
01295
01296 G4double bm=boost4M.mag();
01297 G4double factor=(theMomentum.vect()*boost4M.vect()/(boost4M.e()+bm)-theMomentum.e())/bm;
01298 theMomentum.setE(theMomentum.dot(boost4M)/bm);
01299 theMomentum.setVect(factor*boost4M.vect() + theMomentum.vect());
01300 }
01301
01302
01303 G4QParton* G4QHadron::BuildSeaQuark(G4bool isAntiQuark, G4int aPDGCode)
01304 {
01305 if (isAntiQuark) aPDGCode*=-1;
01306 G4QParton* result = new G4QParton(aPDGCode);
01307 result->SetPosition(GetPosition());
01308 G4ThreeVector aPtVector = GaussianPt(sigmaPt, DBL_MAX);
01309 G4LorentzVector a4Momentum(aPtVector, 0);
01310 result->Set4Momentum(a4Momentum);
01311 return result;
01312 }
01313
01314
01315 G4double* G4QHadron::RandomX(G4int nPart)
01316 {
01317 G4double* x = 0;
01318 if(nPart<2)
01319 {
01320 G4cout<<"-Warning-G4QHadron::RandomX: nPart="<<nPart<<" < 2"<<G4endl;
01321 return x;
01322 }
01323 x = new G4double[nPart];
01324 G4int nP1=nPart-1;
01325 x[0]=G4UniformRand();
01326 for(G4int i=1; i<nP1; ++i)
01327 {
01328 G4double r=G4UniformRand();
01329 G4int j=0;
01330 for( ; j<i; ++j) if(r < x[j])
01331 {
01332 for(G4int k=i; k>j; --k) x[k]=x[k-1];
01333 x[j]=r;
01334 break;
01335 }
01336 if(j==i) x[i]=r;
01337 }
01338 x[nP1]=1.;
01339 for(G4int i=nP1; i>0; --i) x[i]-=x[i-1];
01340 return x;
01341 }
01342
01343
01344 G4double G4QHadron::SampleCHIPSX(G4double anXtot, G4int nSea)
01345 {
01346 G4double ns_value=nSea;
01347 if(nSea<1 || anXtot<=0.) G4cout<<"-Warning-G4QHad::SCX:N="<<nSea<<",tX="<<anXtot<<G4endl;
01348 if(nSea<2) return anXtot;
01349 return anXtot*(1.-std::pow(G4UniformRand(),1./ns_value));
01350 }
01351
01352
01353 void G4QHadron::GetValenceQuarkFlavors(G4QParton* &Parton1, G4QParton* &Parton2)
01354 {
01355
01356 G4int aEnd=0;
01357 G4int bEnd=0;
01358 G4int HadronEncoding = GetPDGCode();
01359 if(!(GetBaryonNumber())) SplitMeson(HadronEncoding, &aEnd, &bEnd);
01360 else SplitBaryon(HadronEncoding, &aEnd, &bEnd);
01361
01362 Parton1 = new G4QParton(aEnd);
01363 Parton1->SetPosition(GetPosition());
01364
01365
01366
01367
01368
01369
01370
01371
01372
01373 Parton2 = new G4QParton(bEnd);
01374 Parton2->SetPosition(GetPosition());
01375
01376
01377
01378
01379
01380
01381
01382
01383
01384
01385
01386
01387 Parton2->SetColour(-(Parton1->GetColour()));
01388
01389
01390
01391
01392
01393
01394 if ( std::abs(Parton1->GetSpinZ() + Parton2->GetSpinZ()) > GetSpin())
01395 Parton2->SetSpinZ(-(Parton2->GetSpinZ()));
01396
01397
01398
01399
01400
01401
01402
01403
01404
01405 }
01406
01407 G4bool G4QHadron::SplitMeson(G4int PDGcode, G4int* aEnd, G4int* bEnd)
01408 {
01409 G4bool result = true;
01410 G4int absPDGcode = std::abs(PDGcode);
01411 if (absPDGcode >= 1000) return false;
01412 if(absPDGcode == 22 || absPDGcode == 111)
01413 {
01414 G4int it=1;
01415 if(G4UniformRand()<.5) it++;
01416 *aEnd = it;
01417 *bEnd = -it;
01418 }
01419 else if(absPDGcode == 130 || absPDGcode == 310)
01420 {
01421 G4int it=1;
01422 if(G4UniformRand()<.5) it=-1;
01423 *aEnd = it;
01424 if(it>0) *bEnd = -3;
01425 else *bEnd = 3;
01426 }
01427 else
01428 {
01429 G4int heavy = absPDGcode/100;
01430 G4int light = (absPDGcode%100)/10;
01431 G4int anti = 1 - 2*(std::max(heavy, light)%2);
01432 if (PDGcode < 0 ) anti = -anti;
01433 heavy *= anti;
01434 light *= -anti;
01435 if ( anti < 0) G4SwapObj(&heavy, &light);
01436 *aEnd = heavy;
01437 *bEnd = light;
01438 }
01439 return result;
01440 }
01441
01442 G4bool G4QHadron::SplitBaryon(G4int PDGcode, G4int* quark, G4int* diQuark)
01443 {
01444 static const G4double r2=.5;
01445 static const G4double r3=1./3.;
01446 static const G4double d3=2./3.;
01447 static const G4double r4=1./4.;
01448 static const G4double r6=1./6.;
01449 static const G4double r12=1./12.;
01450
01451 std::pair<G4int,G4int> qdq[5];
01452 G4double prb[5];
01453 G4int nc=0;
01454 G4int aPDGcode=std::abs(PDGcode);
01455 if(aPDGcode==2212)
01456 {
01457 nc=3;
01458 qdq[0]=make_pair(2203, 1); prb[0]=r3;
01459 qdq[1]=make_pair(2103, 2); prb[1]=r6;
01460 qdq[2]=make_pair(2101, 2); prb[2]=r2;
01461 }
01462 else if(aPDGcode==2112)
01463 {
01464 nc=3;
01465 qdq[0]=make_pair(2103, 1); prb[0]=r6;
01466 qdq[1]=make_pair(2101, 1); prb[1]=r2;
01467 qdq[2]=make_pair(1103, 2); prb[2]=r3;
01468 }
01469 else if(aPDGcode%10<3)
01470 {
01471 if(aPDGcode==3122)
01472 {
01473 nc=5;
01474 qdq[0]=make_pair(2103, 3); prb[0]=r3;
01475 qdq[1]=make_pair(3203, 1); prb[1]=r4;
01476 qdq[2]=make_pair(3201, 1); prb[2]=r12;
01477 qdq[3]=make_pair(3103, 2); prb[3]=r4;
01478 qdq[4]=make_pair(3101, 2); prb[4]=r12;
01479 }
01480 else if(aPDGcode==3222)
01481 {
01482 nc=3;
01483 qdq[0]=make_pair(2203, 3); prb[0]=r3;
01484 qdq[1]=make_pair(3203, 2); prb[1]=r6;
01485 qdq[2]=make_pair(3201, 2); prb[2]=r2;
01486 }
01487 else if(aPDGcode==3212)
01488 {
01489 nc=5;
01490 qdq[0]=make_pair(2103, 3); prb[0]=r3;
01491 qdq[1]=make_pair(3203, 1); prb[1]=r12;
01492 qdq[2]=make_pair(3201, 1); prb[2]=r4;
01493 qdq[3]=make_pair(3103, 2); prb[3]=r12;
01494 qdq[4]=make_pair(3101, 2); prb[4]=r4;
01495 }
01496 else if(aPDGcode==3112)
01497 {
01498 nc=3;
01499 qdq[0]=make_pair(1103, 3); prb[0]=r3;
01500 qdq[1]=make_pair(3103, 1); prb[1]=r6;
01501 qdq[2]=make_pair(3101, 1); prb[2]=r2;
01502 }
01503 else if(aPDGcode==3312)
01504 {
01505 nc=3;
01506 qdq[0]=make_pair(3103, 3); prb[0]=r6;
01507 qdq[1]=make_pair(3101, 3); prb[1]=r2;
01508 qdq[2]=make_pair(3303, 1); prb[2]=r3;
01509 }
01510 else if(aPDGcode==3322)
01511 {
01512 nc=3;
01513 qdq[0]=make_pair(3203, 3); prb[0]=r6;
01514 qdq[1]=make_pair(3201, 3); prb[1]=r2;
01515 qdq[2]=make_pair(3303, 2); prb[2]=r3;
01516 }
01517 else return false;
01518 }
01519 else
01520 {
01521 if(aPDGcode==3334)
01522 {
01523 nc=1;
01524 qdq[0]=make_pair(3303, 3); prb[0]=1.;
01525 }
01526 else if(aPDGcode==2224)
01527 {
01528 nc=1;
01529 qdq[0]=make_pair(2203, 2); prb[0]=1.;
01530 }
01531 else if(aPDGcode==2214)
01532 {
01533 nc=2;
01534 qdq[0]=make_pair(2203, 1); prb[0]=r3;
01535 qdq[1]=make_pair(2103, 2); prb[1]=d3;
01536 }
01537 else if(aPDGcode==2114)
01538 {
01539 nc=2;
01540 qdq[0]=make_pair(1103, 2); prb[0]=r3;
01541 qdq[1]=make_pair(2103, 1); prb[1]=d3;
01542 }
01543 else if(aPDGcode==1114)
01544 {
01545 nc=1;
01546 qdq[0]=make_pair(1103, 1); prb[0]=1.;
01547 }
01548 else if(aPDGcode==3224)
01549 {
01550 nc=2;
01551 qdq[0]=make_pair(2203, 3); prb[0]=r3;
01552 qdq[1]=make_pair(3203, 2); prb[1]=d3;
01553 }
01554 else if(aPDGcode==3214)
01555 {
01556 nc=3;
01557 qdq[0]=make_pair(2103, 3); prb[0]=r3;
01558 qdq[1]=make_pair(3203, 1); prb[1]=r3;
01559 qdq[2]=make_pair(3103, 2); prb[2]=r3;
01560 }
01561 else if(aPDGcode==3114)
01562 {
01563 nc=2;
01564 qdq[0]=make_pair(1103, 3); prb[0]=r3;
01565 qdq[1]=make_pair(3103, 1); prb[1]=d3;
01566 }
01567 else if(aPDGcode==3324)
01568 {
01569 nc=2;
01570 qdq[0]=make_pair(3203, 3); prb[0]=r3;
01571 qdq[1]=make_pair(3303, 2); prb[1]=d3;
01572 }
01573 else if(aPDGcode==3314)
01574 {
01575 nc=2;
01576 qdq[0]=make_pair(3103, 3); prb[0]=d3;
01577 qdq[1]=make_pair(3303, 1); prb[1]=r3;
01578 }
01579 else return false;
01580 }
01581 G4double random = G4UniformRand();
01582 G4double sum = 0.;
01583 for(G4int i=0; i<nc; i++)
01584 {
01585 sum += prb[i];
01586 if(sum>random)
01587 {
01588 if(PDGcode>0)
01589 {
01590 *diQuark= qdq[i].first;
01591 *quark = qdq[i].second;
01592 }
01593 else
01594 {
01595 *diQuark= -qdq[i].second;
01596 *quark = -qdq[i].first;
01597 }
01598 break;
01599 }
01600 }
01601 return true;
01602 }
01603
01604
01605 G4ThreeVector G4QHadron::GaussianPt(G4double widthSquare, G4double maxPtSquare)
01606 {
01607 G4double R=0.;
01608 while ((R = -widthSquare*std::log(G4UniformRand())) > maxPtSquare){}
01609 R = std::sqrt(R);
01610 G4double phi = twopi*G4UniformRand();
01611 return G4ThreeVector(R*std::cos(phi), R*std::sin(phi), 0.);
01612 }
01613
01614 G4QParton* G4QHadron::GetNextParton()
01615 {
01616 if(Color.size()==0) return 0;
01617 G4QParton* result = Color.back();
01618 Color.pop_back();
01619 return result;
01620 }
01621
01622 G4QParton* G4QHadron::GetNextAntiParton()
01623 {
01624 if(AntiColor.size() == 0) return 0;
01625 G4QParton* result = AntiColor.front();
01626 AntiColor.pop_front();
01627 return result;
01628 }
01629
01630
01631 G4QPartonPair* G4QHadron::SplitInTwoPartons()
01632 {
01633 if(std::abs(GetBaryonNumber())>1)
01634 {
01635 G4cerr<<"***G4QHadron::SplitInTwoPartons: Can not split QC="<<valQ<< G4endl;
01636 G4Exception("G4QFragmentation::ChooseX:","72",FatalException,"NeitherMesonNorBaryon");
01637 }
01638 std::pair<G4int,G4int> PP = valQ.MakePartonPair();
01639 return new G4QPartonPair(new G4QParton(PP.first), new G4QParton(PP.second));
01640 }