[819] | 1 | // |
---|
| 2 | // ******************************************************************** |
---|
| 3 | // * License and Disclaimer * |
---|
| 4 | // * * |
---|
| 5 | // * The Geant4 software is copyright of the Copyright Holders of * |
---|
| 6 | // * the Geant4 Collaboration. It is provided under the terms and * |
---|
| 7 | // * conditions of the Geant4 Software License, included in the file * |
---|
| 8 | // * LICENSE and available at http://cern.ch/geant4/license . These * |
---|
| 9 | // * include a list of copyright holders. * |
---|
| 10 | // * * |
---|
| 11 | // * Neither the authors of this software system, nor their employing * |
---|
| 12 | // * institutes,nor the agencies providing financial support for this * |
---|
| 13 | // * work make any representation or warranty, express or implied, * |
---|
| 14 | // * regarding this software system or assume any liability for its * |
---|
| 15 | // * use. Please see the license in the file LICENSE and URL above * |
---|
| 16 | // * for the full disclaimer and the limitation of liability. * |
---|
| 17 | // * * |
---|
| 18 | // * This code implementation is the result of the scientific and * |
---|
| 19 | // * technical work of the GEANT4 collaboration. * |
---|
| 20 | // * By using, copying, modifying or distributing the software (or * |
---|
| 21 | // * any work based on the software) you agree to acknowledge its * |
---|
| 22 | // * use in resulting scientific publications, and indicate your * |
---|
| 23 | // * acceptance of all terms of the Geant4 Software license. * |
---|
| 24 | // ******************************************************************** |
---|
| 25 | // |
---|
| 26 | // |
---|
| 27 | // $Id: G4QHadron.cc,v 1.52.2.1 2008/04/23 14:47:23 gcosmo Exp $ |
---|
| 28 | // GEANT4 tag $Name: geant4-09-01-patch-02 $ |
---|
| 29 | // |
---|
| 30 | // ---------------- G4QHadron ---------------- |
---|
| 31 | // by Mikhail Kossov, Sept 1999. |
---|
| 32 | // class for Quasmon initiated Hadrons generated by CHIPS Model |
---|
| 33 | // ------------------------------------------------------------------ |
---|
| 34 | // |
---|
| 35 | //#define debug |
---|
| 36 | //#define pdebug |
---|
| 37 | //#define sdebug |
---|
| 38 | //#define ppdebug |
---|
| 39 | |
---|
| 40 | #include "G4QHadron.hh" |
---|
| 41 | #include <cmath> |
---|
| 42 | using namespace std; |
---|
| 43 | |
---|
| 44 | G4double G4QHadron::alpha = -0.5; // changing rapidity distribution for all |
---|
| 45 | G4double G4QHadron::beta = 2.5; // changing rapidity distribution for projectile region |
---|
| 46 | G4double G4QHadron::theMinPz = 70.*MeV; // Can be from 14 to 140 MeV |
---|
| 47 | G4double G4QHadron::StrangeSuppress = 0.48; // ? M.K. |
---|
| 48 | G4double G4QHadron::sigmaPt = 1.7*GeV; // Can be 0 |
---|
| 49 | G4double G4QHadron::widthOfPtSquare = 0.01*GeV*GeV; // ? M.K. |
---|
| 50 | G4double G4QHadron::minTransverseMass = 1.*keV; // ? M.K. |
---|
| 51 | |
---|
| 52 | G4QHadron::G4QHadron() : theQPDG(0),theMomentum(0.,0.,0.,0.),valQ(0,0,0,0,0,0),nFragm(0), |
---|
| 53 | thePosition(0.,0.,0.),theCollisionCount(0),isSplit(false),Direction(true), |
---|
| 54 | Color(),AntiColor(),bindE(0.),formTime(0.) {} |
---|
| 55 | |
---|
| 56 | G4QHadron::G4QHadron(G4LorentzVector p) : theQPDG(0),theMomentum(p),valQ(0,0,0,0,0,0), |
---|
| 57 | nFragm(0),thePosition(0.,0.,0.),theCollisionCount(0),isSplit(false),Direction(true), |
---|
| 58 | Color(),AntiColor(),bindE(0.),formTime(0.) {} |
---|
| 59 | |
---|
| 60 | // For Chipolino or Quasmon doesn't make any sense |
---|
| 61 | G4QHadron::G4QHadron(G4int PDGCode, G4LorentzVector p) : theQPDG(PDGCode),theMomentum(p), |
---|
| 62 | nFragm(0),thePosition(0.,0.,0.),theCollisionCount(0),isSplit(false),Direction(true), |
---|
| 63 | Color(),AntiColor(),bindE(0.),formTime(0.) |
---|
| 64 | { |
---|
| 65 | #ifdef debug |
---|
| 66 | G4cout<<"G4QHadron must be created with PDG="<<PDGCode<<", 4M="<<p<<G4endl; |
---|
| 67 | #endif |
---|
| 68 | if(GetQCode()>-1) |
---|
| 69 | { |
---|
| 70 | if(theMomentum.e()==0.) theMomentum.setE(theQPDG.GetMass()); |
---|
| 71 | valQ=theQPDG.GetQuarkContent(); |
---|
| 72 | } |
---|
| 73 | else if(PDGCode>80000000) DefineQC(PDGCode); |
---|
| 74 | else G4cerr<<"***G4QHadron:(P) PDG="<<PDGCode<<", use other constructor"<<G4endl; |
---|
| 75 | #ifdef debug |
---|
| 76 | G4cout<<"G4QHadron is created with QCode="<<GetQCode()<<", QC="<<valQ<<G4endl; |
---|
| 77 | #endif |
---|
| 78 | } |
---|
| 79 | |
---|
| 80 | // For Chipolino or Quasmon doesn't make any sense |
---|
| 81 | G4QHadron::G4QHadron(G4QPDGCode QPDG, G4LorentzVector p) : theQPDG(QPDG),theMomentum(p), |
---|
| 82 | nFragm(0),thePosition(0.,0.,0.),theCollisionCount(0),isSplit(false),Direction(true), |
---|
| 83 | Color(),AntiColor(),bindE(0.),formTime(0.) |
---|
| 84 | { |
---|
| 85 | if(theQPDG.GetQCode()>-1) |
---|
| 86 | { |
---|
| 87 | if(theMomentum.e()==0.) theMomentum.setE(theQPDG.GetMass()); |
---|
| 88 | valQ=theQPDG.GetQuarkContent(); |
---|
| 89 | } |
---|
| 90 | else |
---|
| 91 | { |
---|
| 92 | G4int cPDG=theQPDG.GetPDGCode(); |
---|
| 93 | if(cPDG>80000000) DefineQC(cPDG); |
---|
| 94 | else G4cerr<<"***G4QHadr:(QP) PDG="<<cPDG<<" use other constructor"<<G4endl; |
---|
| 95 | } |
---|
| 96 | } |
---|
| 97 | |
---|
| 98 | // Make sense Chipolino or Quasmon |
---|
| 99 | G4QHadron::G4QHadron(G4QContent QC, G4LorentzVector p): theQPDG(0),theMomentum(p),valQ(QC), |
---|
| 100 | nFragm(0),thePosition(0.,0.,0.),theCollisionCount(0),isSplit(false),Direction(true), |
---|
| 101 | Color(),AntiColor(),bindE(0.),formTime(0.) |
---|
| 102 | { |
---|
| 103 | G4int curPDG=valQ.GetSPDGCode(); |
---|
| 104 | if(curPDG==10&&valQ.GetBaryonNumber()>0) curPDG=valQ.GetZNSPDGCode(); |
---|
| 105 | if(curPDG&&curPDG!=10) theQPDG.SetPDGCode(curPDG); |
---|
| 106 | else theQPDG.InitByQCont(QC); |
---|
| 107 | } |
---|
| 108 | |
---|
| 109 | G4QHadron::G4QHadron(G4int PDGCode, G4double aMass, G4QContent QC) : |
---|
| 110 | theQPDG(PDGCode),theMomentum(0.,0.,0., aMass),valQ(QC),nFragm(0),thePosition(0.,0.,0.), |
---|
| 111 | theCollisionCount(0),isSplit(false),Direction(true),Color(),AntiColor(),bindE(0.), |
---|
| 112 | formTime(0.) |
---|
| 113 | {} |
---|
| 114 | |
---|
| 115 | G4QHadron::G4QHadron(G4QPDGCode QPDG, G4double aMass, G4QContent QC) : |
---|
| 116 | theQPDG(QPDG),theMomentum(0.,0.,0., aMass),valQ(QC),nFragm(0),thePosition(0.,0.,0.), |
---|
| 117 | theCollisionCount(0),isSplit(false),Direction(true),Color(),AntiColor(),bindE(0.), |
---|
| 118 | formTime(0.) |
---|
| 119 | {} |
---|
| 120 | |
---|
| 121 | G4QHadron::G4QHadron(G4int PDGCode, G4LorentzVector p, G4QContent QC) : |
---|
| 122 | theQPDG(PDGCode),theMomentum(p),valQ(QC),nFragm(0),thePosition(0.,0.,0.), |
---|
| 123 | theCollisionCount(0),isSplit(false),Direction(true),Color(),AntiColor(),bindE(0.), |
---|
| 124 | formTime(0.) |
---|
| 125 | {} |
---|
| 126 | |
---|
| 127 | G4QHadron::G4QHadron(G4QPDGCode QPDG, G4LorentzVector p, G4QContent QC) : |
---|
| 128 | theQPDG(QPDG),theMomentum(p),valQ(QC),nFragm(0),thePosition(0.,0.,0.), |
---|
| 129 | theCollisionCount(0),isSplit(false),Direction(true),Color(),AntiColor(),bindE(0.), |
---|
| 130 | formTime(0.) |
---|
| 131 | {} |
---|
| 132 | |
---|
| 133 | G4QHadron::G4QHadron(G4QParticle* pPart, G4double maxM) : |
---|
| 134 | theQPDG(pPart->GetQPDG()),theMomentum(0.,0.,0.,0.),nFragm(0),thePosition(0.,0.,0.), |
---|
| 135 | theCollisionCount(0),isSplit(false),Direction(true),Color(),AntiColor(),bindE(0.), |
---|
| 136 | formTime(0.) |
---|
| 137 | { |
---|
| 138 | #ifdef debug |
---|
| 139 | G4cout<<"G4QHadron is created & randomized with maxM="<<maxM<<G4endl; |
---|
| 140 | #endif |
---|
| 141 | G4int PDGCode = theQPDG.GetPDGCode(); |
---|
| 142 | if(PDGCode<2)G4cerr<<"***G4QHadron:(M) PDGC="<<PDGCode<<" use other constructor"<<G4endl; |
---|
| 143 | valQ=theQPDG.GetQuarkContent(); |
---|
| 144 | theMomentum.setE(RandomizeMass(pPart, maxM)); |
---|
| 145 | } |
---|
| 146 | |
---|
| 147 | G4QHadron::G4QHadron(const G4QHadron& right) |
---|
| 148 | { |
---|
| 149 | theMomentum = right.theMomentum; |
---|
| 150 | theQPDG = right.theQPDG; |
---|
| 151 | valQ = right.valQ; |
---|
| 152 | nFragm = right.nFragm; |
---|
| 153 | thePosition = right.thePosition; |
---|
| 154 | theCollisionCount = 0; |
---|
| 155 | isSplit = false; |
---|
| 156 | Direction = right.Direction; |
---|
| 157 | bindE = right.bindE; |
---|
| 158 | formTime = right.formTime; |
---|
| 159 | } |
---|
| 160 | |
---|
| 161 | G4QHadron::G4QHadron(const G4QHadron* right) |
---|
| 162 | { |
---|
| 163 | theMomentum = right->theMomentum; |
---|
| 164 | theQPDG = right->theQPDG; |
---|
| 165 | valQ = right->valQ; |
---|
| 166 | nFragm = right->nFragm; |
---|
| 167 | thePosition = right->thePosition; |
---|
| 168 | theCollisionCount = 0; |
---|
| 169 | isSplit = false; |
---|
| 170 | Direction = right->Direction; |
---|
| 171 | bindE = right->bindE; |
---|
| 172 | formTime = right->formTime; |
---|
| 173 | } |
---|
| 174 | |
---|
| 175 | G4QHadron::G4QHadron(const G4QHadron* right, G4int C, G4ThreeVector P, G4LorentzVector M) |
---|
| 176 | { |
---|
| 177 | theMomentum = M; |
---|
| 178 | theQPDG = right->theQPDG; |
---|
| 179 | valQ = right->valQ; |
---|
| 180 | nFragm = right->nFragm; |
---|
| 181 | thePosition = P; |
---|
| 182 | theCollisionCount = C; |
---|
| 183 | isSplit = false; |
---|
| 184 | Direction = right->Direction; |
---|
| 185 | bindE = right->bindE; |
---|
| 186 | formTime = right->formTime; |
---|
| 187 | } |
---|
| 188 | |
---|
| 189 | const G4QHadron& G4QHadron::operator=(const G4QHadron &right) |
---|
| 190 | { |
---|
| 191 | if(this != &right) // Beware of self assignment |
---|
| 192 | { |
---|
| 193 | theMomentum = right.theMomentum; |
---|
| 194 | theQPDG = right.theQPDG; |
---|
| 195 | valQ = right.valQ; |
---|
| 196 | nFragm = right.nFragm; |
---|
| 197 | thePosition = right.thePosition; |
---|
| 198 | theCollisionCount = 0; |
---|
| 199 | isSplit = false; |
---|
| 200 | Direction = right.Direction; |
---|
| 201 | bindE = right.bindE; |
---|
| 202 | } |
---|
| 203 | return *this; |
---|
| 204 | } |
---|
| 205 | |
---|
| 206 | G4QHadron::~G4QHadron() |
---|
| 207 | { |
---|
| 208 | std::list<G4QParton*>::iterator ipos = Color.begin(); |
---|
| 209 | std::list<G4QParton*>::iterator epos = Color.end(); |
---|
| 210 | for( ; ipos != epos; ipos++) {delete [] *ipos;} |
---|
| 211 | Color.clear(); |
---|
| 212 | |
---|
| 213 | ipos = AntiColor.begin(); |
---|
| 214 | epos = AntiColor.end(); |
---|
| 215 | for( ; ipos != epos; ipos++) {delete [] *ipos;} |
---|
| 216 | AntiColor.clear(); |
---|
| 217 | } |
---|
| 218 | |
---|
| 219 | // Define quark content of the particle with a particular PDG Code |
---|
| 220 | void G4QHadron::DefineQC(G4int PDGCode) |
---|
| 221 | { |
---|
| 222 | //G4cout<<"G4QHadron::DefineQC is called with PDGCode="<<PDGCode<<G4endl; |
---|
| 223 | G4int szn=PDGCode-90000000; |
---|
| 224 | G4int ds=0; |
---|
| 225 | G4int dz=0; |
---|
| 226 | G4int dn=0; |
---|
| 227 | if(szn<-100000) |
---|
| 228 | { |
---|
| 229 | G4int ns=(-szn)/1000000+1; |
---|
| 230 | szn+=ns*1000000; |
---|
| 231 | ds+=ns; |
---|
| 232 | } |
---|
| 233 | else if(szn<-100) |
---|
| 234 | { |
---|
| 235 | G4int nz=(-szn)/1000+1; |
---|
| 236 | szn+=nz*1000; |
---|
| 237 | dz+=nz; |
---|
| 238 | } |
---|
| 239 | else if(szn<0) |
---|
| 240 | { |
---|
| 241 | G4int nn=-szn; |
---|
| 242 | szn=0; |
---|
| 243 | dn+=nn; |
---|
| 244 | } |
---|
| 245 | G4int sz =szn/1000; |
---|
| 246 | G4int n =szn%1000; |
---|
| 247 | if(n>700) |
---|
| 248 | { |
---|
| 249 | n-=1000; |
---|
| 250 | dz--; |
---|
| 251 | } |
---|
| 252 | G4int z =sz%1000-dz; |
---|
| 253 | if(z>700) |
---|
| 254 | { |
---|
| 255 | z-=1000; |
---|
| 256 | ds--; |
---|
| 257 | } |
---|
| 258 | G4int Sq =sz/1000-ds; |
---|
| 259 | G4int zns=z+n+Sq; |
---|
| 260 | G4int Dq=n+zns; |
---|
| 261 | G4int Uq=z+zns; |
---|
| 262 | if (Dq<0&&Uq<0&&Sq<0)valQ=G4QContent(0 ,0 ,0 ,-Dq,-Uq,-Sq); |
---|
| 263 | else if (Uq<0&&Sq<0) valQ=G4QContent(Dq,0 ,0 ,0 ,-Uq,-Sq); |
---|
| 264 | else if (Dq<0&&Sq<0) valQ=G4QContent(0 ,Uq,0 ,-Dq,0 ,-Sq); |
---|
| 265 | else if (Dq<0&&Uq<0) valQ=G4QContent(0 ,0 ,Sq,-Dq,-Uq,0 ); |
---|
| 266 | else if (Uq<0) valQ=G4QContent(Dq,0 ,Sq,0 ,-Uq,0 ); |
---|
| 267 | else if (Sq<0) valQ=G4QContent(Dq,Uq,0 ,0 ,0 ,-Sq); |
---|
| 268 | else if (Dq<0) valQ=G4QContent(0 ,Uq,Sq,-Dq,0 ,0 ); |
---|
| 269 | else valQ=G4QContent(Dq,Uq,Sq,0 ,0 ,0 ); |
---|
| 270 | } |
---|
| 271 | |
---|
| 272 | // Redefine a Hadron with a new PDGCode |
---|
| 273 | void G4QHadron::SetQPDG(const G4QPDGCode& newQPDG) |
---|
| 274 | { |
---|
| 275 | theQPDG = newQPDG; |
---|
| 276 | G4int PDG= newQPDG.GetPDGCode(); |
---|
| 277 | G4int Q = newQPDG.GetQCode(); |
---|
| 278 | #ifdef debug |
---|
| 279 | G4cout<<"G4QHadron::SetQPDG is called with PDGCode="<<PDG<<", QCode="<<Q<<G4endl; |
---|
| 280 | #endif |
---|
| 281 | if (Q>-1) valQ=theQPDG.GetQuarkContent(); |
---|
| 282 | else if(PDG>80000000) DefineQC(PDG); |
---|
| 283 | else |
---|
| 284 | { |
---|
| 285 | G4cerr<<"***G4QHadron::SetQPDG: QPDG="<<newQPDG<<G4endl; |
---|
| 286 | throw G4QException("***G4QHadron::SetQPDG: Impossible QPDG Probably a Chipolino"); |
---|
| 287 | } |
---|
| 288 | } |
---|
| 289 | |
---|
| 290 | // Decay of Hadron In2Particles f&s, f is in respect to the direction of HadronMomentumDir |
---|
| 291 | G4bool G4QHadron::RelDecayIn2(G4LorentzVector& f4Mom, G4LorentzVector& s4Mom, |
---|
| 292 | G4LorentzVector& dir, G4double maxCost, G4double minCost) |
---|
| 293 | {// =================================================================== |
---|
| 294 | G4double fM2 = f4Mom.m2(); |
---|
| 295 | G4double fM = sqrt(fM2); // Mass of the 1st Hadron |
---|
| 296 | G4double sM2 = s4Mom.m2(); |
---|
| 297 | G4double sM = sqrt(sM2); // Mass of the 2nd Hadron |
---|
| 298 | G4double iM2 = theMomentum.m2(); |
---|
| 299 | G4double iM = sqrt(iM2); // Mass of the decaying hadron |
---|
| 300 | G4double vP = theMomentum.rho(); // Momentum of the decaying hadron |
---|
| 301 | G4double dE = theMomentum.e(); // Energy of the decaying hadron |
---|
| 302 | if(dE<vP) |
---|
| 303 | { |
---|
| 304 | G4cerr<<"***G4QHad::RelDecIn2: Tachionic 4-mom="<<theMomentum<<", E-p="<<dE-vP<<G4endl; |
---|
| 305 | G4double accuracy=.000001*vP; |
---|
| 306 | G4double emodif=std::fabs(dE-vP); |
---|
| 307 | //if(emodif<accuracy) |
---|
| 308 | //{ |
---|
| 309 | G4cerr<<"G4QHadron::RelDecIn2: *Boost* E-p shift is corrected to "<<emodif<<G4endl; |
---|
| 310 | theMomentum.setE(vP+emodif+.01*accuracy); |
---|
| 311 | //} |
---|
| 312 | } |
---|
| 313 | G4ThreeVector ltb = theMomentum.boostVector();// Boost vector for backward Lorentz Trans. |
---|
| 314 | G4ThreeVector ltf = -ltb; // Boost vector for forward Lorentz Trans. |
---|
| 315 | G4LorentzVector cdir = dir; // A copy to make a transformation to CMS |
---|
| 316 | #ifdef ppdebug |
---|
| 317 | if(cdir.e()+.001<cdir.rho()) G4cerr<<"*G4QH::RDIn2:*Boost* cd4M="<<cdir<<",e-p=" |
---|
| 318 | <<cdir.e()-cdir.rho()<<G4endl; |
---|
| 319 | #endif |
---|
| 320 | cdir.boost(ltf); // Direction transpormed to CMS of the Momentum |
---|
| 321 | G4ThreeVector vdir = cdir.vect(); // 3-Vector of the direction-particle |
---|
| 322 | #ifdef ppdebug |
---|
| 323 | G4cout<<"G4QHad::RelDI2:dir="<<dir<<",ltf="<<ltf<<",cdir="<<cdir<<",vdir="<<vdir<<G4endl; |
---|
| 324 | #endif |
---|
| 325 | G4ThreeVector vx(0.,0.,1.); // Ort in the direction of the reference particle |
---|
| 326 | G4ThreeVector vy(0.,1.,0.); // First ort orthogonal to the direction |
---|
| 327 | G4ThreeVector vz(1.,0.,0.); // Second ort orthoganal to the direction |
---|
| 328 | if(vdir.mag2() > 0.) // the refference particle isn't at rest in CMS |
---|
| 329 | { |
---|
| 330 | vx = vdir.unit(); // Ort in the direction of the reference particle |
---|
| 331 | #ifdef ppdebug |
---|
| 332 | G4cout<<"G4QH::RelDecIn2:Vx="<<vx<<",M="<<theMomentum<<",d="<<dir<<",c="<<cdir<<G4endl; |
---|
| 333 | #endif |
---|
| 334 | G4ThreeVector vv= vx.orthogonal(); // Not normed orthogonal vector (!) |
---|
| 335 | vy = vv.unit(); // First ort orthogonal to the direction |
---|
| 336 | vz = vx.cross(vy); // Second ort orthoganal to the direction |
---|
| 337 | } |
---|
| 338 | #ifdef ppdebug |
---|
| 339 | G4cout<<"G4QHad::RelDecIn2:iM="<<iM<<"=>fM="<<fM<<"+sM="<<sM<<",ob="<<vx<<vy<<vz<<G4endl; |
---|
| 340 | #endif |
---|
| 341 | if(maxCost> 1.) maxCost= 1.; |
---|
| 342 | if(minCost<-1.) minCost=-1.; |
---|
| 343 | if(maxCost<-1.) maxCost=-1.; |
---|
| 344 | if(minCost> 1.) minCost= 1.; |
---|
| 345 | if(minCost> maxCost) minCost=maxCost; |
---|
| 346 | if(fabs(iM-fM-sM)<.00000001) |
---|
| 347 | { |
---|
| 348 | G4double fR=fM/iM; |
---|
| 349 | G4double sR=sM/iM; |
---|
| 350 | f4Mom=fR*theMomentum; |
---|
| 351 | s4Mom=sR*theMomentum; |
---|
| 352 | return true; |
---|
| 353 | } |
---|
| 354 | else if (iM+.001<fM+sM || iM==0.) |
---|
| 355 | {//@@ Later on make a quark content check for the decay |
---|
| 356 | G4cerr<<"***G4QH::RelDecIn2: fM="<<fM<<"+sM="<<sM<<">iM="<<iM<<",d="<<iM-fM-sM<<G4endl; |
---|
| 357 | return false; |
---|
| 358 | } |
---|
| 359 | G4double d2 = iM2-fM2-sM2; |
---|
| 360 | G4double p2 = (d2*d2/4.-fM2*sM2)/iM2; // Decay momentum(^2) in CMS of Quasmon |
---|
| 361 | if(p2<0.) |
---|
| 362 | { |
---|
| 363 | #ifdef ppdebug |
---|
| 364 | G4cout<<"**G4QH:RDIn2:p2="<<p2<<"<0,d2^2="<<d2*d2/4.<<"<4*fM2*sM2="<<4*fM2*sM2<<G4endl; |
---|
| 365 | #endif |
---|
| 366 | p2=0.; |
---|
| 367 | } |
---|
| 368 | G4double p = sqrt(p2); |
---|
| 369 | G4double ct = maxCost; |
---|
| 370 | if(maxCost>minCost) |
---|
| 371 | { |
---|
| 372 | G4double dcost=maxCost-minCost; |
---|
| 373 | ct = minCost+dcost*G4UniformRand(); |
---|
| 374 | } |
---|
| 375 | G4double phi= twopi*G4UniformRand(); // @@ Change 360.*deg to M_TWOPI (?) |
---|
| 376 | G4double ps=0.; |
---|
| 377 | if(fabs(ct)<1.) ps = p * sqrt(1.-ct*ct); |
---|
| 378 | else |
---|
| 379 | { |
---|
| 380 | #ifdef ppdebug |
---|
| 381 | G4cout<<"**G4QH::RDIn2:ct="<<ct<<",mac="<<maxCost<<",mic="<<minCost<<G4endl; |
---|
| 382 | //throw G4QException("***G4QHadron::RDIn2: bad cos(theta)"); |
---|
| 383 | #endif |
---|
| 384 | if(ct>1.) ct=1.; |
---|
| 385 | if(ct<-1.) ct=-1.; |
---|
| 386 | } |
---|
| 387 | G4ThreeVector pVect=(ps*sin(phi))*vz+(ps*cos(phi))*vy+p*ct*vx; |
---|
| 388 | #ifdef ppdebug |
---|
| 389 | G4cout<<"G4QH::RelDIn2:ct="<<ct<<",p="<<p<<",ps="<<ps<<",ph="<<phi<<",v="<<pVect<<G4endl; |
---|
| 390 | #endif |
---|
| 391 | |
---|
| 392 | f4Mom.setVect(pVect); |
---|
| 393 | f4Mom.setE(sqrt(fM2+p2)); |
---|
| 394 | s4Mom.setVect((-1)*pVect); |
---|
| 395 | s4Mom.setE(sqrt(sM2+p2)); |
---|
| 396 | |
---|
| 397 | #ifdef ppdebug |
---|
| 398 | G4cout<<"G4QHadr::RelDecIn2:p2="<<p2<<",v="<<ltb<<",f4M="<<f4Mom<<" + s4M="<<s4Mom<<" = " |
---|
| 399 | <<f4Mom+s4Mom<<", M="<<iM<<G4endl; |
---|
| 400 | #endif |
---|
| 401 | if(f4Mom.e()+.001<f4Mom.rho())G4cerr<<"*G4QH::RDIn2:*Boost* f4M="<<f4Mom<<",e-p=" |
---|
| 402 | <<f4Mom.e()-f4Mom.rho()<<G4endl; |
---|
| 403 | f4Mom.boost(ltb); // Lor.Trans. of 1st hadron back to LS |
---|
| 404 | if(s4Mom.e()+.001<s4Mom.rho())G4cerr<<"*G4QH::RDIn2:*Boost* s4M="<<s4Mom<<",e-p=" |
---|
| 405 | <<s4Mom.e()-s4Mom.rho()<<G4endl; |
---|
| 406 | s4Mom.boost(ltb); // Lor.Trans. of 2nd hadron back to LS |
---|
| 407 | #ifdef ppdebug |
---|
| 408 | G4cout<<"G4QHadron::RelDecayIn2:Output, f4Mom="<<f4Mom<<" + s4Mom="<<s4Mom<<" = " |
---|
| 409 | <<f4Mom+s4Mom<<", d4M="<<theMomentum-f4Mom-s4Mom<<G4endl; |
---|
| 410 | #endif |
---|
| 411 | return true; |
---|
| 412 | } // End of "RelDecayIn2" |
---|
| 413 | |
---|
| 414 | // Decay of Hadron In2Particles f&s, f w/r/to dN/dO [cp>0: ~cost^cp, cp<0: ~(1-cost)^(-cp)] |
---|
| 415 | G4bool G4QHadron::CopDecayIn2(G4LorentzVector& f4Mom, G4LorentzVector& s4Mom, |
---|
| 416 | G4LorentzVector& dir, G4double cosp) |
---|
| 417 | {// =================================================================== |
---|
| 418 | G4double fM2 = f4Mom.m2(); |
---|
| 419 | G4double fM = sqrt(fM2); // Mass of the 1st Hadron |
---|
| 420 | G4double sM2 = s4Mom.m2(); |
---|
| 421 | G4double sM = sqrt(sM2); // Mass of the 2nd Hadron |
---|
| 422 | G4double iM2 = theMomentum.m2(); |
---|
| 423 | G4double iM = sqrt(iM2); // Mass of the decaying hadron |
---|
| 424 | G4double vP = theMomentum.rho(); // Momentum of the decaying hadron |
---|
| 425 | G4double dE = theMomentum.e(); // Energy of the decaying hadron |
---|
| 426 | G4bool neg=false; // Negative (backward) distribution of t |
---|
| 427 | if(cosp<0) |
---|
| 428 | { |
---|
| 429 | cosp=-cosp; |
---|
| 430 | neg=true; |
---|
| 431 | } |
---|
| 432 | if(dE<vP) |
---|
| 433 | { |
---|
| 434 | G4cerr<<"***G4QHad::CopDecIn2: Tachionic 4-mom="<<theMomentum<<", E-p="<<dE-vP<<G4endl; |
---|
| 435 | G4double accuracy=.000001*vP; |
---|
| 436 | G4double emodif=std::fabs(dE-vP); |
---|
| 437 | //if(emodif<accuracy) |
---|
| 438 | //{ |
---|
| 439 | G4cerr<<"G4QHadron::CopDecIn2: *Boost* E-p shift is corrected to "<<emodif<<G4endl; |
---|
| 440 | theMomentum.setE(vP+emodif+.01*accuracy); |
---|
| 441 | //} |
---|
| 442 | } |
---|
| 443 | G4ThreeVector ltb = theMomentum.boostVector();// Boost vector for backward Lorentz Trans. |
---|
| 444 | G4ThreeVector ltf = -ltb; // Boost vector for forward Lorentz Trans. |
---|
| 445 | G4LorentzVector cdir = dir; // A copy to make a transformation to CMS |
---|
| 446 | #ifdef ppdebug |
---|
| 447 | if(cdir.e()+.001<cdir.rho()) G4cerr<<"*G4QH::RDIn2:*Boost* cd4M="<<cdir<<",e-p=" |
---|
| 448 | <<cdir.e()-cdir.rho()<<G4endl; |
---|
| 449 | #endif |
---|
| 450 | cdir.boost(ltf); // Direction transpormed to CMS of the Momentum |
---|
| 451 | G4ThreeVector vdir = cdir.vect(); // 3-Vector of the direction-particle |
---|
| 452 | #ifdef ppdebug |
---|
| 453 | G4cout<<"G4QHad::CopDI2:dir="<<dir<<",ltf="<<ltf<<",cdir="<<cdir<<",vdir="<<vdir<<G4endl; |
---|
| 454 | #endif |
---|
| 455 | G4ThreeVector vx(0.,0.,1.); // Ort in the direction of the reference particle |
---|
| 456 | G4ThreeVector vy(0.,1.,0.); // First ort orthogonal to the direction |
---|
| 457 | G4ThreeVector vz(1.,0.,0.); // Second ort orthoganal to the direction |
---|
| 458 | if(vdir.mag2() > 0.) // the refference particle isn't at rest in CMS |
---|
| 459 | { |
---|
| 460 | vx = vdir.unit(); // Ort in the direction of the reference particle |
---|
| 461 | #ifdef ppdebug |
---|
| 462 | G4cout<<"G4QH::CopDecIn2:Vx="<<vx<<",M="<<theMomentum<<",d="<<dir<<",c="<<cdir<<G4endl; |
---|
| 463 | #endif |
---|
| 464 | G4ThreeVector vv= vx.orthogonal(); // Not normed orthogonal vector (!) |
---|
| 465 | vy = vv.unit(); // First ort orthogonal to the direction |
---|
| 466 | vz = vx.cross(vy); // Second ort orthoganal to the direction |
---|
| 467 | } |
---|
| 468 | #ifdef ppdebug |
---|
| 469 | G4cout<<"G4QHad::CopDecIn2:iM="<<iM<<"=>fM="<<fM<<"+sM="<<sM<<",ob="<<vx<<vy<<vz<<G4endl; |
---|
| 470 | #endif |
---|
| 471 | if(fabs(iM-fM-sM)<.00000001) |
---|
| 472 | { |
---|
| 473 | G4double fR=fM/iM; |
---|
| 474 | G4double sR=sM/iM; |
---|
| 475 | f4Mom=fR*theMomentum; |
---|
| 476 | s4Mom=sR*theMomentum; |
---|
| 477 | return true; |
---|
| 478 | } |
---|
| 479 | else if (iM+.001<fM+sM || iM==0.) |
---|
| 480 | {//@@ Later on make a quark content check for the decay |
---|
| 481 | G4cerr<<"***G4QH::CopDecIn2: fM="<<fM<<"+sM="<<sM<<">iM="<<iM<<",d="<<iM-fM-sM<<G4endl; |
---|
| 482 | return false; |
---|
| 483 | } |
---|
| 484 | G4double d2 = iM2-fM2-sM2; |
---|
| 485 | G4double p2 = (d2*d2/4.-fM2*sM2)/iM2; // Decay momentum(^2) in CMS of Quasmon |
---|
| 486 | if(p2<0.) |
---|
| 487 | { |
---|
| 488 | #ifdef ppdebug |
---|
| 489 | G4cout<<"*G4QH:CopDI2:p2="<<p2<<"<0,d4/4="<<d2*d2/4.<<"<4*fM2*sM2="<<4*fM2*sM2<<G4endl; |
---|
| 490 | #endif |
---|
| 491 | p2=0.; |
---|
| 492 | } |
---|
| 493 | G4double p = sqrt(p2); |
---|
| 494 | G4double ct = 0; |
---|
| 495 | G4double rn = pow(G4UniformRand(),cosp+1.); |
---|
| 496 | if(neg) ct = rn+rn-1.; // More backward than forward |
---|
| 497 | else ct = 1.-rn-rn; // More forward than backward |
---|
| 498 | // |
---|
| 499 | G4double phi= twopi*G4UniformRand(); // @@ Change 360.*deg to M_TWOPI (?) |
---|
| 500 | G4double ps=0.; |
---|
| 501 | if(fabs(ct)<1.) ps = p * sqrt(1.-ct*ct); |
---|
| 502 | else |
---|
| 503 | { |
---|
| 504 | #ifdef ppdebug |
---|
| 505 | G4cout<<"**G4QH::CopDecayIn2:ct="<<ct<<",mac="<<maxCost<<",mic="<<minCost<<G4endl; |
---|
| 506 | //throw G4QException("***G4QHadron::RDIn2: bad cos(theta)"); |
---|
| 507 | #endif |
---|
| 508 | if(ct>1.) ct=1.; |
---|
| 509 | if(ct<-1.) ct=-1.; |
---|
| 510 | } |
---|
| 511 | G4ThreeVector pVect=(ps*sin(phi))*vz+(ps*cos(phi))*vy+p*ct*vx; |
---|
| 512 | #ifdef ppdebug |
---|
| 513 | G4cout<<"G4QH::CopDIn2:ct="<<ct<<",p="<<p<<",ps="<<ps<<",ph="<<phi<<",v="<<pVect<<G4endl; |
---|
| 514 | #endif |
---|
| 515 | |
---|
| 516 | f4Mom.setVect(pVect); |
---|
| 517 | f4Mom.setE(sqrt(fM2+p2)); |
---|
| 518 | s4Mom.setVect((-1)*pVect); |
---|
| 519 | s4Mom.setE(sqrt(sM2+p2)); |
---|
| 520 | |
---|
| 521 | #ifdef ppdebug |
---|
| 522 | G4cout<<"G4QHadr::CopDecIn2:p2="<<p2<<",v="<<ltb<<",f4M="<<f4Mom<<" + s4M="<<s4Mom<<" = " |
---|
| 523 | <<f4Mom+s4Mom<<", M="<<iM<<G4endl; |
---|
| 524 | #endif |
---|
| 525 | if(f4Mom.e()+.001<f4Mom.rho())G4cerr<<"*G4QH::RDIn2:*Boost* f4M="<<f4Mom<<",e-p=" |
---|
| 526 | <<f4Mom.e()-f4Mom.rho()<<G4endl; |
---|
| 527 | f4Mom.boost(ltb); // Lor.Trans. of 1st hadron back to LS |
---|
| 528 | if(s4Mom.e()+.001<s4Mom.rho())G4cerr<<"*G4QH::RDIn2:*Boost* s4M="<<s4Mom<<",e-p=" |
---|
| 529 | <<s4Mom.e()-s4Mom.rho()<<G4endl; |
---|
| 530 | s4Mom.boost(ltb); // Lor.Trans. of 2nd hadron back to LS |
---|
| 531 | #ifdef ppdebug |
---|
| 532 | G4cout<<"G4QHadron::CopDecayIn2:Output, f4Mom="<<f4Mom<<" + s4Mom="<<s4Mom<<" = " |
---|
| 533 | <<f4Mom+s4Mom<<", d4M="<<theMomentum-f4Mom-s4Mom<<G4endl; |
---|
| 534 | #endif |
---|
| 535 | return true; |
---|
| 536 | } // End of "CopDecayIn2" |
---|
| 537 | |
---|
| 538 | // Decay of the Hadron in 2 particles (f + s) |
---|
| 539 | G4bool G4QHadron::DecayIn2(G4LorentzVector& f4Mom, G4LorentzVector& s4Mom) |
---|
| 540 | {// =================================================================== |
---|
| 541 | G4double fM2 = f4Mom.m2(); |
---|
| 542 | if(fM2<0.) fM2=0.; |
---|
| 543 | G4double fM = sqrt(fM2); // Mass of the 1st Hadron |
---|
| 544 | G4double sM2 = s4Mom.m2(); |
---|
| 545 | if(sM2<0.) sM2=0.; |
---|
| 546 | G4double sM = sqrt(sM2); // Mass of the 2nd Hadron |
---|
| 547 | G4double iM2 = theMomentum.m2(); |
---|
| 548 | if(iM2<0.) iM2=0.; |
---|
| 549 | G4double iM = sqrt(iM2); // Mass of the decaying hadron |
---|
| 550 | #ifdef debug |
---|
| 551 | G4cout<<"G4QHadron::DecIn2: iM="<<iM<<" => fM="<<fM<<" + sM="<<sM<<" = "<<fM+sM<<G4endl; |
---|
| 552 | #endif |
---|
| 553 | //@@ Later on make a quark content check for the decay |
---|
| 554 | if (fabs(iM-fM-sM)<.0000001) |
---|
| 555 | { |
---|
| 556 | G4double fR=fM/iM; |
---|
| 557 | G4double sR=sM/iM; |
---|
| 558 | f4Mom=fR*theMomentum; |
---|
| 559 | s4Mom=sR*theMomentum; |
---|
| 560 | return true; |
---|
| 561 | } |
---|
| 562 | else if (iM+.001<fM+sM || iM==0.) |
---|
| 563 | { |
---|
| 564 | #ifdef pdebug |
---|
| 565 | G4cerr<<"***G4QHadron::DecayIn2*** fM="<<fM<<" + sM="<<sM<<"="<<fM+sM<<" > iM="<<iM |
---|
| 566 | <<", d="<<iM-fM-sM<<G4endl; |
---|
| 567 | #endif |
---|
| 568 | return false; |
---|
| 569 | } |
---|
| 570 | |
---|
| 571 | G4double d2 = iM2-fM2-sM2; |
---|
| 572 | G4double p2 = (d2*d2/4.-fM2*sM2)/iM2; // Decay momentum(^2) in CMS of Quasmon |
---|
| 573 | if (p2<0.) |
---|
| 574 | { |
---|
| 575 | #ifdef debug |
---|
| 576 | G4cerr<<"***G4QH::DI2:p2="<<p2<<"<0,d2^2="<<d2*d2/4.<<"<4*fM2*sM2="<<4*fM2*sM2<<G4endl; |
---|
| 577 | #endif |
---|
| 578 | p2=0.; |
---|
| 579 | } |
---|
| 580 | G4double p = sqrt(p2); |
---|
| 581 | G4double ct = 1.-2*G4UniformRand(); |
---|
| 582 | #ifdef debug |
---|
| 583 | G4cout<<"G4QHadron::DecayIn2: ct="<<ct<<", p="<<p<<G4endl; |
---|
| 584 | #endif |
---|
| 585 | G4double phi= twopi*G4UniformRand(); // @@ Change 360.*deg to M_TWOPI (?) |
---|
| 586 | G4double ps = p * sqrt(1.-ct*ct); |
---|
| 587 | G4ThreeVector pVect(ps*sin(phi),ps*cos(phi),p*ct); |
---|
| 588 | |
---|
| 589 | f4Mom.setVect(pVect); |
---|
| 590 | f4Mom.setE(sqrt(fM2+p2)); |
---|
| 591 | s4Mom.setVect((-1)*pVect); |
---|
| 592 | s4Mom.setE(sqrt(sM2+p2)); |
---|
| 593 | |
---|
| 594 | if(theMomentum.e()<theMomentum.rho()) |
---|
| 595 | { |
---|
| 596 | G4cerr<<"*G4QH::DecIn2:*Boost* 4M="<<theMomentum<<",e-p=" |
---|
| 597 | <<theMomentum.e()-theMomentum.rho()<<G4endl; |
---|
| 598 | //throw G4QException("G4QHadron::DecayIn2: Decay of particle with zero mass") |
---|
| 599 | theMomentum.setE(1.0000001*theMomentum.rho()); |
---|
| 600 | } |
---|
| 601 | G4double vP = theMomentum.rho(); // Momentum of the decaying hadron |
---|
| 602 | G4double dE = theMomentum.e(); // Energy of the decaying hadron |
---|
| 603 | if(dE<vP) |
---|
| 604 | { |
---|
| 605 | G4cerr<<"***G4QHad::RelDecIn2: Tachionic 4-mom="<<theMomentum<<", E-p="<<dE-vP<<G4endl; |
---|
| 606 | G4double accuracy=.000001*vP; |
---|
| 607 | G4double emodif=std::fabs(dE-vP); |
---|
| 608 | if(emodif<accuracy) |
---|
| 609 | { |
---|
| 610 | G4cerr<<"G4QHadron::DecayIn2: *Boost* E-p shift is corrected to "<<emodif<<G4endl; |
---|
| 611 | theMomentum.setE(vP+emodif+.01*accuracy); |
---|
| 612 | } |
---|
| 613 | } |
---|
| 614 | G4ThreeVector ltb = theMomentum.boostVector(); // Boost vector for backward Lor.Trans. |
---|
| 615 | #ifdef pdebug |
---|
| 616 | G4cout<<"G4QHadron::DecIn2:LorTrans v="<<ltb<<",f4Mom="<<f4Mom<<",s4Mom="<<s4Mom<<G4endl; |
---|
| 617 | #endif |
---|
| 618 | if(f4Mom.e()+.001<f4Mom.rho())G4cerr<<"*G4QH::DecIn2:*Boost* f4M="<<f4Mom<<G4endl; |
---|
| 619 | f4Mom.boost(ltb); // Lor.Trans. of 1st hadron back to LS |
---|
| 620 | if(s4Mom.e()+.001<s4Mom.rho())G4cerr<<"*G4QH::DecIn2:*Boost* s4M="<<s4Mom<<G4endl; |
---|
| 621 | s4Mom.boost(ltb); // Lor.Trans. of 2nd hadron back to LS |
---|
| 622 | #ifdef pdebug |
---|
| 623 | G4cout<<"G4QHadron::DecayIn2: ROOT OUTPUT f4Mom="<<f4Mom<<", s4Mom="<<s4Mom<<G4endl; |
---|
| 624 | #endif |
---|
| 625 | return true; |
---|
| 626 | } // End of "DecayIn2" |
---|
| 627 | |
---|
| 628 | // Correction for the Hadron + fr decay in case of the new corM mass of the Hadron |
---|
| 629 | G4bool G4QHadron::CorMDecayIn2(G4double corM, G4LorentzVector& fr4Mom) |
---|
| 630 | {// =============================================================== |
---|
| 631 | G4double fM = fr4Mom.m(); // Mass of the Fragment |
---|
| 632 | G4LorentzVector comp=theMomentum+fr4Mom; // 4Mom of the decaying compound system |
---|
| 633 | G4double iM = comp.m(); // mass of the decaying compound system |
---|
| 634 | #ifdef pdebug |
---|
| 635 | G4cout<<"G4QH::CMDIn2: iM="<<iM<<comp<<"=>fM="<<fM<<"+corM="<<corM<<"="<<fM+corM<<G4endl; |
---|
| 636 | #endif |
---|
| 637 | G4double dE=iM-fM-corM; |
---|
| 638 | //@@ Later on make a quark content check for the decay |
---|
| 639 | if (fabs(dE)<.001) |
---|
| 640 | { |
---|
| 641 | G4double fR=fM/iM; |
---|
| 642 | G4double cR=corM/iM; |
---|
| 643 | fr4Mom=fR*comp; |
---|
| 644 | theMomentum=cR*comp; |
---|
| 645 | return true; |
---|
| 646 | } |
---|
| 647 | else if (dE<-.001 || iM==0.) |
---|
| 648 | { |
---|
| 649 | G4cerr<<"***G4QH::CorMDIn2***fM="<<fM<<" + cM="<<corM<<" > iM="<<iM<<",d="<<dE<<G4endl; |
---|
| 650 | return false; |
---|
| 651 | } |
---|
| 652 | G4double corM2= corM*corM; |
---|
| 653 | G4double fM2 = fM*fM; |
---|
| 654 | G4double iM2 = iM*iM; |
---|
| 655 | G4double d2 = iM2-fM2-corM2; |
---|
| 656 | G4double p2 = (d2*d2/4.-fM2*corM2)/iM2; // Decay momentum(^2) in CMS of Quasmon |
---|
| 657 | if (p2<0.) |
---|
| 658 | { |
---|
| 659 | #ifdef pdebug |
---|
| 660 | G4cerr<<"**G4QH::CMDI2:p2="<<p2<<"<0,d="<<d2*d2/4.<<"<4*fM2*hM2="<<4*fM2*corM2<<G4endl; |
---|
| 661 | #endif |
---|
| 662 | p2=0.; |
---|
| 663 | } |
---|
| 664 | G4double p = sqrt(p2); |
---|
| 665 | if(comp.e()<comp.rho())G4cerr<<"*G4QH::CorMDecayIn2:*Boost* comp4M="<<comp<<",e-p=" |
---|
| 666 | <<comp.e()-comp.rho()<<G4endl; |
---|
| 667 | G4ThreeVector ltb = comp.boostVector(); // Boost vector for backward Lor.Trans. |
---|
| 668 | G4ThreeVector ltf = -ltb; // Boost vector for forward Lorentz Trans. |
---|
| 669 | G4LorentzVector cm4Mom=fr4Mom; // Copy of fragment 4Mom to transform to CMS |
---|
| 670 | if(cm4Mom.e()<cm4Mom.rho()) |
---|
| 671 | { |
---|
| 672 | G4cerr<<"*G4QH::CorMDecIn2:*Boost* c4M="<<cm4Mom<<G4endl; |
---|
| 673 | cm4Mom.setE(1.0000001*cm4Mom.rho()); |
---|
| 674 | } |
---|
| 675 | cm4Mom.boost(ltf); // Now it is in CMS (Forward Lor.Trans.) |
---|
| 676 | G4double pfx= cm4Mom.px(); |
---|
| 677 | G4double pfy= cm4Mom.py(); |
---|
| 678 | G4double pfz= cm4Mom.pz(); |
---|
| 679 | G4double pt2= pfx*pfx+pfy*pfy; |
---|
| 680 | G4double tx=0.; |
---|
| 681 | G4double ty=0.; |
---|
| 682 | if(pt2<=0.) |
---|
| 683 | { |
---|
| 684 | G4double phi= 360.*deg*G4UniformRand(); // @@ Change 360.*deg to M_TWOPI (?) |
---|
| 685 | tx=sin(phi); |
---|
| 686 | ty=cos(phi); |
---|
| 687 | } |
---|
| 688 | else |
---|
| 689 | { |
---|
| 690 | G4double pt=sqrt(pt2); |
---|
| 691 | tx=pfx/pt; |
---|
| 692 | ty=pfy/pt; |
---|
| 693 | } |
---|
| 694 | G4double pc2=pt2+pfz*pfz; |
---|
| 695 | G4double ct=0.; |
---|
| 696 | if(pc2<=0.) |
---|
| 697 | { |
---|
| 698 | G4double rnd= G4UniformRand(); |
---|
| 699 | ct=1.-rnd-rnd; |
---|
| 700 | } |
---|
| 701 | else |
---|
| 702 | { |
---|
| 703 | G4double pc=sqrt(pc2); |
---|
| 704 | ct=pfz/pc; |
---|
| 705 | } |
---|
| 706 | #ifdef debug |
---|
| 707 | G4cout<<"G4QHadron::CorMDecayIn2: ct="<<ct<<", p="<<p<<G4endl; |
---|
| 708 | #endif |
---|
| 709 | G4double ps = p * sqrt(1.-ct*ct); |
---|
| 710 | G4ThreeVector pVect(ps*tx,ps*ty,p*ct); |
---|
| 711 | fr4Mom.setVect(pVect); |
---|
| 712 | fr4Mom.setE(sqrt(fM2+p2)); |
---|
| 713 | theMomentum.setVect((-1)*pVect); |
---|
| 714 | theMomentum.setE(sqrt(corM2+p2)); |
---|
| 715 | #ifdef pdebug |
---|
| 716 | G4LorentzVector dif2=comp-fr4Mom-theMomentum; |
---|
| 717 | G4cout<<"G4QH::CorMDIn2:c="<<comp<<"-f="<<fr4Mom<<"-4M="<<theMomentum<<"="<<dif2<<G4endl; |
---|
| 718 | #endif |
---|
| 719 | if(fr4Mom.e()+.001<fr4Mom.rho())G4cerr<<"*G4QH::CorMDecIn2:*Boost*fr4M="<<fr4Mom<<G4endl; |
---|
| 720 | fr4Mom.boost(ltb); // Lor.Trans. of the Fragment back to LS |
---|
| 721 | if(theMomentum.e()<theMomentum.rho()) |
---|
| 722 | { |
---|
| 723 | G4cerr<<"*G4QH::CMDI2:4="<<theMomentum<<G4endl; |
---|
| 724 | theMomentum.setE(1.0000001*theMomentum.rho()); |
---|
| 725 | } |
---|
| 726 | theMomentum.boost(ltb); // Lor.Trans. of the Hadron back to LS |
---|
| 727 | #ifdef pdebug |
---|
| 728 | G4LorentzVector dif3=comp-fr4Mom-theMomentum; |
---|
| 729 | G4cout<<"G4QH::CorMDecIn2:OUTPUT:f4M="<<fr4Mom<<",h4M="<<theMomentum<<"d="<<dif3<<G4endl; |
---|
| 730 | #endif |
---|
| 731 | return true; |
---|
| 732 | } // End of "CorMDecayIn2" |
---|
| 733 | |
---|
| 734 | |
---|
| 735 | // Fragment fr4Mom louse energy corE and transfer it to This Hadron |
---|
| 736 | G4bool G4QHadron::CorEDecayIn2(G4double corE, G4LorentzVector& fr4Mom) |
---|
| 737 | {// =============================================================== |
---|
| 738 | G4double fE = fr4Mom.m(); // Energy of the Fragment |
---|
| 739 | #ifdef pdebug |
---|
| 740 | G4cout<<"G4QH::CorEDecIn2:fE="<<fE<<fr4Mom<<">corE="<<corE<<",h4M="<<theMomentum<<G4endl; |
---|
| 741 | #endif |
---|
| 742 | if (fE+.001<=corE) |
---|
| 743 | { |
---|
| 744 | #ifdef pdebug |
---|
| 745 | G4cerr<<"***G4QHadron::CorEDecIn2*** fE="<<fE<<"<corE="<<corE<<", d="<<corE-fE<<G4endl; |
---|
| 746 | #endif |
---|
| 747 | return false; |
---|
| 748 | } |
---|
| 749 | G4double fM2=fr4Mom.m2(); // Squared Mass of the Fragment |
---|
| 750 | if(fM2<0.) fM2=0.; |
---|
| 751 | G4double iPx=fr4Mom.px(); // Initial Px of the Fragment |
---|
| 752 | G4double iPy=fr4Mom.py(); // Initial Py of the Fragment |
---|
| 753 | G4double iPz=fr4Mom.pz(); // Initial Pz of the Fragment |
---|
| 754 | G4double fP2=iPx*iPx+iPy*iPy+iPz*iPz; // Initial Squared 3-momentum of the Fragment |
---|
| 755 | G4double finE = fE - corE; // Final energy of the fragment |
---|
| 756 | G4double rP = sqrt((finE*finE-fM2)/fP2); // Reduction factor for the momentum |
---|
| 757 | G4double fPx=iPx*rP; |
---|
| 758 | G4double fPy=iPy*rP; |
---|
| 759 | G4double fPz=iPz*rP; |
---|
| 760 | fr4Mom= G4LorentzVector(fPx,fPy,fPz,finE); |
---|
| 761 | G4double Px=theMomentum.px()+iPx-fPx; |
---|
| 762 | G4double Py=theMomentum.py()+iPy-fPy; |
---|
| 763 | G4double Pz=theMomentum.pz()+iPz-fPz; |
---|
| 764 | G4double mE=theMomentum.e(); |
---|
| 765 | ///////////G4double mM2=theMomentum.m2(); |
---|
| 766 | theMomentum= G4LorentzVector(Px,Py,Pz,mE+corE); |
---|
| 767 | #ifdef pdebug |
---|
| 768 | G4double difF=fr4Mom.m2()-fM2; |
---|
| 769 | G4cout<<"G4QH::CorEDecIn2: dF="<<difF<<",out:"<<theMomentum<<fr4Mom<<G4endl; |
---|
| 770 | #endif |
---|
| 771 | return true; |
---|
| 772 | } // End of "CorEDecayIn2" |
---|
| 773 | |
---|
| 774 | // Decay of the hadron in 3 particles i=>r+s+t |
---|
| 775 | G4bool G4QHadron::DecayIn3 |
---|
| 776 | (G4LorentzVector& f4Mom, G4LorentzVector& s4Mom, G4LorentzVector& t4Mom) |
---|
| 777 | {// ==================================================================================== |
---|
| 778 | #ifdef debug |
---|
| 779 | G4cout<<"G4QH::DIn3:"<<theMomentum<<"=>pf="<<f4Mom<<"+ps="<<s4Mom<<"+pt="<<t4Mom<<G4endl; |
---|
| 780 | #endif |
---|
| 781 | G4double iM = theMomentum.m(); // Mass of the decaying hadron |
---|
| 782 | G4double fM = f4Mom.m(); // Mass of the 1st hadron |
---|
| 783 | G4double sM = s4Mom.m(); // Mass of the 2nd hadron |
---|
| 784 | G4double tM = t4Mom.m(); // Mass of the 3rd hadron |
---|
| 785 | G4double eps = 0.001; // Accuracy of the split condition |
---|
| 786 | if (fabs(iM-fM-sM-tM)<=eps) |
---|
| 787 | { |
---|
| 788 | G4double fR=fM/iM; |
---|
| 789 | G4double sR=sM/iM; |
---|
| 790 | G4double tR=tM/iM; |
---|
| 791 | f4Mom=fR*theMomentum; |
---|
| 792 | s4Mom=sR*theMomentum; |
---|
| 793 | t4Mom=tR*theMomentum; |
---|
| 794 | return true; |
---|
| 795 | } |
---|
| 796 | if (iM+eps<fM+sM+tM) |
---|
| 797 | { |
---|
| 798 | G4cout<<"***G4QHadron::DecayIn3:fM="<<fM<<" + sM="<<sM<<" + tM="<<tM<<" > iM="<<iM |
---|
| 799 | <<",d="<<iM-fM-sM-tM<<G4endl; |
---|
| 800 | return false; |
---|
| 801 | } |
---|
| 802 | G4double fM2 = fM*fM; |
---|
| 803 | G4double sM2 = sM*sM; |
---|
| 804 | G4double tM2 = tM*tM; |
---|
| 805 | G4double iM2 = iM*iM; |
---|
| 806 | G4double m13sBase=(iM-sM)*(iM-sM)-(fM+tM)*(fM+tM); |
---|
| 807 | G4double m12sMin =(fM+sM)*(fM+sM); |
---|
| 808 | G4double m12sBase=(iM-tM)*(iM-tM)-m12sMin; |
---|
| 809 | G4double rR = 0.; |
---|
| 810 | G4double rnd= 1.; |
---|
| 811 | #ifdef debug |
---|
| 812 | G4int tr = 0; //@@ Comment if "cout" below is skiped @@ |
---|
| 813 | #endif |
---|
| 814 | G4double m12s = 0.; // Fake definition before the Loop |
---|
| 815 | while (rnd > rR) |
---|
| 816 | { |
---|
| 817 | m12s = m12sMin + m12sBase*G4UniformRand(); |
---|
| 818 | G4double e1=m12s+fM2-sM2; |
---|
| 819 | G4double e2=iM2-m12s-tM2; |
---|
| 820 | G4double four12=4*m12s; |
---|
| 821 | G4double m13sRange=0.; |
---|
| 822 | G4double dif=(e1*e1-four12*fM2)*(e2*e2-four12*tM2); |
---|
| 823 | if(dif<0.) |
---|
| 824 | { |
---|
| 825 | #ifdef debug |
---|
| 826 | if(dif<-.01) G4cerr<<"*G4QHadron::DecayIn3:iM="<<iM<<",tM="<<tM<<",sM="<<sM<<",fM=" |
---|
| 827 | <<fM<<",m12(s+f)="<<sqrt(m12s)<<", d="<<iM-fM-sM-tM<<G4endl; |
---|
| 828 | #endif |
---|
| 829 | } |
---|
| 830 | else m13sRange=sqrt(dif)/m12s; |
---|
| 831 | rR = m13sRange/m13sBase; |
---|
| 832 | rnd= G4UniformRand(); |
---|
| 833 | #ifdef debug |
---|
| 834 | G4cout<<"G4QHadron::DecayIn3: try to decay #"<<++tr<<", rR="<<rR<<",rnd="<<rnd<<G4endl; |
---|
| 835 | #endif |
---|
| 836 | } |
---|
| 837 | G4double m12 = sqrt(m12s); // Mass of the H1+H2 system |
---|
| 838 | G4LorentzVector dh4Mom(0.,0.,0.,m12); |
---|
| 839 | |
---|
| 840 | if(!DecayIn2(t4Mom,dh4Mom)) |
---|
| 841 | { |
---|
| 842 | G4cerr<<"***G4QHadron::DecayIn3: Exception1"<<G4endl; |
---|
| 843 | //throw G4QException("G4QHadron::DecayIn3(): DecayIn2 did not succeed"); |
---|
| 844 | return false; |
---|
| 845 | } |
---|
| 846 | #ifdef debug |
---|
| 847 | G4cout<<"G4QHadron::DecayIn3: Now the last decay of m12="<<dh4Mom.m()<<G4endl; |
---|
| 848 | #endif |
---|
| 849 | if(!G4QHadron(dh4Mom).DecayIn2(f4Mom,s4Mom)) |
---|
| 850 | { |
---|
| 851 | G4cerr<<"***G4QHadron::DecayIn3: Error in DecayIn2 -> Exception2"<<G4endl; |
---|
| 852 | //throw G4QException("G4QHadron::DecayIn3(): DecayIn2 did not succeed"); |
---|
| 853 | return false; |
---|
| 854 | } |
---|
| 855 | return true; |
---|
| 856 | } // End of DecayIn3 |
---|
| 857 | |
---|
| 858 | // Relative Decay of the hadron in 3 particles i=>f+s+t (t is with respect to minC<ct<maxC) |
---|
| 859 | G4bool G4QHadron::RelDecayIn3(G4LorentzVector& f4Mom, G4LorentzVector& s4Mom, |
---|
| 860 | G4LorentzVector& t4Mom, G4LorentzVector& dir, |
---|
| 861 | G4double maxCost, G4double minCost) |
---|
| 862 | {// ==================================================================================== |
---|
| 863 | #ifdef debug |
---|
| 864 | G4cout<<"G4QH::RelDIn3:"<<theMomentum<<"=>f="<<f4Mom<<"+s="<<s4Mom<<"+t="<<t4Mom<<G4endl; |
---|
| 865 | #endif |
---|
| 866 | G4double iM = theMomentum.m(); // Mass of the decaying hadron |
---|
| 867 | G4double fM = f4Mom.m(); // Mass of the 1st hadron |
---|
| 868 | G4double sM = s4Mom.m(); // Mass of the 2nd hadron |
---|
| 869 | G4double tM = t4Mom.m(); // Mass of the 3rd hadron |
---|
| 870 | G4double eps = 0.001; // Accuracy of the split condition |
---|
| 871 | if (fabs(iM-fM-sM-tM)<=eps) |
---|
| 872 | { |
---|
| 873 | G4double fR=fM/iM; |
---|
| 874 | G4double sR=sM/iM; |
---|
| 875 | G4double tR=tM/iM; |
---|
| 876 | f4Mom=fR*theMomentum; |
---|
| 877 | s4Mom=sR*theMomentum; |
---|
| 878 | t4Mom=tR*theMomentum; |
---|
| 879 | return true; |
---|
| 880 | } |
---|
| 881 | if (iM+eps<fM+sM+tM) |
---|
| 882 | { |
---|
| 883 | G4cout<<"***G4QHadron::RelDecayIn3:fM="<<fM<<" + sM="<<sM<<" + tM="<<tM<<" > iM="<<iM |
---|
| 884 | <<",d="<<iM-fM-sM-tM<<G4endl; |
---|
| 885 | return false; |
---|
| 886 | } |
---|
| 887 | G4double fM2 = fM*fM; |
---|
| 888 | G4double sM2 = sM*sM; |
---|
| 889 | G4double tM2 = tM*tM; |
---|
| 890 | G4double iM2 = iM*iM; |
---|
| 891 | G4double m13sBase=(iM-sM)*(iM-sM)-(fM+tM)*(fM+tM); |
---|
| 892 | G4double m12sMin =(fM+sM)*(fM+sM); |
---|
| 893 | G4double m12sBase=(iM-tM)*(iM-tM)-m12sMin; |
---|
| 894 | G4double rR = 0.; |
---|
| 895 | G4double rnd= 1.; |
---|
| 896 | #ifdef debug |
---|
| 897 | G4int tr = 0; //@@ Comment if "cout" below is skiped @@ |
---|
| 898 | #endif |
---|
| 899 | G4double m12s = 0.; // Fake definition before the Loop |
---|
| 900 | while (rnd > rR) |
---|
| 901 | { |
---|
| 902 | m12s = m12sMin + m12sBase*G4UniformRand(); |
---|
| 903 | G4double e1=m12s+fM2-sM2; |
---|
| 904 | G4double e2=iM2-m12s-tM2; |
---|
| 905 | G4double four12=4*m12s; |
---|
| 906 | G4double m13sRange=0.; |
---|
| 907 | G4double dif=(e1*e1-four12*fM2)*(e2*e2-four12*tM2); |
---|
| 908 | if(dif<0.) |
---|
| 909 | { |
---|
| 910 | #ifdef debug |
---|
| 911 | if(dif<-.01) G4cerr<<"G4QHadron::RelDecayIn3:iM="<<iM<<",tM="<<tM<<",sM="<<sM<<",fM=" |
---|
| 912 | <<fM<<",m12(s+f)="<<sqrt(m12s)<<", d="<<iM-fM-sM-tM<<G4endl; |
---|
| 913 | #endif |
---|
| 914 | } |
---|
| 915 | else m13sRange=sqrt(dif)/m12s; |
---|
| 916 | rR = m13sRange/m13sBase; |
---|
| 917 | rnd= G4UniformRand(); |
---|
| 918 | #ifdef debug |
---|
| 919 | G4cout<<"G4QHadron::RelDecayIn3: try decay #"<<++tr<<", rR="<<rR<<",rnd="<<rnd<<G4endl; |
---|
| 920 | #endif |
---|
| 921 | } |
---|
| 922 | G4double m12 = sqrt(m12s); // Mass of the H1+H2 system |
---|
| 923 | G4LorentzVector dh4Mom(0.,0.,0.,m12); |
---|
| 924 | |
---|
| 925 | if(!RelDecayIn2(t4Mom,dh4Mom,dir,maxCost,minCost)) |
---|
| 926 | { |
---|
| 927 | G4cerr<<"***G4QHadron::RelDecayIn3: Exception1"<<G4endl; |
---|
| 928 | //throw G4QException("G4QHadron::DecayIn3(): DecayIn2 did not succeed"); |
---|
| 929 | return false; |
---|
| 930 | } |
---|
| 931 | #ifdef debug |
---|
| 932 | G4cout<<"G4QHadron::RelDecayIn3: Now the last decay of m12="<<dh4Mom.m()<<G4endl; |
---|
| 933 | #endif |
---|
| 934 | if(!G4QHadron(dh4Mom).DecayIn2(f4Mom,s4Mom)) |
---|
| 935 | { |
---|
| 936 | G4cerr<<"***G4QHadron::RelDecayIn3: Error in DecayIn2 -> Exception2"<<G4endl; |
---|
| 937 | //throw G4QException("G4QHadron::DecayIn3(): DecayIn2 did not succeed"); |
---|
| 938 | return false; |
---|
| 939 | } |
---|
| 940 | return true; |
---|
| 941 | } // End of RelDecayIn3 |
---|
| 942 | |
---|
| 943 | // Relative Decay of hadron in 3: i=>f+s+t. dN/dO [cp>0:~cost^cp, cp<0:~(1-cost)^(-cp)] |
---|
| 944 | G4bool G4QHadron::CopDecayIn3(G4LorentzVector& f4Mom, G4LorentzVector& s4Mom, |
---|
| 945 | G4LorentzVector& t4Mom, G4LorentzVector& dir, G4double cosp) |
---|
| 946 | {// ==================================================================================== |
---|
| 947 | #ifdef debug |
---|
| 948 | G4cout<<"G4QH::CopDIn3:"<<theMomentum<<"=>f="<<f4Mom<<"+s="<<s4Mom<<"+t="<<t4Mom<<G4endl; |
---|
| 949 | #endif |
---|
| 950 | G4double iM = theMomentum.m(); // Mass of the decaying hadron |
---|
| 951 | G4double fM = f4Mom.m(); // Mass of the 1st hadron |
---|
| 952 | G4double sM = s4Mom.m(); // Mass of the 2nd hadron |
---|
| 953 | G4double tM = t4Mom.m(); // Mass of the 3rd hadron |
---|
| 954 | G4double eps = 0.001; // Accuracy of the split condition |
---|
| 955 | if (fabs(iM-fM-sM-tM)<=eps) |
---|
| 956 | { |
---|
| 957 | G4double fR=fM/iM; |
---|
| 958 | G4double sR=sM/iM; |
---|
| 959 | G4double tR=tM/iM; |
---|
| 960 | f4Mom=fR*theMomentum; |
---|
| 961 | s4Mom=sR*theMomentum; |
---|
| 962 | t4Mom=tR*theMomentum; |
---|
| 963 | return true; |
---|
| 964 | } |
---|
| 965 | if (iM+eps<fM+sM+tM) |
---|
| 966 | { |
---|
| 967 | G4cout<<"***G4QHadron::CopDecayIn3:fM="<<fM<<" + sM="<<sM<<" + tM="<<tM<<" > iM="<<iM |
---|
| 968 | <<",d="<<iM-fM-sM-tM<<G4endl; |
---|
| 969 | return false; |
---|
| 970 | } |
---|
| 971 | G4double fM2 = fM*fM; |
---|
| 972 | G4double sM2 = sM*sM; |
---|
| 973 | G4double tM2 = tM*tM; |
---|
| 974 | G4double iM2 = iM*iM; |
---|
| 975 | G4double m13sBase=(iM-sM)*(iM-sM)-(fM+tM)*(fM+tM); |
---|
| 976 | G4double m12sMin =(fM+sM)*(fM+sM); |
---|
| 977 | G4double m12sBase=(iM-tM)*(iM-tM)-m12sMin; |
---|
| 978 | G4double rR = 0.; |
---|
| 979 | G4double rnd= 1.; |
---|
| 980 | #ifdef debug |
---|
| 981 | G4int tr = 0; //@@ Comment if "cout" below is skiped @@ |
---|
| 982 | #endif |
---|
| 983 | G4double m12s = 0.; // Fake definition before the Loop |
---|
| 984 | while (rnd > rR) |
---|
| 985 | { |
---|
| 986 | m12s = m12sMin + m12sBase*G4UniformRand(); |
---|
| 987 | G4double e1=m12s+fM2-sM2; |
---|
| 988 | G4double e2=iM2-m12s-tM2; |
---|
| 989 | G4double four12=4*m12s; |
---|
| 990 | G4double m13sRange=0.; |
---|
| 991 | G4double dif=(e1*e1-four12*fM2)*(e2*e2-four12*tM2); |
---|
| 992 | if(dif<0.) |
---|
| 993 | { |
---|
| 994 | #ifdef debug |
---|
| 995 | if(dif<-.01) G4cerr<<"G4QHadron::CopDecayIn3:iM="<<iM<<",tM="<<tM<<",sM="<<sM<<",fM=" |
---|
| 996 | <<fM<<",m12(s+f)="<<sqrt(m12s)<<", d="<<iM-fM-sM-tM<<G4endl; |
---|
| 997 | #endif |
---|
| 998 | } |
---|
| 999 | else m13sRange=sqrt(dif)/m12s; |
---|
| 1000 | rR = m13sRange/m13sBase; |
---|
| 1001 | rnd= G4UniformRand(); |
---|
| 1002 | #ifdef debug |
---|
| 1003 | G4cout<<"G4QHadron::CopDecayIn3: try decay #"<<++tr<<", rR="<<rR<<",rnd="<<rnd<<G4endl; |
---|
| 1004 | #endif |
---|
| 1005 | } |
---|
| 1006 | G4double m12 = sqrt(m12s); // Mass of the H1+H2 system |
---|
| 1007 | G4LorentzVector dh4Mom(0.,0.,0.,m12); |
---|
| 1008 | |
---|
| 1009 | if(!CopDecayIn2(t4Mom,dh4Mom,dir,cosp)) |
---|
| 1010 | { |
---|
| 1011 | G4cerr<<"***G4QHadron::CopDecayIn3: Exception1"<<G4endl; |
---|
| 1012 | //throw G4QException("G4QHadron::DecayIn3(): DecayIn2 did not succeed"); |
---|
| 1013 | return false; |
---|
| 1014 | } |
---|
| 1015 | #ifdef debug |
---|
| 1016 | G4cout<<"G4QHadron::DecayIn3: Now the last decay of m12="<<dh4Mom.m()<<G4endl; |
---|
| 1017 | #endif |
---|
| 1018 | if(!G4QHadron(dh4Mom).DecayIn2(f4Mom,s4Mom)) |
---|
| 1019 | { |
---|
| 1020 | G4cerr<<"***G4QHadron::CopDecayIn3: Error in DecayIn2 -> Exception2"<<G4endl; |
---|
| 1021 | //throw G4QException("G4QHadron::DecayIn3(): DecayIn2 did not succeed"); |
---|
| 1022 | return false; |
---|
| 1023 | } |
---|
| 1024 | return true; |
---|
| 1025 | } // End of CopDecayIn3 |
---|
| 1026 | |
---|
| 1027 | // Randomize particle mass taking into account the width |
---|
| 1028 | G4double G4QHadron::RandomizeMass(G4QParticle* pPart, G4double maxM) |
---|
| 1029 | // =========================================================== |
---|
| 1030 | { |
---|
| 1031 | G4double meanM = theQPDG.GetMass(); |
---|
| 1032 | G4double width = theQPDG.GetWidth()/2.; |
---|
| 1033 | #ifdef debug |
---|
| 1034 | G4cout<<"G4QHadron::RandomizeMass: meanM="<<meanM<<", halfWidth="<<width<<G4endl; |
---|
| 1035 | #endif |
---|
| 1036 | if(maxM<meanM-3*width) |
---|
| 1037 | { |
---|
| 1038 | #ifdef debug |
---|
| 1039 | G4cout<<"***G4QH::RandM:m=0 maxM="<<maxM<<"<meanM="<<meanM<<"-3*halfW="<<width<<G4endl; |
---|
| 1040 | #endif |
---|
| 1041 | return 0.; |
---|
| 1042 | } |
---|
| 1043 | ///////////////G4double theMass = 0.; |
---|
| 1044 | if(width==0.) |
---|
| 1045 | { |
---|
| 1046 | #ifdef debug |
---|
| 1047 | if(meanM>maxM) G4cerr<<"***G4QHadron::RandM:Stable m="<<meanM<<">maxM="<<maxM<<G4endl; |
---|
| 1048 | #endif |
---|
| 1049 | return meanM; |
---|
| 1050 | //return 0.; |
---|
| 1051 | } |
---|
| 1052 | else if(width<0.) |
---|
| 1053 | { |
---|
| 1054 | G4cerr<<"***G4QHadron::RandM: width="<<width<<"<0,PDGC="<<theQPDG.GetPDGCode()<<G4endl; |
---|
| 1055 | throw G4QException("G4QHadron::RandomizeMass: with the width of the Hadron < 0."); |
---|
| 1056 | } |
---|
| 1057 | G4double minM = pPart->MinMassOfFragm(); |
---|
| 1058 | if(minM>maxM) |
---|
| 1059 | { |
---|
| 1060 | #ifdef debug |
---|
| 1061 | G4cout<<"***G4QHadron::RandomizeMass:for PDG="<<theQPDG.GetPDGCode()<<" minM="<<minM |
---|
| 1062 | <<" > maxM="<<maxM<<G4endl; |
---|
| 1063 | #endif |
---|
| 1064 | return 0.; |
---|
| 1065 | } |
---|
| 1066 | //Now calculate the Breit-Wigner distribution with two cuts |
---|
| 1067 | G4double v1=atan((minM-meanM)/width); |
---|
| 1068 | G4double v2=atan((maxM-meanM)/width); |
---|
| 1069 | G4double dv=v2-v1; |
---|
| 1070 | #ifdef debug |
---|
| 1071 | G4cout<<"G4QHadr::RandM:Mi="<<minM<<",i="<<v1<<",Ma="<<maxM<<",a="<<v2<<","<<dv<<G4endl; |
---|
| 1072 | #endif |
---|
| 1073 | return meanM+width*tan(v1+dv*G4UniformRand()); |
---|
| 1074 | } |
---|
| 1075 | |
---|
| 1076 | // Split hadron in partons |
---|
| 1077 | void G4QHadron::SplitUp() |
---|
| 1078 | { |
---|
| 1079 | if (IsSplit()) return; |
---|
| 1080 | Splitting(); |
---|
| 1081 | if (Color.empty()) return; |
---|
| 1082 | if (GetSoftCollisionCount() == 0) |
---|
| 1083 | { |
---|
| 1084 | // Diffractive splitting: take the particle definition and get the partons |
---|
| 1085 | G4QParton* Left = 0; |
---|
| 1086 | G4QParton* Right = 0; |
---|
| 1087 | GetValenceQuarkFlavors(Left, Right); |
---|
| 1088 | Left->SetPosition(GetPosition()); |
---|
| 1089 | Right->SetPosition(GetPosition()); |
---|
| 1090 | |
---|
| 1091 | G4LorentzVector HadronMom = Get4Momentum(); |
---|
| 1092 | //G4cout<<"DSU 1 - "<<HadronMom<<G4endl; |
---|
| 1093 | |
---|
| 1094 | // momenta of string ends |
---|
| 1095 | G4double pt2 = HadronMom.perp2(); |
---|
| 1096 | G4double transverseMass2 = HadronMom.plus()*HadronMom.minus(); |
---|
| 1097 | G4double maxAvailMomentum2 = sqr(std::sqrt(transverseMass2) - std::sqrt(pt2)); |
---|
| 1098 | G4ThreeVector pt(minTransverseMass, minTransverseMass, 0); |
---|
| 1099 | if(maxAvailMomentum2/widthOfPtSquare>0.01) |
---|
| 1100 | pt=GaussianPt(widthOfPtSquare, maxAvailMomentum2); |
---|
| 1101 | //G4cout<<"DSU 1.1 - "<<maxAvailMomentum2<<", pt="<<pt<<G4endl; |
---|
| 1102 | |
---|
| 1103 | G4LorentzVector LeftMom(pt, 0.); |
---|
| 1104 | G4LorentzVector RightMom; |
---|
| 1105 | RightMom.setPx(HadronMom.px() - pt.x()); |
---|
| 1106 | RightMom.setPy(HadronMom.py() - pt.y()); |
---|
| 1107 | //G4cout<<"DSU 2: Right4M="<<RightMom<<", Left4M= "<<LeftMom<<G4endl; |
---|
| 1108 | |
---|
| 1109 | G4double Local1 = HadronMom.minus() + (RightMom.perp2() - LeftMom.perp2())/HadronMom.plus(); |
---|
| 1110 | G4double Local2 = std::sqrt(std::max(0., sqr(Local1) - 4.*RightMom.perp2()*HadronMom.minus()/HadronMom.plus())); |
---|
| 1111 | //G4cout<<"DSU 3: L1="<< Local1 <<", L2="<<Local2<<G4endl; |
---|
| 1112 | |
---|
| 1113 | if (Direction) Local2 = -Local2; |
---|
| 1114 | G4double RightMinus = 0.5*(Local1 + Local2); |
---|
| 1115 | G4double LeftMinus = HadronMom.minus() - RightMinus; |
---|
| 1116 | //G4cout<<"DSU 4: Rm="<<RightMinus<<", Lm="<<LeftMinus<<" "<<HadronMom.minus()<<G4endl; |
---|
| 1117 | |
---|
| 1118 | G4double LeftPlus = LeftMom.perp2()/LeftMinus; |
---|
| 1119 | G4double RightPlus = HadronMom.plus() - LeftPlus; |
---|
| 1120 | //G4cout<<"DSU 5: Rp="<<RightPlus<<", Lp="<<LeftPlus<<G4endl; |
---|
| 1121 | LeftMom.setPz(0.5*(LeftPlus - LeftMinus)); |
---|
| 1122 | LeftMom.setE (0.5*(LeftPlus + LeftMinus)); |
---|
| 1123 | RightMom.setPz(0.5*(RightPlus - RightMinus)); |
---|
| 1124 | RightMom.setE (0.5*(RightPlus + RightMinus)); |
---|
| 1125 | //G4cout<<"DSU 6: Left4M="<<LeftMom<<", Right4M="<<RightMom<<G4endl; |
---|
| 1126 | Left->Set4Momentum(LeftMom); |
---|
| 1127 | Right->Set4Momentum(RightMom); |
---|
| 1128 | Color.push_back(Left); |
---|
| 1129 | AntiColor.push_back(Right); |
---|
| 1130 | } |
---|
| 1131 | else |
---|
| 1132 | { |
---|
| 1133 | // Soft hadronization splitting: sample transversal momenta for sea and valence quarks |
---|
| 1134 | G4double phi, pts; |
---|
| 1135 | G4double SumPy = 0.; |
---|
| 1136 | G4double SumPx = 0.; |
---|
| 1137 | G4ThreeVector Pos = GetPosition(); |
---|
| 1138 | G4int nSeaPair = GetSoftCollisionCount()-1; |
---|
| 1139 | |
---|
| 1140 | // here the condition,to ensure viability of splitting, also in cases |
---|
| 1141 | // where difractive excitation occured together with soft scattering. |
---|
| 1142 | // G4double LightConeMomentum = (Direction)? Get4Momentum().plus() : Get4Momentum().minus(); |
---|
| 1143 | // G4double Xmin = theMinPz/LightConeMomentum; |
---|
| 1144 | G4double Xmin = theMinPz/( Get4Momentum().e() - GetMass() ); |
---|
| 1145 | while(Xmin>=1-(2*nSeaPair+1)*Xmin) Xmin*=0.95; |
---|
| 1146 | |
---|
| 1147 | G4int aSeaPair; |
---|
| 1148 | for (aSeaPair = 0; aSeaPair < nSeaPair; aSeaPair++) |
---|
| 1149 | { |
---|
| 1150 | // choose quark flavour, d:u:s = 1:1:(1/StrangeSuppress-2) |
---|
| 1151 | G4int aPDGCode = 1 + (G4int)(G4UniformRand()/StrangeSuppress); |
---|
| 1152 | |
---|
| 1153 | // BuildSeaQuark() determines quark spin, isospin and colour |
---|
| 1154 | // via parton-constructor G4QParton(aPDGCode) |
---|
| 1155 | G4QParton* aParton = BuildSeaQuark(false, aPDGCode); |
---|
| 1156 | |
---|
| 1157 | // G4cout << "G4QGSMSplitableHadron::SoftSplitUp()" << G4endl; |
---|
| 1158 | |
---|
| 1159 | // G4cout << "Parton 1: " |
---|
| 1160 | // << " PDGcode: " << aPDGCode |
---|
| 1161 | // << " - Name: " << aParton->GetDefinition()->GetParticleName() |
---|
| 1162 | // << " - Type: " << aParton->GetDefinition()->GetParticleType() |
---|
| 1163 | // << " - Spin-3: " << aParton->GetSpinZ() |
---|
| 1164 | // << " - Colour: " << aParton->GetColour() << G4endl; |
---|
| 1165 | |
---|
| 1166 | // save colour a spin-3 for anti-quark |
---|
| 1167 | G4int firstPartonColour = aParton->GetColour(); |
---|
| 1168 | G4double firstPartonSpinZ = aParton->GetSpinZ(); |
---|
| 1169 | |
---|
| 1170 | SumPx += aParton->Get4Momentum().px(); |
---|
| 1171 | SumPy += aParton->Get4Momentum().py(); |
---|
| 1172 | Color.push_back(aParton); |
---|
| 1173 | |
---|
| 1174 | // create anti-quark |
---|
| 1175 | aParton = BuildSeaQuark(true, aPDGCode); |
---|
| 1176 | aParton->SetSpinZ(-firstPartonSpinZ); |
---|
| 1177 | aParton->SetColour(-firstPartonColour); |
---|
| 1178 | |
---|
| 1179 | // G4cout << "Parton 2: " |
---|
| 1180 | // << " PDGcode: " << -aPDGCode |
---|
| 1181 | // << " - Name: " << aParton->GetDefinition()->GetParticleName() |
---|
| 1182 | // << " - Type: " << aParton->GetDefinition()->GetParticleType() |
---|
| 1183 | // << " - Spin-3: " << aParton->GetSpinZ() |
---|
| 1184 | // << " - Colour: " << aParton->GetColour() << G4endl; |
---|
| 1185 | // G4cerr << "------------" << G4endl; |
---|
| 1186 | |
---|
| 1187 | SumPx += aParton->Get4Momentum().px(); |
---|
| 1188 | SumPy += aParton->Get4Momentum().py(); |
---|
| 1189 | AntiColor.push_back(aParton); |
---|
| 1190 | } |
---|
| 1191 | // Valence quark |
---|
| 1192 | G4QParton* pColorParton = 0; |
---|
| 1193 | G4QParton* pAntiColorParton = 0; |
---|
| 1194 | GetValenceQuarkFlavors(pColorParton, pAntiColorParton); |
---|
| 1195 | G4int ColorEncoding = pColorParton->GetPDGCode(); |
---|
| 1196 | G4int AntiColorEncoding = pAntiColorParton->GetPDGCode(); |
---|
| 1197 | |
---|
| 1198 | pts = sigmaPt*std::sqrt(-std::log(G4UniformRand())); |
---|
| 1199 | phi = twopi*G4UniformRand(); |
---|
| 1200 | G4double Px = pts*std::cos(phi); |
---|
| 1201 | G4double Py = pts*std::sin(phi); |
---|
| 1202 | SumPx += Px; |
---|
| 1203 | SumPy += Py; |
---|
| 1204 | |
---|
| 1205 | if (ColorEncoding < 0) // use particle definition |
---|
| 1206 | { |
---|
| 1207 | G4LorentzVector ColorMom(-SumPx, -SumPy, 0, 0); |
---|
| 1208 | pColorParton->Set4Momentum(ColorMom); |
---|
| 1209 | G4LorentzVector AntiColorMom(Px, Py, 0, 0); |
---|
| 1210 | pAntiColorParton->Set4Momentum(AntiColorMom); |
---|
| 1211 | } |
---|
| 1212 | else |
---|
| 1213 | { |
---|
| 1214 | G4LorentzVector ColorMom(Px, Py, 0, 0); |
---|
| 1215 | pColorParton->Set4Momentum(ColorMom); |
---|
| 1216 | G4LorentzVector AntiColorMom(-SumPx, -SumPy, 0, 0); |
---|
| 1217 | pAntiColorParton->Set4Momentum(AntiColorMom); |
---|
| 1218 | } |
---|
| 1219 | Color.push_back(pColorParton); |
---|
| 1220 | AntiColor.push_back(pAntiColorParton); |
---|
| 1221 | |
---|
| 1222 | // Sample X |
---|
| 1223 | G4int nAttempt = 0; |
---|
| 1224 | G4double SumX = 0; |
---|
| 1225 | G4double aBeta = beta; |
---|
| 1226 | G4double ColorX, AntiColorX; |
---|
| 1227 | G4double HPWtest = 0; |
---|
| 1228 | G4int aPDG=std::abs(GetPDGCode()); |
---|
| 1229 | if (aPDG ==211 || aPDG == 22 || aPDG == 111) aBeta = 1.; |
---|
| 1230 | else if (aPDG == 321) aBeta = 0.; |
---|
| 1231 | else G4cout<<"-Warning-G4QHadron::SplitUp: wrong PDG="<<GetPDGCode()<<G4endl; |
---|
| 1232 | do |
---|
| 1233 | { |
---|
| 1234 | SumX = 0; |
---|
| 1235 | nAttempt++; |
---|
| 1236 | G4int NumberOfUnsampledSeaQuarks = 2*nSeaPair; |
---|
| 1237 | G4double beta1 = beta; |
---|
| 1238 | if (std::abs(ColorEncoding) <= 1000 && std::abs(AntiColorEncoding) <= 1000) beta1 = 1.; //... in a meson |
---|
| 1239 | ColorX = SampleX(Xmin, NumberOfUnsampledSeaQuarks, 2*nSeaPair, aBeta); |
---|
| 1240 | HPWtest = ColorX; |
---|
| 1241 | while (ColorX < Xmin || ColorX > 1.|| 1. - ColorX <= Xmin) { |
---|
| 1242 | ; // Possible dead loop? Don't know why this loop is here - DHW |
---|
| 1243 | } |
---|
| 1244 | Color.back()->SetX(SumX = ColorX);// this is the valenz quark. |
---|
| 1245 | |
---|
| 1246 | std::list<G4QParton*>::iterator icolor = Color.begin(); |
---|
| 1247 | std::list<G4QParton*>::iterator ecolor = Color.end(); |
---|
| 1248 | std::list<G4QParton*>::iterator ianticolor = AntiColor.begin(); |
---|
| 1249 | std::list<G4QParton*>::iterator eanticolor = AntiColor.end(); |
---|
| 1250 | for ( ; icolor != ecolor && ianticolor != eanticolor; ++icolor, ++ianticolor) |
---|
| 1251 | { |
---|
| 1252 | NumberOfUnsampledSeaQuarks--; |
---|
| 1253 | ColorX = SampleX(Xmin, NumberOfUnsampledSeaQuarks, 2*nSeaPair, aBeta); |
---|
| 1254 | (*icolor)->SetX(ColorX); |
---|
| 1255 | SumX += ColorX; |
---|
| 1256 | NumberOfUnsampledSeaQuarks--; |
---|
| 1257 | AntiColorX = SampleX(Xmin, NumberOfUnsampledSeaQuarks, 2*nSeaPair, aBeta); |
---|
| 1258 | (*ianticolor)->SetX(AntiColorX); // the 'sea' partons |
---|
| 1259 | SumX += AntiColorX; |
---|
| 1260 | if (1. - SumX <= Xmin) break; |
---|
| 1261 | } |
---|
| 1262 | } while (1. - SumX <= Xmin); |
---|
| 1263 | AntiColor.back()->SetX(1.0 - SumX); // the di-quark takes the rest, then go to momentum |
---|
| 1264 | // and here is the bug ;-) @@@@@@@@@@@@@ |
---|
| 1265 | if(getenv("debug_QGSMSplitableHadron") )G4cout << "particle energy at split = "<<Get4Momentum().t()<<G4endl; |
---|
| 1266 | G4double lightCone = ((!Direction) ? Get4Momentum().minus() : Get4Momentum().plus()); |
---|
| 1267 | // lightCone -= 0.5*Get4Momentum().m(); |
---|
| 1268 | // hpw testing @@@@@ lightCone = 2.*Get4Momentum().t(); |
---|
| 1269 | if(getenv("debug_QGSMSplitableHadron") )G4cout << "Light cone = "<<lightCone<<G4endl; |
---|
| 1270 | std::list<G4QParton*>::iterator icolor = Color.begin(); |
---|
| 1271 | std::list<G4QParton*>::iterator ecolor = Color.end(); |
---|
| 1272 | std::list<G4QParton*>::iterator ianticolor = AntiColor.begin(); |
---|
| 1273 | std::list<G4QParton*>::iterator eanticolor = AntiColor.end(); |
---|
| 1274 | for ( ; icolor != ecolor && ianticolor != eanticolor; ++icolor, ++ianticolor) |
---|
| 1275 | { |
---|
| 1276 | (*icolor)->DefineMomentumInZ(lightCone, Direction); |
---|
| 1277 | (*ianticolor)->DefineMomentumInZ(lightCone, Direction); |
---|
| 1278 | } |
---|
| 1279 | //G4cout <<G4endl<<"XSAMPLE "<<HPWtest<<G4endl; |
---|
| 1280 | return; |
---|
| 1281 | } |
---|
| 1282 | } // End of SplitUp |
---|
| 1283 | |
---|
| 1284 | // Boost hadron 4-momentum, using Boost Lorentz vector |
---|
| 1285 | void G4QHadron::Boost(const G4LorentzVector& boost4M) |
---|
| 1286 | { |
---|
| 1287 | // see CERNLIB short writeup U101 for the algorithm |
---|
| 1288 | G4double bm=boost4M.mag(); |
---|
| 1289 | G4double factor=(theMomentum.vect()*boost4M.vect()/(boost4M.e()+bm)-theMomentum.e())/bm; |
---|
| 1290 | theMomentum.setE(theMomentum.dot(boost4M)/bm); |
---|
| 1291 | theMomentum.setVect(factor*boost4M.vect() + theMomentum.vect()); |
---|
| 1292 | } // End of Boost |
---|
| 1293 | |
---|
| 1294 | // Build one (?M.K.) sea-quark |
---|
| 1295 | G4QParton* G4QHadron::BuildSeaQuark(G4bool isAntiQuark, G4int aPDGCode) |
---|
| 1296 | { |
---|
| 1297 | if (isAntiQuark) aPDGCode*=-1; |
---|
| 1298 | G4QParton* result = new G4QParton(aPDGCode); |
---|
| 1299 | result->SetPosition(GetPosition()); |
---|
| 1300 | G4ThreeVector aPtVector = GaussianPt(sigmaPt, DBL_MAX); |
---|
| 1301 | G4LorentzVector a4Momentum(aPtVector, 0); |
---|
| 1302 | result->Set4Momentum(a4Momentum); |
---|
| 1303 | return result; |
---|
| 1304 | } // End of BuildSeaQuark |
---|
| 1305 | |
---|
| 1306 | G4double G4QHadron::SampleX(G4double anXmin, G4int nSea, G4int totalSea, G4double aBeta) |
---|
| 1307 | { |
---|
| 1308 | G4double result; |
---|
| 1309 | G4double x1, x2; |
---|
| 1310 | G4double ymax = 0; |
---|
| 1311 | for(G4int ii=0; ii<100; ii++) // @@ 100 is hardwired ? M.K. |
---|
| 1312 | { |
---|
| 1313 | G4double y = std::pow(1./G4double(ii), alpha); |
---|
| 1314 | y*=std::pow(std::pow(1.-anXmin-totalSea*anXmin,alpha+1)-std::pow(anXmin,alpha+1),nSea); |
---|
| 1315 | y*=std::pow(1.-anXmin-totalSea*anXmin, aBeta+1) - std::pow(anXmin, aBeta+1); |
---|
| 1316 | if(y>ymax) ymax = y; |
---|
| 1317 | } |
---|
| 1318 | G4double y; |
---|
| 1319 | G4double xMax=1.-(totalSea+1.)*anXmin; |
---|
| 1320 | if(anXmin > xMax) |
---|
| 1321 | { |
---|
| 1322 | G4cerr<<"***G4QHadron::SampleX: anXmin="<<anXmin<<" > xMax="<<xMax<<", nSea="<<nSea |
---|
| 1323 | <<", totSea="<<totalSea<<G4endl; |
---|
| 1324 | G4Exception("G4QHadron::SampleX:","72",FatalException,"TooBigXValue"); |
---|
| 1325 | } |
---|
| 1326 | do |
---|
| 1327 | { |
---|
| 1328 | x1 = CLHEP::RandFlat::shoot(anXmin, xMax); |
---|
| 1329 | y = std::pow(x1, alpha); |
---|
| 1330 | y*=std::pow(std::pow(1.-x1-totalSea*anXmin,alpha+1) - std::pow(anXmin, alpha+1), nSea); |
---|
| 1331 | y*=std::pow(1.-x1-totalSea*anXmin, aBeta+1) - std::pow(anXmin, aBeta+1); |
---|
| 1332 | x2 = ymax*G4UniformRand(); |
---|
| 1333 | } while(x2>y); |
---|
| 1334 | result = x1; |
---|
| 1335 | return result; |
---|
| 1336 | } // End of SampleX |
---|
| 1337 | |
---|
| 1338 | |
---|
| 1339 | void G4QHadron::GetValenceQuarkFlavors(G4QParton* &Parton1, G4QParton* &Parton2) |
---|
| 1340 | { |
---|
| 1341 | // Note! convention aEnd = q or (qq)bar and bEnd = qbar or qq. |
---|
| 1342 | G4int aEnd=0; |
---|
| 1343 | G4int bEnd=0; |
---|
| 1344 | G4int HadronEncoding = GetPDGCode(); |
---|
| 1345 | if(!(GetBaryonNumber())) SplitMeson(HadronEncoding,&aEnd,&bEnd); |
---|
| 1346 | else SplitBaryon(HadronEncoding, &aEnd, &bEnd); |
---|
| 1347 | |
---|
| 1348 | Parton1 = new G4QParton(aEnd); |
---|
| 1349 | Parton1->SetPosition(GetPosition()); |
---|
| 1350 | |
---|
| 1351 | // G4cerr << "G4QGSMSplitableHadron::GetValenceQuarkFlavors()" << G4endl; |
---|
| 1352 | // G4cerr << "Parton 1: " |
---|
| 1353 | // << " PDGcode: " << aEnd |
---|
| 1354 | // << " - Name: " << Parton1->GetDefinition()->GetParticleName() |
---|
| 1355 | // << " - Type: " << Parton1->GetDefinition()->GetParticleType() |
---|
| 1356 | // << " - Spin-3: " << Parton1->GetSpinZ() |
---|
| 1357 | // << " - Colour: " << Parton1->GetColour() << G4endl; |
---|
| 1358 | |
---|
| 1359 | Parton2 = new G4QParton(bEnd); |
---|
| 1360 | Parton2->SetPosition(GetPosition()); |
---|
| 1361 | |
---|
| 1362 | // G4cerr << "Parton 2: " |
---|
| 1363 | // << " PDGcode: " << bEnd |
---|
| 1364 | // << " - Name: " << Parton2->GetDefinition()->GetParticleName() |
---|
| 1365 | // << " - Type: " << Parton2->GetDefinition()->GetParticleType() |
---|
| 1366 | // << " - Spin-3: " << Parton2->GetSpinZ() |
---|
| 1367 | // << " - Colour: " << Parton2->GetColour() << G4endl; |
---|
| 1368 | // G4cerr << "... now checking for color and spin conservation - yielding: " << G4endl; |
---|
| 1369 | |
---|
| 1370 | // colour of parton 1 choosen at random by G4QParton(aEnd) |
---|
| 1371 | // colour of parton 2 is the opposite: |
---|
| 1372 | |
---|
| 1373 | Parton2->SetColour(-(Parton1->GetColour())); |
---|
| 1374 | |
---|
| 1375 | // isospin-3 of both partons is handled by G4Parton(PDGCode) |
---|
| 1376 | |
---|
| 1377 | // spin-3 of parton 1 and 2 choosen at random by G4QParton(aEnd) |
---|
| 1378 | // spin-3 of parton 2 may be constrained by spin of original particle: |
---|
| 1379 | |
---|
| 1380 | if ( std::abs(Parton1->GetSpinZ() + Parton2->GetSpinZ()) > GetSpin()) |
---|
| 1381 | { |
---|
| 1382 | Parton2->SetSpinZ(-(Parton2->GetSpinZ())); |
---|
| 1383 | } |
---|
| 1384 | |
---|
| 1385 | // G4cerr << "Parton 2: " |
---|
| 1386 | // << " PDGcode: " << bEnd |
---|
| 1387 | // << " - Name: " << Parton2->GetDefinition()->GetParticleName() |
---|
| 1388 | // << " - Type: " << Parton2->GetDefinition()->GetParticleType() |
---|
| 1389 | // << " - Spin-3: " << Parton2->GetSpinZ() |
---|
| 1390 | // << " - Colour: " << Parton2->GetColour() << G4endl; |
---|
| 1391 | // G4cerr << "------------" << G4endl; |
---|
| 1392 | |
---|
| 1393 | } // End of GetValenceQuarkFlavors |
---|
| 1394 | |
---|
| 1395 | G4bool G4QHadron::SplitMeson(G4int PDGcode, G4int* aEnd, G4int* bEnd) |
---|
| 1396 | { |
---|
| 1397 | G4bool result = true; |
---|
| 1398 | G4int absPDGcode = std::abs(PDGcode); |
---|
| 1399 | if (absPDGcode >= 1000) return false; |
---|
| 1400 | if(absPDGcode == 22) |
---|
| 1401 | { |
---|
| 1402 | G4int it=1; |
---|
| 1403 | if(G4UniformRand()<.5) it++; |
---|
| 1404 | *aEnd = it; |
---|
| 1405 | *bEnd = -it; |
---|
| 1406 | } |
---|
| 1407 | else |
---|
| 1408 | { |
---|
| 1409 | G4int heavy = absPDGcode/100; |
---|
| 1410 | G4int light = (absPDGcode%100)/10; |
---|
| 1411 | G4int anti = 1 - 2*(std::max(heavy, light)%2); |
---|
| 1412 | if (PDGcode < 0 ) anti = -anti; |
---|
| 1413 | heavy *= anti; |
---|
| 1414 | light *= -anti; |
---|
| 1415 | if ( anti < 0) G4SwapObj(&heavy, &light); |
---|
| 1416 | *aEnd = heavy; |
---|
| 1417 | *bEnd = light; |
---|
| 1418 | } |
---|
| 1419 | return result; |
---|
| 1420 | } |
---|
| 1421 | |
---|
| 1422 | G4bool G4QHadron::SplitBaryon(G4int PDGcode, G4int* quark, G4int* diQuark) |
---|
| 1423 | { |
---|
| 1424 | static const G4double r2=.5; |
---|
| 1425 | static const G4double r3=1./3.; |
---|
| 1426 | static const G4double d3=2./3.; |
---|
| 1427 | static const G4double r4=1./4.; |
---|
| 1428 | static const G4double r6=1./6.; |
---|
| 1429 | static const G4double r12=1./12.; |
---|
| 1430 | // |
---|
| 1431 | std::pair<G4int,G4int> qdq[5]; |
---|
| 1432 | G4double prb[5]; |
---|
| 1433 | G4int nc=0; |
---|
| 1434 | G4int aPDGcode=std::abs(PDGcode); |
---|
| 1435 | if(aPDGcode==2212) // ==> Proton |
---|
| 1436 | { |
---|
| 1437 | nc=3; |
---|
| 1438 | qdq[0]=make_pair(2203, 1); prb[0]=r3; // uu_1, d |
---|
| 1439 | qdq[1]=make_pair(2103, 2); prb[1]=r6; // ud_1, u |
---|
| 1440 | qdq[2]=make_pair(2101, 2); prb[2]=r2; // ud_0, u |
---|
| 1441 | } |
---|
| 1442 | else if(aPDGcode==2112) // ==> Neutron |
---|
| 1443 | { |
---|
| 1444 | nc=3; |
---|
| 1445 | qdq[0]=make_pair(2103, 1); prb[0]=r6; // ud_1, d |
---|
| 1446 | qdq[1]=make_pair(2101, 1); prb[1]=r2; // ud_0, d |
---|
| 1447 | qdq[2]=make_pair(1103, 2); prb[2]=r3; // dd_1, u |
---|
| 1448 | } |
---|
| 1449 | else if(aPDGcode%10<3) // ==> Spin 1/2 Hyperons |
---|
| 1450 | { |
---|
| 1451 | if(aPDGcode==3122) |
---|
| 1452 | { |
---|
| 1453 | nc=5; |
---|
| 1454 | qdq[0]=make_pair(2103, 3); prb[0]=r3; // ud_1, s |
---|
| 1455 | qdq[1]=make_pair(3203, 1); prb[1]=r4; // su_1, d |
---|
| 1456 | qdq[2]=make_pair(3201, 1); prb[2]=r12; // su_0, d |
---|
| 1457 | qdq[3]=make_pair(3103, 2); prb[3]=r4; // sd_1, u |
---|
| 1458 | qdq[4]=make_pair(3101, 2); prb[4]=r12; // sd_0, u |
---|
| 1459 | } |
---|
| 1460 | else if(aPDGcode==3222) |
---|
| 1461 | { |
---|
| 1462 | nc=3; |
---|
| 1463 | qdq[0]=make_pair(2203, 3); prb[0]=r3; |
---|
| 1464 | qdq[1]=make_pair(3203, 2); prb[1]=r6; |
---|
| 1465 | qdq[2]=make_pair(3201, 2); prb[2]=r2; |
---|
| 1466 | } |
---|
| 1467 | else if(aPDGcode==3212) |
---|
| 1468 | { |
---|
| 1469 | nc=5; |
---|
| 1470 | qdq[0]=make_pair(2103, 3); prb[0]=r3; |
---|
| 1471 | qdq[1]=make_pair(3203, 1); prb[1]=r12; |
---|
| 1472 | qdq[2]=make_pair(3201, 1); prb[2]=r4; |
---|
| 1473 | qdq[3]=make_pair(3103, 2); prb[3]=r12; |
---|
| 1474 | qdq[4]=make_pair(3101, 2); prb[4]=r4; |
---|
| 1475 | } |
---|
| 1476 | else if(aPDGcode==3112) |
---|
| 1477 | { |
---|
| 1478 | nc=3; |
---|
| 1479 | qdq[0]=make_pair(1103, 3); prb[0]=r3; |
---|
| 1480 | qdq[1]=make_pair(3103, 1); prb[1]=r6; |
---|
| 1481 | qdq[2]=make_pair(3101, 1); prb[2]=r2; |
---|
| 1482 | } |
---|
| 1483 | else if(aPDGcode==3312) |
---|
| 1484 | { |
---|
| 1485 | nc=3; |
---|
| 1486 | qdq[0]=make_pair(3103, 3); prb[0]=r6; |
---|
| 1487 | qdq[1]=make_pair(3101, 3); prb[1]=r2; |
---|
| 1488 | qdq[2]=make_pair(3303, 1); prb[2]=r3; |
---|
| 1489 | } |
---|
| 1490 | else if(aPDGcode==3322) |
---|
| 1491 | { |
---|
| 1492 | nc=3; |
---|
| 1493 | qdq[0]=make_pair(3203, 3); prb[0]=r6; |
---|
| 1494 | qdq[1]=make_pair(3201, 3); prb[1]=r2; |
---|
| 1495 | qdq[2]=make_pair(3303, 2); prb[2]=r3; |
---|
| 1496 | } |
---|
| 1497 | else return false; |
---|
| 1498 | } |
---|
| 1499 | else // ==> Spin 3/2 Baryons (Spin>3/2 is ERROR) |
---|
| 1500 | { |
---|
| 1501 | if(aPDGcode==3334) |
---|
| 1502 | { |
---|
| 1503 | nc=1; |
---|
| 1504 | qdq[0]=make_pair(3303, 3); prb[0]=1.; |
---|
| 1505 | } |
---|
| 1506 | else if(aPDGcode==2224) |
---|
| 1507 | { |
---|
| 1508 | nc=1; |
---|
| 1509 | qdq[0]=make_pair(2203, 2); prb[0]=1.; |
---|
| 1510 | } |
---|
| 1511 | else if(aPDGcode==2214) |
---|
| 1512 | { |
---|
| 1513 | nc=2; |
---|
| 1514 | qdq[0]=make_pair(2203, 1); prb[0]=r3; |
---|
| 1515 | qdq[1]=make_pair(2103, 2); prb[1]=d3; |
---|
| 1516 | } |
---|
| 1517 | else if(aPDGcode==2114) |
---|
| 1518 | { |
---|
| 1519 | nc=2; |
---|
| 1520 | qdq[0]=make_pair(2103, 1); prb[0]=d3; |
---|
| 1521 | qdq[1]=make_pair(2103, 2); prb[1]=r3; |
---|
| 1522 | } |
---|
| 1523 | else if(aPDGcode==1114) |
---|
| 1524 | { |
---|
| 1525 | nc=1; |
---|
| 1526 | qdq[0]=make_pair(1103, 1); prb[0]=1.; |
---|
| 1527 | } |
---|
| 1528 | else if(aPDGcode==3224) |
---|
| 1529 | { |
---|
| 1530 | nc=2; |
---|
| 1531 | qdq[0]=make_pair(2203, 3); prb[0]=r3; |
---|
| 1532 | qdq[1]=make_pair(3203, 2); prb[1]=d3; |
---|
| 1533 | } |
---|
| 1534 | else if(aPDGcode==3214) |
---|
| 1535 | { |
---|
| 1536 | nc=3; |
---|
| 1537 | qdq[0]=make_pair(2103, 3); prb[0]=r3; |
---|
| 1538 | qdq[1]=make_pair(3203, 1); prb[1]=r3; |
---|
| 1539 | qdq[2]=make_pair(3103, 2); prb[2]=r3; |
---|
| 1540 | } |
---|
| 1541 | else if(aPDGcode==3114) |
---|
| 1542 | { |
---|
| 1543 | nc=2; |
---|
| 1544 | qdq[0]=make_pair(1103, 3); prb[0]=r3; |
---|
| 1545 | qdq[1]=make_pair(3103, 1); prb[1]=d3; |
---|
| 1546 | } |
---|
| 1547 | else if(aPDGcode==3324) |
---|
| 1548 | { |
---|
| 1549 | nc=2; |
---|
| 1550 | qdq[0]=make_pair(3203, 3); prb[0]=r3; |
---|
| 1551 | qdq[1]=make_pair(3303, 2); prb[1]=d3; |
---|
| 1552 | } |
---|
| 1553 | else if(aPDGcode==3314) |
---|
| 1554 | { |
---|
| 1555 | nc=2; |
---|
| 1556 | qdq[0]=make_pair(3103, 3); prb[0]=d3; |
---|
| 1557 | qdq[1]=make_pair(3303, 1); prb[1]=r3; |
---|
| 1558 | } |
---|
| 1559 | else return false; |
---|
| 1560 | } |
---|
| 1561 | G4double random = G4UniformRand(); |
---|
| 1562 | G4double sum = 0.; |
---|
| 1563 | for(G4int i=0; i<nc; i++) |
---|
| 1564 | { |
---|
| 1565 | sum += prb[i]; |
---|
| 1566 | if(sum>random) |
---|
| 1567 | { |
---|
| 1568 | if(PDGcode>0) |
---|
| 1569 | { |
---|
| 1570 | *diQuark= qdq[i].first; |
---|
| 1571 | *quark = qdq[i].second; |
---|
| 1572 | } |
---|
| 1573 | else |
---|
| 1574 | { |
---|
| 1575 | *diQuark= -qdq[i].first; |
---|
| 1576 | *quark = -qdq[i].second; |
---|
| 1577 | } |
---|
| 1578 | break; |
---|
| 1579 | } |
---|
| 1580 | } |
---|
| 1581 | return true; |
---|
| 1582 | } |
---|
| 1583 | |
---|
| 1584 | G4ThreeVector G4QHadron::GaussianPt(G4double widthSquare, G4double maxPtSquare) |
---|
| 1585 | { |
---|
| 1586 | G4double R=0.; |
---|
| 1587 | while ((R = -widthSquare*std::log(G4UniformRand())) > maxPtSquare) { |
---|
| 1588 | ; |
---|
| 1589 | } |
---|
| 1590 | R = std::sqrt(R); |
---|
| 1591 | G4double phi = twopi*G4UniformRand(); |
---|
| 1592 | return G4ThreeVector(R*std::cos(phi), R*std::sin(phi), 0.); |
---|
| 1593 | } |
---|