#include <G4QIonIonCollision.hh>
Public Member Functions | |
G4QIonIonCollision (G4QNucleus &pNucleus, const G4QNucleus &tNucleus) | |
~G4QIonIonCollision () | |
G4QHadronVector * | Fragment () |
Static Public Member Functions | |
static void | SetParameters (G4int nC, G4double strTens, G4double tubeDens, G4double SigPt) |
Protected Member Functions | |
G4bool | ExciteDiffParticipants (G4QHadron *aPartner, G4QHadron *bPartner) const |
G4bool | ExciteSingDiffParticipants (G4QHadron *aPartner, G4QHadron *bPartner) const |
std::pair< G4int, G4int > | ReducePair (G4int P1, G4int P2) const |
void | Breeder () |
G4bool | IsSingleDiffractive () |
G4int | SumPartonPDG (G4int PDG1, G4int PFG2) const |
G4double | ChooseX (G4double Xmin, G4double Xmax) const |
G4ThreeVector | GaussianPt (G4double widthSquare, G4double maxPtSquare) const |
G4int | AnnihilationOrder (G4int LS, G4int MS, G4int uP, G4int mP, G4int sP, G4int nP) |
void | SwapPartons () |
Definition at line 59 of file G4QIonIonCollision.hh.
G4QIonIonCollision::G4QIonIonCollision | ( | G4QNucleus & | pNucleus, | |
const G4QNucleus & | tNucleus | |||
) |
Definition at line 62 of file G4QIonIonCollision.cc.
References G4QString::Boost(), G4QNucleus::ChooseImpactXandY(), DBL_MAX, G4QNucleus::DeleteNucleons(), G4QPartonPair::DIFFRACTIVE, G4QNucleus::DoLorentzBoost(), G4QNucleus::DoLorentzRotation(), ExciteDiffParticipants(), ExciteSingDiffParticipants(), FatalException, G4cerr, G4cout, G4endl, G4Exception(), G4UniformRand, G4QCHIPSWorld::Get(), G4QParton::Get4Momentum(), G4QString::Get4Momentum(), G4QHadron::Get4Momentum(), G4QNucleus::GetA(), G4QHadron::GetAntiColor(), G4QHadron::GetBaryonNumber(), G4QPartonPair::GetCollisionType(), G4QHadron::GetColor(), G4QProbability::GetCutPomProbability(), G4QNucleus::GetGSMass(), G4QHadron::GetMass(), G4QNucleus::GetN(), G4QHadron::GetNextAntiParton(), G4QNucleus::GetNextNucleon(), G4QHadron::GetNextParton(), G4QNucleus::GetNucleons4Momentum(), G4QNucleus::GetOuterRadius(), G4QParton::GetPDGCode(), G4QHadron::GetPDGCode(), G4QProbability::GetPomDiffProbability(), G4QProbability::GetPomInelProbability(), G4QHadron::GetPosition(), G4QHadron::GetQC(), G4QContent::GetSPDGCode(), G4QNucleus::GetThickness(), G4QParton::GetType(), G4QNucleus::GetZ(), G4QHadron::IncrementCollisionCount(), G4QNucleus::Init3D(), G4QNucleus::InitByPDG(), IsSingleDiffractive(), LT, G4QPDGCode::MakeTwoBaryons(), G4QPartonPair::PROJECTILE, G4QParton::ReduceDiQADiQ(), G4QParton::Set4Momentum(), G4QHadron::Set4Momentum(), G4QInteraction::SetNumberOfDiffractiveCollisions(), G4QInteraction::SetNumberOfDINRCollisions(), G4QInteraction::SetNumberOfSoftCollisions(), G4QParton::SetPDGCode(), G4QInteraction::SetTarget(), G4InuclParticleNames::sm, G4QPartonPair::SOFT, sqr(), G4QNucleus::StartLoop(), G4QNucleus::SubtractNucleon(), SwapPartons(), and G4QPartonPair::TARGET.
00063 { 00064 static const G4double mProt= G4QPDGCode(2212).GetMass(); // Mass of proton 00065 static const G4double mPi0= G4QPDGCode(111).GetMass(); // Mass of Pi0 00066 theWorld= G4QCHIPSWorld::Get(); // Get a pointer to the CHIPS World 00067 theResult = new G4QHadronVector; // Define theResultingHadronVector 00068 G4bool stringsInitted=false; // Strings are initiated 00069 G4LorentzRotation toZ; // Lorentz Transformation to the projectileSystem 00070 G4LorentzVector proj4M=pNucleus.Get4Momentum(); // Projectile nucleus 4-momentum in LS 00071 //G4double projM=pNucleus.GetGSMass(); // Ground state mass of the projectile nucleus 00072 G4double targM=tNucleus.GetGSMass(); // Ground state mass of the target nucleus 00073 #ifdef edebug 00074 G4cout<<"G4QIn::Constr:*Called*,pA="<<pNucleus<<proj4M<<",tA="<<tNucleus<<targM<<G4endl; 00075 #endif 00076 G4int pZ=pNucleus.GetZ(); 00077 G4int pN=pNucleus.GetN(); 00078 G4int pA=pZ+pN; 00079 G4int pPDG=90000000+pZ*1000+pN; 00080 G4int tZ=tNucleus.GetZ(); 00081 G4int tN=tNucleus.GetN(); 00082 G4int tA=tZ+tN; 00083 G4int tPDG=90000000+tZ*1000+tN; 00084 toZ.rotateZ(-proj4M.phi()); 00085 toZ.rotateY(-proj4M.theta()); 00086 G4LorentzVector zProj4M=toZ*proj4M; // Proj 4-momentum in LS rotated to Z axis 00087 pNucleus.Set4Momentum(zProj4M); // Now the projectile nucleus moves along Z axis 00088 #ifdef edebug 00089 G4int totChg=pZ+tZ; // Charge of the projectile+target for the CHECK 00090 G4int totBaN=pA+tA; // Baryon Number of Proj+Targ for CHECK 00091 G4LorentzVector tgLS4M(0.,0.,0.,targM); // Target 4-momentum in LS 00092 G4LorentzVector totLS4M=proj4M+tgLS4M; // Total 4-momentum in LS 00093 G4LorentzVector totZLS4M=zProj4M+tgLS4M;// Total 4-momentum in LS with momentum along Z 00094 G4cout<<"-EMC-G4QIonIonCollision::Constr: tLS4M="<<totLS4M<<",tZLS4M="<<totZLS4M<<G4endl; 00095 // === From nere all consideration is made in the rotated LS frame (proj is along Z) === 00096 #endif 00097 G4LorentzRotation toLab(toZ.inverse()); // Lorentz Transfornation "ZLS"->LS (at the end) 00098 G4QInteractionVector theInteractions; // A vector of interactions (taken from the Body) 00099 G4QPartonPairVector thePartonPairs; // The parton pairs (taken from the Body) 00100 G4int ModelMode=SOFT; // The current model type (taken from the Body) 00101 theTargNucleus=G4QNucleus(tZ,tN); // Create theTargNucleus to Move From LS to CM 00102 theTargNucleus.InitByPDG(tPDG); // Reinit the Nucleus for the new Attempt? 00103 theTargNucleus.Init3D(); // 3D-initialisation(nucleons)ofTheTGNucleusClone 00104 #ifdef debug 00105 G4cout<<"G4QIonIonCollision::Constr:TargNuclMom="<<theTargNucleus.Get4Momentum()<<G4endl; 00106 #endif 00107 theProjNucleus=G4QNucleus(pZ,pN); // Create theProjNucleus to Move From LS to CM 00108 theProjNucleus.InitByPDG(pPDG); // Reinit the Nucleus for the new Attempt? 00109 theProjNucleus.Init3D(); // 3D-initialisation(nucleons)ofThePrNucleusClone 00110 #ifdef debug 00111 G4cout<<"G4QIonIonCollision::Constr:ProjNuclMom="<<theProjNucleus.Get4Momentum()<<G4endl; 00112 #endif 00113 #ifdef edebug 00114 G4LorentzVector sumP1=theProjNucleus.GetNucleons4Momentum();// Sum ofPrNucleons4M inRotLS 00115 G4LorentzVector sumT1=theTargNucleus.GetNucleons4Momentum();// Sum ofTgNucleons4M inRotLS 00116 G4cout<<"-EMC-G4QIonIonCollision::Construct: pNuc4M="<<sumP1<<", tNuc4M="<<sumT1<<G4endl; 00117 #endif 00118 G4ThreeVector theCMVelocity(0.,0.,0.); // Prototype of the "CM" velocity 00119 G4ThreeVector theALVelocity(0.,0.,0.); // Prototype of "CMAntiLab" velocity 00120 // Very important! Here the projectile 4M is recalculated in the ZLS (previously in LS) 00121 proj4M = pNucleus.Get4Momentum(); // 4-mom of theProjectileHadron inZLS 00122 G4double pz_per_projectile = proj4M.pz()/pA; // Momentum per nucleon 00123 G4double e_per_projectile = proj4M.e()/pA+targM/tA;// Add MassOfTargetNucleon 00124 G4double vz = pz_per_projectile/e_per_projectile; // Speed of the "oneNuclenCM" 00125 #ifdef debug 00126 G4cout<<"G4QIonIonCollision::Construct: (КщеЯ)Projectile4M="<<proj4M<<", Vz="<<vz 00127 <<", pA="<<pA<<", pE="<<e_per_projectile<<G4endl; 00128 #endif 00129 theCMVelocity.setZ(vz); // CM (w/r to one nucleon) velosity 00130 theProjNucleus.DoLorentzBoost(-theCMVelocity); // BoostProjNucleus to "CM" 00131 theTargNucleus.DoLorentzBoost(-theCMVelocity); // BoostProjNucleus to "CM" 00132 G4LorentzVector prN4M=theProjNucleus.Get4Momentum(); 00133 G4double avz=prN4M.pz()/prN4M.e(); // Speed of AntiLabSys in CMS 00134 theALVelocity.setZ(avz); // CM (w/r to one nucleon) velosity 00135 #ifdef edebug 00136 G4LorentzVector sumP2=theProjNucleus.GetNucleons4Momentum();// Sum of Nucleons4M in RotCM 00137 G4LorentzVector sumT2=theTargNucleus.GetNucleons4Momentum();// Sum of Nucleons4M in RotCM 00138 G4cout<<"-EMC-G4QIonIonCollision::Construct: Boosted, vCM="<<vz<<", vAL="<<avz<<", tN4M=" 00139 <<sumT2<<", pN4M="<<sumP2<<G4endl; 00140 #endif 00141 G4int maxCuts = 7; // @@ Can be reduced for low energies 00142 // 00143 // >>----------->> Find collisions meeting collision conditions and the NN interaction XS 00144 // 00145 G4double outerRadius = theProjNucleus.GetOuterRadius()+theTargNucleus.GetOuterRadius(); 00146 G4QProbability theProbability(2212); // *** proj/targ nucleons *** 00147 // Clean up all previous interactions and reset the counters 00148 #ifdef debug 00149 G4cout<<"G4QIonIonCollision::Construct: theIntSize="<<theInteractions.size()<<G4endl; 00150 #endif 00151 G4int attCnt=-1; 00152 G4int maxAtt=27; 00153 // From here a loop over nucleons of the projectile Nucleus (@@ Z-sorting isn't implem!!) 00154 // 00155 G4QHadron* pNucleon; // Proto of the Projectile Nucleon pointer 00156 G4QHadron* tNucleon; // Proto of the Target Nucleon pointer 00157 G4QNucleus curProjNucleus; 00158 G4QNucleus curTargNucleus; 00159 while(!theInteractions.size() && ++attCnt < maxAtt)// TillAtLeastOneInteractionIsCreated 00160 { 00161 // std::for_each(theInteractions.begin(),theInteractions.end(), DeleteQInteraction()); 00162 // Here we need to clean up the ProjNucleon and the TargNucleon in the Interactions ! 00163 G4int nint=theInteractions.size(); // For the 1st attempt should be zero 00164 for(G4int ni=0; ni < nint; ++ni) 00165 { 00166 G4QInteraction* curInt = theInteractions[ni]; 00167 delete curInt->GetProjectile(); 00168 delete curInt->GetTarget(); 00169 delete curInt; 00170 } 00171 theInteractions.clear(); 00172 // Now we need to make a copy of the projectile and the target nuclei with 3D info 00173 if(attCnt) // This is not theFirstAttempt->Clean 00174 { 00175 curProjNucleus.DeleteNucleons(); 00176 curTargNucleus.DeleteNucleons(); 00177 } 00178 curProjNucleus = theProjNucleus; 00179 curTargNucleus = theTargNucleus; 00180 // choose random impact parameter 00181 std::pair<G4double, G4double> theImpactParameter; 00182 theImpactParameter = curProjNucleus.ChooseImpactXandY(outerRadius); 00183 G4double impactX = theImpactParameter.first; 00184 G4double impactY = theImpactParameter.second; 00185 #ifdef debug 00186 G4cout<<"G4QIonIonCollision::Con:At#="<<attCnt<<",X="<<impactX<<",Y="<<impactY<<G4endl; 00187 #endif 00188 curProjNucleus.StartLoop(); // Prepare Loop ovder nucleons 00189 G4int pnCnt=0; 00190 G4int pnA=curProjNucleus.GetA(); 00191 #ifdef debu 00192 G4cout<<"G4QIonIonCollision::Constr: Before the WHILE over pNucleons,pA="<<pnA<<G4endl; 00193 #endif 00194 while((pNucleon=curProjNucleus.GetNextNucleon()) && pnCnt < pnA) // @@ can be for LOOP? 00195 { 00196 ++pnCnt; 00197 G4QInteractionVector curInteractions; // A temporary vector of interactions 00198 G4LorentzVector pNuc4M=pNucleon->Get4Momentum(); 00199 G4ThreeVector pNucR=pNucleon->GetPosition(); // Position of the pNucleon WRTo pNucl 00200 G4double pImpX = impactX+pNucR.x(); 00201 G4double pImpY = impactY+pNucR.y(); 00202 G4double prEn=proj4M.e(); // For mesons 00203 G4int proB=pNucleus.GetBaryonNumber(); 00204 if (proB>0) prEn-=pNucleus.GetMass(); // For baryons 00205 else if(proB<0) prEn+=mProt; // For anti-baryons 00206 #ifdef debug 00207 G4double impactUsed = 0.; 00208 G4cout<<"G4QIonIonCollision::Construct: estimated energy, prEn="<<prEn<<G4endl; 00209 #endif 00210 G4int totalCuts = 0; 00211 // @@ This is a fake (random) s calculation @@ can be done inside the TARG-while 00212 G4int tnA=curTargNucleus.GetA(); 00213 G4double pImp=std::sqrt(pImpX*pImpX+pImpY*pImpY); 00214 G4double eflen=curTargNucleus.GetThickness(pImp); // EffectiveLength 00215 maxEn=eflen*stringTension; // maxAbsorbedEnergy in IndUnits=MeV 00216 maxNuc=static_cast<int>(eflen*tubeDensity+0.5);// max #0f involved nuclear nucleons 00217 #ifdef debug 00218 G4cout<<"G4QIonIonCollision::Con:pE="<<prEn<<" < mE="<<maxEn<<",mN="<<maxNuc<<G4endl; 00219 #endif 00220 if(prEn < maxEn) // Create DIN interaction and go out 00221 { 00222 G4QHadron* aProjectile = new G4QHadron(*pNucleon);// Copy selected PrNuc for String 00223 curTargNucleus.StartLoop(); // Initialize newSelection forNucleons 00224 tNucleon=curTargNucleus.GetNextNucleon(); // Select a nucleon 00225 G4QHadron* aTarget = new G4QHadron(*tNucleon);// Copy selected nucleon for String 00226 G4QInteraction* anInteraction = new G4QInteraction(aProjectile); 00227 anInteraction->SetTarget(aTarget); 00228 anInteraction->SetNumberOfDINRCollisions(1); // Consider this interaction as DINR 00229 curInteractions.push_back(anInteraction); //--> now the Interaction is not empty 00230 curProjNucleus.DoLorentzBoost(-theALVelocity);// Boost theResPrNucleus toRotAntiLab 00231 curProjNucleus.SubtractNucleon(pNucleon); // Pointer to the used ProjNucleon 00232 curProjNucleus.DoLorentzBoost(theALVelocity);// Boost theResProjNucleus back to CM 00233 curTargNucleus.DoLorentzBoost(theCMVelocity);// Boost theResTgNucleus to RotatedLS 00234 curTargNucleus.SubtractNucleon(tNucleon); // Pointer to the used TargNucleon 00235 curTargNucleus.DoLorentzBoost(-theCMVelocity);// Boost theResTargNucleus back to CM 00236 #ifdef debug 00237 G4cout<<"G4QIonIonCollision::Construct: DINR interaction is created"<<G4endl; 00238 #endif 00239 break; // Break the WHILE of the pNucleons 00240 } 00241 // LOOP over nuclei of the target nucleus to select collisions 00242 curTargNucleus.StartLoop(); // To get the same nucleon 00243 G4int tnCnt = 0; 00244 #ifdef debu 00245 G4cout<<"G4QIonIonCollision::Constr: Before WHILE over tNucleons, tA="<<tnA<<G4endl; 00246 #endif 00247 while((tNucleon=curTargNucleus.GetNextNucleon()) && tnCnt<tnA && totalCuts<maxCuts) 00248 { 00249 ++tnCnt; 00250 G4LorentzVector tNuc4M=tNucleon->Get4Momentum(); 00251 #ifdef debug 00252 G4cout<<"G4QIonIonCollision::Constr: OuterR="<<outerRadius<<", mC="<<maxCuts 00253 <<", pA="<<curProjNucleus<<", tA="<<curTargNucleus<<G4endl; 00254 #endif 00255 // Check the reaction threshold 00256 G4double s_value = (tNuc4M + pNuc4M).mag2(); // Squared CM Energy of compound 00257 G4double ThresholdMass = pNucleon->GetMass() + tNucleon->GetMass(); 00258 #ifdef debug 00259 G4cout<<"G4QInel::Constr: s="<<s_value<<", ThreshM="<<sqr(ThresholdMass)<<G4endl; 00260 #endif 00261 ModelMode = SOFT; // NOT-Diffractive hadronization 00262 if (s_value < 0.) // At ThP=0 is impossible(virtNucl) 00263 { 00264 G4cerr<<"*G4QInelast::Constr:s="<<s_value<<",pN4M="<<pNuc4M<<",tN4M="<<tNuc4M<<G4endl; 00265 G4Exception("G4QIonIonCollision::Construct:","72",FatalException,"LowEn(NegS)"); 00266 } 00267 if (s_value < sqr(ThresholdMass)) // --> Only diffractive interaction 00268 { 00269 #ifdef debug 00270 G4cout<<"G4QIonIonCollision::Constr: ***OnlyDiffraction***, ThM="<<ThresholdMass 00271 <<">sqrt(s)="<<std::sqrt(s_value)<<" -> only Diffraction is possible"<<G4endl; 00272 #endif 00273 ModelMode = DIFFRACTIVE; 00274 } 00275 G4ThreeVector tNucR=tNucleon->GetPosition(); // Position of the tNucleon WRTo tNucl 00276 G4double dImpX = pImpX-tNucR.x(); 00277 G4double dImpY = pImpY-tNucR.y(); 00278 G4double Distance2=dImpX*dImpX+dImpY*dImpY; 00279 #ifdef sdebug 00280 G4cout<<"G4QIonIonCollision::Construct: s="<<s_value<<", D2="<<Distance2<<G4endl; 00281 #endif 00282 // Needs to be moved to Probability class @@ 00283 if(s_value<=10000.) 00284 { 00285 G4cout<<"-Warning-G4QIonIonCollision::Construct: s < .01 GeV^2, p4M=" 00286 <<pNucleon->Get4Momentum()<<",t4M="<<tNucleon->Get4Momentum()<<G4endl; 00287 continue; // skip the rest of the targetNucleons 00288 } 00289 G4double Probability = theProbability.GetPomInelProbability(s_value, Distance2);// P_INEL 00290 // test for inelastic collision 00291 #ifdef sdebug 00292 G4cout<<"G4QIonIonCollision::Construct: Probubility="<<Probability<<G4endl; 00293 #endif 00294 G4double rndNumber = G4UniformRand(); // For the printing purpose 00295 // ModelMode = DIFFRACTIVE; 00296 #ifdef sdebug 00297 G4cout<<"G4QIonIonCollision::Construct: NLOOP prob="<<Probability<<", rndm=" 00298 <<rndNumber<<", d="<<std::sqrt(Distance2)<<G4endl; 00299 #endif 00300 if (Probability > rndNumber) // Inelastic (diffractive or soft) interaction (JOB) 00301 { 00302 G4QHadron* aProjectile = new G4QHadron(*pNucleon);// Copy selected pNuc to String 00303 G4QHadron* aTarget = new G4QHadron(*tNucleon);// Copy selected tNucleon to String 00304 #ifdef edebug 00305 G4cout<<"--->EMC-->G4QIonIonCollision::Construct: TargNucleon is filled, 4M/PDG=" 00306 <<aTarget->Get4Momentum()<<aTarget->GetPDGCode()<<G4endl; 00307 #endif 00308 // Now the energy of the nucleons must be updated in CMS 00309 curProjNucleus.DoLorentzBoost(-theALVelocity);// Boost theResNucleus toRotatedLS 00310 curProjNucleus.SubtractNucleon(pNucleon); // Pointer to the used nucleon 00311 curProjNucleus.DoLorentzBoost(theALVelocity); // Boost theResNucleus back to CM 00312 curTargNucleus.DoLorentzBoost(theCMVelocity); // Boost theResNucleus toRotatedLS 00313 curTargNucleus.SubtractNucleon(tNucleon); // Pointer to the used nucleon 00314 curTargNucleus.DoLorentzBoost(-theCMVelocity);// Boost theResNucleus back to CM 00315 if((theProbability.GetPomDiffProbability(s_value,Distance2)/Probability > 00316 G4UniformRand() && ModelMode==SOFT ) || ModelMode==DIFFRACTIVE) 00317 { 00318 // ------------->> diffractive interaction @@ IsSingleDiffractive called once 00319 if(IsSingleDiffractive()) ExciteSingDiffParticipants(aProjectile, aTarget); 00320 else ExciteDiffParticipants(aProjectile, aTarget); 00321 G4QInteraction* anInteraction = new G4QInteraction(aProjectile); 00322 anInteraction->SetTarget(aTarget); 00323 anInteraction->SetNumberOfDiffractiveCollisions(1); // Why not increment? M.K.? 00324 curInteractions.push_back(anInteraction);//--> now theInteractions not empty 00325 // @@ Why not breake the NLOOP, if only one diffractive can happend? 00326 totalCuts++; // UpdateOfNucleons in't necessary 00327 #ifdef debug 00328 G4cout<<"G4QIonIonCollision::Constr:NLOOP DiffInteract,tC="<<totalCuts<<G4endl; 00329 #endif 00330 } 00331 else 00332 { 00333 // -------------->> nondiffractive = soft interaction 00334 // sample nCut+1 (cut Pomerons) pairs of strings can be produced 00335 G4int nCut; // Result in a chosen number of cuts 00336 G4double* running = new G4double[nCutMax];// @@ This limits the max cuts 00337 for(nCut = 0; nCut < nCutMax; nCut++) // Calculates multiCut probabilities 00338 { 00339 running[nCut]= theProbability.GetCutPomProbability(s_value, Distance2, nCut+1); 00340 if(nCut) running[nCut] += running[nCut-1];// Sum up with the previous one 00341 } 00342 G4double random = running[nCutMax-1]*G4UniformRand(); 00343 for(nCut = 0; nCut < nCutMax; nCut++) if(running[nCut] > random) break; // tNuc 00344 delete [] running; 00345 #ifdef debug 00346 G4cout<<"G4QIonIonCollision::Construct: NLOOP-Soft Chosen nCut="<<nCut<<G4endl; 00347 #endif 00348 // @@ If nCut>0 interaction with a few nucleons is possible 00349 // @@ nCut is found with big efforts and now nCut=0 ?? M.K. ?? !! 00350 //nCut = 0; // @@ in original code ?? @@ 00351 aTarget->IncrementCollisionCount(nCut+1); // @@ What about multyNucleon target? 00352 aProjectile->IncrementCollisionCount(nCut+1); 00353 G4QInteraction* anInteraction = new G4QInteraction(aProjectile); 00354 anInteraction->SetTarget(aTarget); 00355 anInteraction->SetNumberOfSoftCollisions(nCut+1); 00356 curInteractions.push_back(anInteraction); // Now curInteractions are not empty 00357 totalCuts += nCut+1; 00358 #ifdef debug 00359 G4cout<<"G4QIonIonCollision::Constr:NLOOP SoftInteract,tC="<<totalCuts<<G4endl; 00360 impactUsed=Distance2; 00361 #endif 00362 } 00363 }// Probability selection 00364 } // End of While over target nucleons 00365 G4int nCurInt=curInteractions.size(); 00366 for(G4int ni=0; ni < nCurInt; ++ni) theInteractions.push_back(curInteractions[ni]); 00367 curInteractions.clear(); // Probably, not necessary... 00368 } // End of WHILE over projectile nucleons 00369 } // End of WHILE over attempts to find at least one nucleus-nucleus interaction 00370 theProjNucleus.DeleteNucleons(); 00371 theTargNucleus.DeleteNucleons(); 00372 theProjNucleus = curProjNucleus; 00373 theTargNucleus = curTargNucleus; 00374 curProjNucleus.DeleteNucleons(); 00375 curTargNucleus.DeleteNucleons(); 00376 G4int nInt=theInteractions.size(); 00377 #ifdef debug 00378 G4cout<<"G4QIonIonCollision::Constr: #ofInteractions = "<<nInt<<", #OfDINR = " 00379 <<theInteractions[0]->GetNumberOfDINRCollisions()<<G4endl; 00380 #endif 00381 if(!nInt || (nInt==1 && theInteractions[0]->GetNumberOfDINRCollisions()==1)) // QFreeInel 00382 { 00383 G4QHadron* aTarget; 00384 G4QHadron* aProjectile; 00385 if(nInt) // Take Targ/Proj from the Interaction 00386 { 00387 aTarget=theInteractions[0]->GetTarget(); 00388 aProjectile=theInteractions[0]->GetProjectile(); 00389 delete theInteractions[0]; 00390 theInteractions.clear(); 00391 } 00392 else // Create a new target nucleon (?) 00393 { 00394 theProjNucleus.StartLoop(); // To get the same nucleon from projec 00395 pNucleon=theProjNucleus.GetNextNucleon(); // Get theNucleon to create projectNuc 00396 aProjectile = new G4QHadron(*pNucleon); // Copy selected pNucleon for String 00397 theProjNucleus.DoLorentzBoost(-theALVelocity); // Boost theResProjNucleus toRotatedLS 00398 theProjNucleus.SubtractNucleon(pNucleon); // Pointer to SelProjNucleon to delete 00399 theProjNucleus.DoLorentzBoost(theALVelocity); // Boost theResProNucleus back to CMS 00400 theTargNucleus.StartLoop(); // To get the same nucleon from target 00401 tNucleon=theTargNucleus.GetNextNucleon(); // Get theNucleon to create targetNucl 00402 aTarget = new G4QHadron(*tNucleon); // Copy selected tNucleon for String 00403 theTargNucleus.DoLorentzBoost(theCMVelocity); // Boost theResTargNucleus toRotatedLS 00404 theTargNucleus.SubtractNucleon(tNucleon); // Pointer to SelTargNucleon to delete 00405 theTargNucleus.DoLorentzBoost(-theCMVelocity); // Boost theResTargNucleus back to CMS 00406 } 00407 G4QContent QQC=aTarget->GetQC()+aProjectile->GetQC(); // QContent of the compound 00408 G4LorentzVector Q4M=aTarget->Get4Momentum()+aProjectile->Get4Momentum(); // 4-mom of Q 00409 delete aTarget; 00410 delete aProjectile; 00411 // @@ 4-Mom should be converted to LS for the TargQuasmon and to AL for the ProjQuasmon 00412 Q4M.boost(theCMVelocity); 00413 Q4M=toLab*Q4M; 00414 G4Quasmon* stringQuasmon = new G4Quasmon(QQC, Q4M); 00415 theQuasmons.push_back(stringQuasmon); 00416 theTargNucleus.DoLorentzBoost(theCMVelocity); // BoostTheResidualNucleus toRotatedLS 00417 theTargNucleus.DoLorentzRotation(toLab); // Recove Z-direction in LS ("LS"->LS) 00418 return; 00419 } 00420 // 00421 // ------------------ now build the parton pairs for the strings ------------------ 00422 // 00423 for(G4int i=0; i<nInt; i++) 00424 { 00425 theInteractions[i]->SplitHadrons(); 00426 #ifdef edebug 00427 G4QHadron* projH=theInteractions[i]->GetProjectile(); // Projectile of theInteraction 00428 G4QHadron* targH=theInteractions[i]->GetTarget(); // Target of the Interaction 00429 G4LorentzVector pSP(0.,0.,0.,0.); // Sum of parton's 4mom's for proj 00430 G4LorentzVector tSP(0.,0.,0.,0.); // Sum of parton's 4mom's for proj 00431 std::list<G4QParton*> projCP=projH->GetColor(); // Pointers to proj Color-partons 00432 std::list<G4QParton*> projAC=projH->GetAntiColor();// PointersTo projAntiColorPartons 00433 std::list<G4QParton*> targCP=targH->GetColor(); // Pointers to targ Color-partons 00434 std::list<G4QParton*> targAC=targH->GetAntiColor();// PointersTo targAntiColorPartons 00435 std::list<G4QParton*>::iterator picp = projCP.begin(); 00436 std::list<G4QParton*>::iterator pecp = projCP.end(); 00437 std::list<G4QParton*>::iterator piac = projAC.begin(); 00438 std::list<G4QParton*>::iterator peac = projAC.end(); 00439 std::list<G4QParton*>::iterator ticp = targCP.begin(); 00440 std::list<G4QParton*>::iterator tecp = targCP.end(); 00441 std::list<G4QParton*>::iterator tiac = targAC.begin(); 00442 std::list<G4QParton*>::iterator teac = targAC.end(); 00443 for(; picp!=pecp&& piac!=peac&& ticp!=tecp&& tiac!=teac; ++picp,++piac,++ticp,++tiac) 00444 { 00445 pSP+=(*picp)->Get4Momentum(); 00446 pSP+=(*piac)->Get4Momentum(); 00447 tSP+=(*ticp)->Get4Momentum(); 00448 tSP+=(*tiac)->Get4Momentum(); 00449 } 00450 G4cout<<"-EMC-G4QIonIonCollision::Construct: Interaction#"<<i<<",dP4M=" 00451 <<projH->Get4Momentum()-pSP<<",dT4M="<<targH->Get4Momentum()-tSP<<G4endl; 00452 #endif 00453 } 00454 // 00455 // >>........>> make soft collisions (ordering is vital) 00456 // 00457 G4QInteractionVector::iterator it; 00458 #ifdef debug 00459 G4cout<<"G4QIonIonCollision::Constr: Creation ofSoftCollisionPartonPair STARTS"<<G4endl; 00460 #endif 00461 G4bool rep=true; 00462 while(rep && theInteractions.size()) 00463 { 00464 for(it = theInteractions.begin(); it != theInteractions.end(); ++it) 00465 { 00466 G4QInteraction* anIniteraction = *it; 00467 G4QPartonPair* aPair=0; 00468 G4int nSoftCollisions = anIniteraction->GetNumberOfSoftCollisions(); 00469 #ifdef debug 00470 G4cout<<"G4QIonIonCollision::Construct: #0f SoftCollisions ="<<nSoftCollisions<<G4endl; 00471 #endif 00472 if (nSoftCollisions) 00473 { 00474 G4QHadron* pProjectile = anIniteraction->GetProjectile(); 00475 G4QHadron* pTarget = anIniteraction->GetTarget(); 00476 for (G4int j = 0; j < nSoftCollisions; j++) 00477 { 00478 aPair = new G4QPartonPair(pTarget->GetNextParton(), 00479 pProjectile->GetNextAntiParton(), 00480 G4QPartonPair::SOFT, G4QPartonPair::TARGET); 00481 thePartonPairs.push_back(aPair); // A target pair (Why TAGRET?) 00482 aPair = new G4QPartonPair(pProjectile->GetNextParton(), 00483 pTarget->GetNextAntiParton(), 00484 G4QPartonPair::SOFT, G4QPartonPair::PROJECTILE); 00485 thePartonPairs.push_back(aPair); // A projectile pair (Why Projectile?) 00486 #ifdef debug 00487 G4cout<<"--->G4QIonIonCollision::Constr: SOFT, 2 parton pairs are filled"<<G4endl; 00488 #endif 00489 } 00490 delete *it; 00491 it=theInteractions.erase(it); // SoftInteractions're converted&erased 00492 if( it != theInteractions.begin() )// To avoid going below begin() (for Windows) 00493 { 00494 it--; 00495 rep=false; 00496 } 00497 else 00498 { 00499 rep=true; 00500 break; 00501 } 00502 } 00503 else rep=false; 00504 } 00505 } 00506 #ifdef debug 00507 G4cout<<"G4QIonIonCollision::Constr: -> Parton pairs for SOFT strings are made"<<G4endl; 00508 #endif 00509 // 00510 // >>.............>> make the rest as the diffractive interactions 00511 // 00512 for(unsigned i = 0; i < theInteractions.size(); i++) // Interactions are reduced bySoft 00513 { 00514 // The double or single diffraction is defined by the presonce of proj/targ partons 00515 G4QInteraction* anIniteraction = theInteractions[i]; 00516 G4QPartonPair* aPartonPair; 00517 #ifdef debug 00518 G4cout<<"G4QIonIonCollision::Constr: CreationOfDiffractivePartonPairs, i="<<i<<G4endl; 00519 #endif 00520 // the projectile diffraction parton pair is created first 00521 G4QHadron* aProjectile = anIniteraction->GetProjectile(); 00522 G4QParton* aParton = aProjectile->GetNextParton(); 00523 if (aParton) 00524 { 00525 aPartonPair = new G4QPartonPair(aParton, aProjectile->GetNextAntiParton(), 00526 G4QPartonPair::DIFFRACTIVE, 00527 G4QPartonPair::PROJECTILE); 00528 thePartonPairs.push_back(aPartonPair); 00529 #ifdef debug 00530 G4cout<<"G4QIonIonCollision::Constr: proj Diffractive PartonPair is filled"<<G4endl; 00531 #endif 00532 } 00533 // then the target diffraction parton pair is created 00534 G4QHadron* aTarget = anIniteraction->GetTarget(); 00535 aParton = aTarget->GetNextParton(); 00536 if (aParton) 00537 { 00538 aPartonPair = new G4QPartonPair(aParton, aTarget->GetNextAntiParton(), 00539 G4QPartonPair::DIFFRACTIVE, G4QPartonPair::TARGET); 00540 thePartonPairs.push_back(aPartonPair); 00541 #ifdef debug 00542 G4cout<<"G4QIonIonCollision::Constr: target Diffractive PartonPair's filled"<<G4endl; 00543 #endif 00544 } 00545 } 00546 #ifdef debug 00547 G4cout<<"G4QIonIonCollision::Construct: DiffractivePartonPairs are created"<<G4endl; 00548 #endif 00549 // 00550 // >>.......>> clean-up Interactions, if necessary 00551 // 00552 std::for_each(theInteractions.begin(),theInteractions.end(), DeleteQInteraction()); 00553 theInteractions.clear(); 00554 #ifdef debug 00555 G4cout<<"G4QIonIonCollision::Construct: Temporary objects are cleaned up"<<G4endl; 00556 #endif 00557 // This function prepares theBoost for transformation of secondaries to LS (-ProjRot!) 00558 theProjNucleus.DoLorentzBoost(theCMVelocity);// Boost theResidualProjNucleus to RotatedLS 00559 theTargNucleus.DoLorentzBoost(theCMVelocity);// Boost theResidualTargNucleus to RotatedLS 00560 // @@ Nucleus isn't completely in LS, it needs the toZ (-ProjRot) rotation to consE/M 00561 #ifdef debug 00562 G4cout<<">>........>>G4QIonIonCollision::Construct: >>..>> Strings are created "<<G4endl; 00563 #endif 00564 G4QPartonPair* aPair; 00565 G4QString* aString=0; 00566 while(thePartonPairs.size()) // @@ At present noDifference in stringBuild (? M.K.) 00567 { 00568 aPair = thePartonPairs.back(); // Get the parton pair 00569 thePartonPairs.pop_back(); // Clean up thePartonPairPointer in the Vector 00570 #ifdef debug 00571 G4cout<<"G4QIonIonCollision::Construct:StringType="<<aPair->GetCollisionType()<<G4endl; 00572 #endif 00573 aString= new G4QString(aPair); 00574 #ifdef debug 00575 G4cout<<"G4QIonIonCollision::Construct:NewString4M="<<aString->Get4Momentum()<<G4endl; 00576 #endif 00577 aString->Boost(theCMVelocity); // ! Strings are moved to ZLS when pushed ! 00578 strings.push_back(aString); 00579 stringsInitted=true; 00580 delete aPair; 00581 } // End of the String Creation LOOP 00582 #ifdef edebug 00583 G4LorentzVector sum=theProjNucleus.Get4Momentum()+theTargNucleus.Get4Momentum(); // in LS 00584 G4int rChg=totChg-theProjNucleus.GetZ()-theTargNucleus.GetZ(); 00585 G4int rBaN=totBaN-theProjNucleus.GetA()-theTargNucleus.GetA(); 00586 G4int nStrings=strings.size(); 00587 G4cout<<"-EMC-G4QIonIonCollision::Constr:#ofString="<<nStrings<<",tNuc4M="<<sum<<G4endl; 00588 for(G4int i=0; i<nStrings; i++) 00589 { 00590 G4QString* prString=strings[i]; 00591 G4LorentzVector strI4M=prString->Get4Momentum(); 00592 sum+=strI4M; 00593 G4int sChg=prString->GetCharge(); 00594 G4int sBaN=prString->GetBaryonNumber(); 00595 G4int LPDG=prString->GetLeftParton()->GetPDGCode(); 00596 G4int RPDG=prString->GetRightParton()->GetPDGCode(); 00597 G4QContent LQC =prString->GetLeftParton()->GetQC(); 00598 G4QContent RQC =prString->GetRightParton()->GetQC(); 00599 rChg-=sChg; 00600 rBaN-=sBaN; 00601 G4cout<<"-EMC-G4QIonIonCollision::Construct: String#"<<i<<", 4M="<<strI4M<<", LPDG=" 00602 <<LPDG<<LQC<<",RPDG="<<RPDG<<RQC<<", Ch="<<sChg<<", BN="<<sBaN<<G4endl; 00603 } 00604 G4cout<<"-EMC-G4QInel::Constr: r4M="<<sum-totZLS4M<<", rC="<<rChg<<", rB="<<rBaN<<G4endl; 00605 #endif 00606 if(!stringsInitted) 00607 { 00608 G4cerr<<"*****G4QIonIonCollision::Construct:**** No strings are created *****"<<G4endl; 00609 G4Exception("G4QIonIonCollision::Construct:","72",FatalException,"No Strings created"); 00610 } 00611 #ifdef debug 00612 G4cout<<"G4QIonIonCollision::Constr:BeforeRotation, #0fStrings="<<strings.size()<<G4endl; 00613 #endif 00614 // 00615 // ---------------- At this point the strings must be already created in "LS" ----------- 00616 // 00617 for(unsigned astring=0; astring < strings.size(); astring++) 00618 strings[astring]->LorentzRotate(toLab); // Recove Z-direction in LS ("LS"->LS) 00619 theProjNucleus.DoLorentzRotation(toLab); // Recove Z-dir in LS ("LS"->LS) for ProjNucleus 00620 theTargNucleus.DoLorentzRotation(toLab); // Recove Z-dir in LS ("LS"->LS) for TargNucleus 00621 // Now everything is in LS system 00622 #ifdef edebug 00623 G4LorentzVector sm=theProjNucleus.Get4Momentum()+theTargNucleus.Get4Momentum();// NucInLS 00624 G4int rCg=totChg-theProjNucleus.GetZ()-theTargNucleus.GetZ(); 00625 G4int rBC=totBaN-theProjNucleus.GetA()-theTargNucleus.GetA(); 00626 G4int nStrs=strings.size(); 00627 G4cout<<"-EMCLS-G4QIonIonCollision::Constr:#ofS="<<nStrings<<",Nuc4M(E=M)="<<sum<<G4endl; 00628 for(G4int i=0; i<nStrs; i++) 00629 { 00630 G4LorentzVector strI4M=strings[i]->Get4Momentum(); 00631 sm+=strI4M; 00632 G4int sChg=strings[i]->GetCharge(); 00633 rCg-=sChg; 00634 G4int sBaN=strings[i]->GetBaryonNumber(); 00635 rBC-=sBaN; 00636 G4cout<<"-EMCLS-G4QInel::Construct:String#"<<i<<",4M="<<strI4M<<strI4M.m()<<",Charge=" 00637 <<sChg<<",BaryN="<<sBaN<<G4endl; 00638 } 00639 G4cout<<"-EMCLS-...G4QInel::Constr: r4M="<<sm-totLS4M<<", rC="<<rCg<<",rB="<<rBC<<G4endl; 00640 #endif 00641 // 00642 // --- Strings are created, but we should try to get rid of negative mass strings ----- 00643 // 00644 SwapPartons(); 00645 // 00646 // --- Strings are created, but we should get rid of too light strings (Mmin+MPi0) ----- 00647 // 00648 G4int problem=0; // 0="no problem", incremented by ASIS 00649 G4QStringVector::iterator ist; 00650 G4bool con=true; 00651 while(con && strings.size()) 00652 { 00653 for(ist = strings.begin(); ist < strings.end(); ++ist) 00654 { 00655 G4bool bad=true; 00656 G4LorentzVector cS4M=(*ist)->Get4Momentum(); 00657 G4double cSM2=cS4M.m2(); // Squared mass of the String 00658 G4QParton* cLeft=(*ist)->GetLeftParton(); 00659 G4QParton* cRight=(*ist)->GetRightParton(); 00660 G4int cLT=cLeft->GetType(); 00661 G4int cRT=cRight->GetType(); 00662 G4int cLPDG=cLeft->GetPDGCode(); 00663 G4int cRPDG=cRight->GetPDGCode(); 00664 G4int aLPDG=0; 00665 G4int aRPDG=0; 00666 if (cLPDG > 7) aLPDG= cLPDG/100; 00667 else if(cLPDG <-7) aLPDG=(-cLPDG)/100; 00668 if (cRPDG > 7) aRPDG= cRPDG/100; 00669 else if(cRPDG <-7) aRPDG=(-cRPDG)/100; 00670 G4int L1=0; 00671 G4int L2=0; 00672 if(aLPDG) 00673 { 00674 L1=aLPDG/10; 00675 L2=aLPDG%10; 00676 } 00677 G4int R1=0; 00678 G4int R2=0; 00679 if(aRPDG) 00680 { 00681 R1=aRPDG/10; 00682 R2=aRPDG%10; 00683 } 00684 G4double cSM=cSM2; 00685 if(cSM2>0.) cSM=std::sqrt(cSM2); 00686 #ifdef debug 00687 G4cout<<"G4QIonIonCollision::Construct: Needs Fusion? cLPDG="<<cLPDG<<",cRPDG="<<cRPDG 00688 <<",cM(cM2 if neg)="<<cSM<<G4endl; 00689 #endif 00690 if(cSM>0.) // Mass can be calculated 00691 { 00692 G4bool single=true; 00693 G4double miM=0.; // Proto of the Min HadronString Mass 00694 if(cLT==2 && cRT==2) 00695 { 00696 if(L1!=R1 && L1!=R2 && L2!=R1 && L2!=R2) // Unreducable DiQ-aDiQ 00697 { 00698 single=false; 00699 G4QPDGCode tmp; 00700 std::pair<G4int,G4int> paB=tmp.MakeTwoBaryons(L1, L2, R1, R2); 00701 miM=G4QPDGCode(paB.first).GetMass()+G4QPDGCode(paB.second).GetMass(); 00702 } 00703 } 00704 if(single) miM=G4QPDGCode((*ist)->GetQC().GetSPDGCode()).GetMass() + mPi0;//MinHSMass 00705 //if(single) miM=G4QPDGCode((*ist)->GetQC().GetSPDGCode()).GetMass();//MinHSMass 00706 #ifdef debug 00707 G4cout<<"G4QInel::Const:*IsItGood? realM="<<std::sqrt(cSM2)<<" > GSM="<<miM<<G4endl; 00708 #endif 00709 if(std::sqrt(cSM2) > miM) bad=false; // String is OK 00710 } 00711 if(bad) // String should be merged with others 00712 { 00713 #ifdef debug 00714 G4cout<<"G4QInel::Const:TryFuse,L1="<<L1<<",L2="<<L2<<",R1="<<R1<<",R2="<<R2<<G4endl; 00715 #endif 00716 G4int cST=cLT+cRT; 00717 G4double excess=-DBL_MAX; // The value to be maximized excess M 00718 G4double maxiM2=-DBL_MAX; // The value to be maximized M2 00719 G4QStringVector::iterator sst; // Selected partner string 00720 G4QStringVector::iterator pst; 00721 G4int sLPDG=0; // selectedLeft (like inStringPartner) 00722 G4int sRPDG=0; // selectedRight(like inStringPartner) 00723 G4int sOrd=0; // selected Order of PartonFusion 00724 G4bool minC=true; // for the case when M2<0 00725 if(cSM2>0.) minC=false; // If M2>0 already don'tSearchFor M2>0 00726 for(pst = strings.begin(); pst < strings.end(); pst++) if(pst != ist) 00727 { 00728 G4LorentzVector pS4M=(*pst)->Get4Momentum()+cS4M; // Summed 4-momentum 00729 G4int nLPDG=0; // new Left (like in theStringPartner) 00730 G4int nRPDG=0; // new Right(like in theStringPartner) 00731 G4double pExcess=-DBL_MAX; // Prototype of the excess 00732 G4double pSM2=pS4M.m2(); // Squared mass of the Fused Strings 00733 #ifdef debug 00734 G4cout<<"->G4QIonIonCollision::Construct: sum4M="<<pS4M<<", M2="<<pSM2<<G4endl; 00735 #endif 00736 //if(pSM2>0.) // The partner can be a candidate 00737 //{ 00738 G4QParton* pLeft=(*pst)->GetLeftParton(); 00739 G4QParton* pRight=(*pst)->GetRightParton(); 00740 G4int pLT=pLeft->GetType(); 00741 G4int pRT=pRight->GetType(); 00742 G4int pLPDG=pLeft->GetPDGCode(); 00743 G4int pRPDG=pRight->GetPDGCode(); 00744 G4int LPDG=0; 00745 G4int RPDG=0; 00746 if (pLPDG > 7) LPDG= pLPDG/100; 00747 else if(pLPDG <-7) LPDG=(-pLPDG)/100; 00748 if (pRPDG > 7) RPDG= pRPDG/100; 00749 else if(pRPDG <-7) RPDG=(-pRPDG)/100; 00750 G4int pL1=0; 00751 G4int pL2=0; 00752 if(LPDG) 00753 { 00754 pL1=LPDG/10; 00755 pL2=LPDG%10; 00756 } 00757 G4int pR1=0; 00758 G4int pR2=0; 00759 if(RPDG) 00760 { 00761 pR1=RPDG/10; 00762 pR2=RPDG%10; 00763 } 00764 G4int pST=pLT+pRT; 00765 #ifdef debug 00766 G4cout<<"G4QInel::Construct: Partner/w pLPDG="<<pLPDG<<", pRPDG="<<pRPDG<<", pM2=" 00767 <<pSM2<<G4endl; 00768 #endif 00769 // Different fromCompactAlrorithm ofStringFusionAfterDecay (no DiQaDiQ reduction) 00770 G4int tf=0; // Type combination flag 00771 G4int af=0; // Annihilatio combination flag 00772 if (cST==2 && pST==2) // QaQ + QaQ = DiQaDiQ (always) 00773 { 00774 tf=1; 00775 af=1; 00776 } 00777 else if(cST==2 && pST==3) // QaQ + QDiQ/aQaDiQ = QDiQ/aQaDiQ (s) 00778 { 00779 tf=2; 00780 if (pLPDG > 7 && 00781 ( (cLPDG<0 && (-cLPDG==pL1 || -cLPDG==pL2 || -cLPDG==pRPDG) ) || 00782 (cRPDG<0 && (-cRPDG==pL1 || -cRPDG==pL2 || -cRPDG==pRPDG) ) 00783 ) 00784 ) af=1; 00785 else if(pRPDG > 7 && 00786 ( (cLPDG<0 && (-cLPDG==pR1 || -cLPDG==pR2 || -cLPDG==pLPDG) ) || 00787 (cRPDG<0 && (-cRPDG==pR1 || -cRPDG==pR2 || -cRPDG==pLPDG) ) 00788 ) 00789 ) af=2; 00790 else if(pLPDG <-7 && 00791 ( (cLPDG>0 && ( cLPDG==pL1 || cLPDG==pL2 || cLPDG==-pRPDG) ) || 00792 (cRPDG>0 && ( cRPDG==pL1 || cRPDG==pL2 || cRPDG==-pRPDG) ) 00793 ) 00794 ) af=3; 00795 else if(pRPDG <-7 && 00796 ( (cLPDG>0 && ( cLPDG==pR1 || cLPDG==pR2 || cLPDG==-pLPDG) ) || 00797 (cRPDG>0 && ( cRPDG==pR1 || cRPDG==pR2 || cRPDG==-pLPDG) ) 00798 ) 00799 ) af=4; 00800 #ifdef debug 00801 else G4cout<<"G4QIonIonCollision::Constr:2(QaQ+QDiQ/aQaDiQ) Can't fuse"<<G4endl; 00802 #endif 00803 } 00804 else if(cST==3 && pST==2) // QDiQ/aQaDiQ + QaQ = QDiQ/aQaDiQ (s) 00805 { 00806 tf=3; 00807 if (cLPDG > 7 && 00808 ( (pLPDG<0 && (-pLPDG==L1 || -pLPDG==L2) ) || 00809 (pRPDG<0 && (-pRPDG==L1 || -pRPDG==L2) ) 00810 ) 00811 ) af=1; 00812 else if(cRPDG > 7 && 00813 ( (pLPDG<0 && (-pLPDG==R1 || -pLPDG==R2) ) || 00814 (pRPDG<0 && (-pRPDG==R1 || -pRPDG==R2) ) 00815 ) 00816 ) af=2; 00817 else if(cLPDG <-7 && 00818 ( (pLPDG>0 && ( pLPDG==L1 || pLPDG==L2) ) || 00819 (pRPDG>0 && ( pRPDG==L1 || pRPDG==L2) ) 00820 ) 00821 ) af=3; 00822 else if(cRPDG <-7 && 00823 ( (pLPDG>0 && ( pLPDG==R1 || pLPDG==R2) ) || 00824 (pRPDG>0 && ( pRPDG==R1 || pRPDG==R2) ) 00825 ) 00826 ) af=4; 00827 #ifdef debug 00828 else G4cout<<"G4QIonIonCollision::Constr:3(QDiQ/aQaDiQ+QaQ) Can't fuse"<<G4endl; 00829 #endif 00830 } 00831 else if(cST==2 && pST==4) // QaQ + aDiQDiQ = QaQ (double) 00832 { 00833 tf=4; 00834 if (pLPDG > 7) // pRPDG <-7 00835 { 00836 if ( (-cLPDG==pL1 || -cLPDG==pL2) && (cRPDG==pR1 || cRPDG==pR2) ) af=1; 00837 else if( (-cRPDG==pL1 || -cRPDG==pL2) && (cLPDG==pR1 || cLPDG==pR2) ) af=2; 00838 } 00839 else if(pRPDG > 7) // pLPDG <-7 00840 { 00841 if ( (-cRPDG==pR1 || -cRPDG==pR2) && (cLPDG==pL1 || cLPDG==pL2) ) af=3; 00842 else if( (-cLPDG==pR1 || -cLPDG==pR2) && (cRPDG==pL1 || cRPDG==pL2) ) af=4; 00843 } 00844 #ifdef debug 00845 else G4cout<<"-G4QIonIonCollision::Constr: 4 (QaQ+aQDiQDiQ) Can't fuse"<<G4endl; 00846 #endif 00847 } 00848 else if(cST==4 && pST==2) // aDiQDiQ + QaQ = QaQ (double) 00849 { 00850 tf=5; 00851 if (cLPDG > 7) // cRPDG<-7 00852 { 00853 if ( (-pLPDG==L1 || -pLPDG==L2) && (pRPDG==R1 || pRPDG==R2) ) af=1; 00854 else if( (-pRPDG==L1 || -pRPDG==L2) && (pLPDG==R1 || pLPDG==R2) ) af=2; 00855 } 00856 else if(cRPDG > 7) // cLPDG<-7 00857 { 00858 if ( (-pRPDG==R1 || -pRPDG==R2) && (pLPDG==L1 || pLPDG==L2) ) af=3; 00859 else if( (-pLPDG==R1 || -pLPDG==R2) && (pRPDG==L1 || pRPDG==L2) ) af=4; 00860 } 00861 #ifdef debug 00862 else G4cout<<"-G4QIonIonCollision::Constr: 5 (aQDiQDiQ+QaQ) Can't fuse"<<G4endl; 00863 #endif 00864 } 00865 else if(cST==3 && pST==3) // QDiQ + aQaDiQ = QaQ (double) 00866 { 00867 tf=6; 00868 if(pLPDG > 7) 00869 { 00870 if (cLPDG<-7 && (-cRPDG==pL1 || -cRPDG==pL2) && (pRPDG==L1 || pRPDG==L2)) 00871 af=1; 00872 else if(cRPDG<-7 && (-cLPDG==pL1 || -cLPDG==pL2) && (pRPDG==R1 || pRPDG==R2)) 00873 af=2; 00874 } 00875 else if(pRPDG > 7) 00876 { 00877 if (cLPDG<-7 && (-cRPDG==pR1 || -cRPDG==pR2) && (pLPDG==L1 || pLPDG==L2)) 00878 af=3; 00879 else if(cRPDG<-7 && (-cLPDG==pR1 || -cLPDG==pR2) && (pLPDG==R1 || pLPDG==R2)) 00880 af=4; 00881 } 00882 else if(cLPDG > 7) 00883 { 00884 if (pLPDG<-7 && (-pRPDG==L1 || -pRPDG==L2) && (cRPDG==pL1 || cRPDG==pL2)) 00885 af=5; 00886 else if(pRPDG<-7 && (-pLPDG==L1 || -pLPDG==L2) && (cRPDG==pR1 || cRPDG==pR2)) 00887 af=6; 00888 } 00889 else if(cRPDG > 7) 00890 { 00891 if (pLPDG<-7 && (-pRPDG==R1 || -pRPDG==R2) && (cLPDG==pL1 || cLPDG==pL2)) 00892 af=7; 00893 else if(pRPDG<-7 && (-pLPDG==R1 || -pLPDG==R2) && (cLPDG==pR1 || cLPDG==pR2)) 00894 af=8; 00895 } 00896 #ifdef debug 00897 else G4cout<<"-G4QIonIonCollision::Constr: 6 (QDiQ+aQaDiQ) Can't fuse"<<G4endl; 00898 #endif 00899 } 00900 #ifdef debug 00901 G4cout<<"G4QIonIonCollision::Constr: **Possibility**, tf="<<tf<<",af="<<af<<G4endl; 00902 #endif 00903 if(tf && af) 00904 { 00905 // Strings can be fused, update the max excess and remember usefull switches 00906 G4int order=0; // LL/RR (1) or LR/RL (2) PartonFusion 00907 switch (tf) 00908 { 00909 case 1: // ------------------------------------> QaQ + QaQ = DiQaDiQ (always) 00910 if (cLPDG > 0 && pLPDG > 0) 00911 { 00912 order= 1; 00913 if (cLPDG > pLPDG) nLPDG=cLPDG*1000+pLPDG*100+1; 00914 else if(cLPDG < pLPDG) nLPDG=pLPDG*1000+cLPDG*100+1; 00915 else nLPDG=pLPDG*1000+cLPDG*100+3; 00916 if (cRPDG < pRPDG) nRPDG=cRPDG*1000+pRPDG*100-1; 00917 else if(cRPDG > pRPDG) nRPDG=pRPDG*1000+cRPDG*100-1; 00918 else nRPDG=pRPDG*1000+cRPDG*100-3; 00919 } 00920 else if(cLPDG < 0 && pLPDG < 0) 00921 { 00922 order= 1; 00923 if (cRPDG > pRPDG) nRPDG=cRPDG*1000+pRPDG*100+1; 00924 else if(cRPDG < pRPDG) nRPDG=pRPDG*1000+cRPDG*100+1; 00925 else nRPDG=pRPDG*1000+cRPDG*100+3; 00926 if (cLPDG < pLPDG) nLPDG=cLPDG*1000+pLPDG*100-1; 00927 else if(cLPDG > pLPDG) nLPDG=pLPDG*1000+cLPDG*100-1; 00928 else nLPDG=pLPDG*1000+cLPDG*100-3; 00929 } 00930 else if(cRPDG > 0 && pLPDG > 0) 00931 { 00932 order=-1; 00933 if (cRPDG > pLPDG) nLPDG=cRPDG*1000+pLPDG*100+1; 00934 else if(cRPDG < pLPDG) nLPDG=pLPDG*1000+cRPDG*100+1; 00935 else nLPDG=pLPDG*1000+cRPDG*100+3; 00936 if (cLPDG < pRPDG) nRPDG=cLPDG*1000+pRPDG*100-1; 00937 else if(cLPDG > pRPDG) nRPDG=pRPDG*1000+cLPDG*100-1; 00938 else nRPDG=pRPDG*1000+cLPDG*100-3; 00939 } 00940 else if(cRPDG < 0 && pLPDG < 0) 00941 { 00942 order=-1; 00943 if (cLPDG > pRPDG) nRPDG=cLPDG*1000+pRPDG*100+1; 00944 else if(cLPDG < pRPDG) nRPDG=pRPDG*1000+cLPDG*100+1; 00945 else nRPDG=pRPDG*1000+cLPDG*100+3; 00946 if (cRPDG < pLPDG) nLPDG=cRPDG*1000+pLPDG*100-1; 00947 else if(cRPDG > pLPDG) nLPDG=pLPDG*1000+cRPDG*100-1; 00948 else nLPDG=pLPDG*1000+cRPDG*100-3; 00949 } 00950 break; 00951 case 2: // ------------------------> QaQ + QDiQ/aQaDiQ = QDiQ/aQaDiQ (single) 00952 switch (af) 00953 { 00954 case 1: // ....................... pLPDG > 7 00955 if(cLPDG < 0) 00956 { 00957 order= 1; 00958 if(-cLPDG==pRPDG) 00959 { 00960 nLPDG=pLPDG; 00961 nRPDG=cRPDG; 00962 } 00963 else 00964 { 00965 if (cRPDG > pRPDG) nRPDG=cRPDG*1000+pRPDG*100+1; 00966 else if(cRPDG < pRPDG) nRPDG=pRPDG*1000+cRPDG*100+1; 00967 else nRPDG=pRPDG*1000+cRPDG*100+3; 00968 if (-cLPDG == pL1) nLPDG=pL2; 00969 else nLPDG=pL1; // -cLPDG == pL2 00970 } 00971 } 00972 else // cRPDG < 0 00973 { 00974 order=-1; 00975 if(-cRPDG==pRPDG) 00976 { 00977 nLPDG=pLPDG; 00978 nRPDG=cLPDG; 00979 } 00980 else 00981 { 00982 if (cLPDG > pRPDG) nRPDG=cLPDG*1000+pRPDG*100+1; 00983 else if(cLPDG < pRPDG) nRPDG=pRPDG*1000+cLPDG*100+1; 00984 else nRPDG=pRPDG*1000+cLPDG*100+3; 00985 if (-cRPDG == pL1) nLPDG=pL2; 00986 else nLPDG=pL1; // -cRPDG == pL2 00987 } 00988 } 00989 break; 00990 case 2: // ....................... pRPDG > 7 00991 if(cLPDG < 0) 00992 { 00993 order=-1; 00994 if(-cLPDG==pLPDG) 00995 { 00996 nLPDG=cRPDG; 00997 nRPDG=pRPDG; 00998 } 00999 else 01000 { 01001 if (cRPDG > pLPDG) nLPDG=cRPDG*1000+pLPDG*100+1; 01002 else if(cRPDG < pLPDG) nLPDG=pLPDG*1000+cRPDG*100+1; 01003 else nLPDG=pLPDG*1000+cRPDG*100+3; 01004 if (-cLPDG == pR1) nRPDG=pR2; 01005 else nRPDG=pR1; // -cLPDG == pR2 01006 } 01007 } 01008 else // cRPDG < 0 01009 { 01010 order= 1; 01011 if(-cRPDG==pLPDG) 01012 { 01013 nLPDG=cLPDG; 01014 nRPDG=pRPDG; 01015 } 01016 else 01017 { 01018 if (cLPDG > pLPDG) nLPDG=cLPDG*1000+pLPDG*100+1; 01019 else if(cLPDG < pLPDG) nLPDG=pLPDG*1000+cLPDG*100+1; 01020 else nLPDG=pLPDG*1000+cLPDG*100+3; 01021 if (-cRPDG == pR1) nRPDG=pR2; 01022 else nRPDG=pR1; // -cRPDG == pR2 01023 } 01024 } 01025 break; 01026 case 3: // ....................... pLPDG <-7 01027 if(cLPDG > 0) 01028 { 01029 order= 1; 01030 if(cLPDG==-pRPDG) 01031 { 01032 nLPDG=pLPDG; 01033 nRPDG=cRPDG; 01034 } 01035 else 01036 { 01037 if (cRPDG < pRPDG) nRPDG=cRPDG*1000+pRPDG*100-1; 01038 else if(cRPDG > pRPDG) nRPDG=pRPDG*1000+cRPDG*100-1; 01039 else nRPDG=pRPDG*1000+cRPDG*100-3; 01040 if ( cLPDG == pL1) nLPDG=-pL2; 01041 else nLPDG=-pL1; // cLPDG == pL2 01042 } 01043 } 01044 else // cRPDG > 0 01045 { 01046 order=-1; 01047 if(cRPDG==-pRPDG) 01048 { 01049 nLPDG=pLPDG; 01050 nRPDG=cLPDG; 01051 } 01052 else 01053 { 01054 if (cLPDG < pRPDG) nRPDG=cLPDG*1000+pRPDG*100-1; 01055 else if(cLPDG > pRPDG) nRPDG=pRPDG*1000+cLPDG*100-1; 01056 else nRPDG=pRPDG*1000+cLPDG*100-3; 01057 if ( cRPDG == pL1) nLPDG=-pL2; 01058 else nLPDG=-pL1; // cRPDG == pL2 01059 } 01060 } 01061 break; 01062 case 4: // ....................... pRPDG <-7 01063 if(cLPDG > 0) 01064 { 01065 order=-1; 01066 if(cLPDG==-pLPDG) 01067 { 01068 nLPDG=cRPDG; 01069 nRPDG=pRPDG; 01070 } 01071 else 01072 { 01073 if (cRPDG < pLPDG) nLPDG=cRPDG*1000+pLPDG*100-1; 01074 else if(cRPDG > pLPDG) nLPDG=pLPDG*1000+cRPDG*100-1; 01075 else nLPDG=pLPDG*1000+cRPDG*100-3; 01076 if ( cLPDG == pR1) nRPDG=-pR2; 01077 else nRPDG=-pR1; // cLPDG == pR2 01078 } 01079 } 01080 else // cRPDG > 0 01081 { 01082 order= 1; 01083 if(cRPDG==-pLPDG) 01084 { 01085 nLPDG=cLPDG; 01086 nRPDG=pRPDG; 01087 } 01088 else 01089 { 01090 if (cLPDG < pLPDG) nLPDG=cLPDG*1000+pLPDG*100-1; 01091 else if(cLPDG > pLPDG) nLPDG=pLPDG*1000+cLPDG*100-1; 01092 else nLPDG=pLPDG*1000+cLPDG*100-3; 01093 if ( cRPDG == pR1) nRPDG=-pR2; 01094 else nRPDG=-pR1; // cRPDG == pR2 01095 } 01096 } 01097 break; 01098 } 01099 break; 01100 case 3: // ------------------------> QDiQ/aQaDiQ + QaQ = QDiQ/aQaDiQ (single) 01101 switch (af) 01102 { 01103 case 1: // ....................... cLPDG > 7 01104 if(pLPDG < 0) 01105 { 01106 order= 1; 01107 if (pRPDG > cRPDG) nRPDG=pRPDG*1000+cRPDG*100+1; 01108 else if(pRPDG < cRPDG) nRPDG=cRPDG*1000+pRPDG*100+1; 01109 else nRPDG=cRPDG*1000+pRPDG*100+3; 01110 if (-pLPDG == L1) nLPDG=L2; 01111 else nLPDG=L1; // -pLPDG == L2 01112 } 01113 else // pRPDG < 0 01114 { 01115 order=-1; 01116 if (pLPDG > cRPDG) nLPDG=pLPDG*1000+cRPDG*100+1; 01117 else if(pLPDG < cRPDG) nLPDG=cRPDG*1000+pLPDG*100+1; 01118 else nLPDG=cRPDG*1000+pLPDG*100+3; 01119 if (-pRPDG == L1) nRPDG=L2; 01120 else nRPDG=L1; // -pRPDG == L2 01121 } 01122 break; 01123 case 2: // ....................... cRPDG > 7 01124 if(pLPDG < 0) 01125 { 01126 order=-1; 01127 if (pRPDG > cLPDG) nRPDG=pRPDG*1000+cLPDG*100+1; 01128 else if(pRPDG < cLPDG) nRPDG=cLPDG*1000+pRPDG*100+1; 01129 else nRPDG=cLPDG*1000+pRPDG*100+3; 01130 if (-pLPDG == R1) nLPDG=R2; 01131 else nLPDG=R1; // -pLPDG == R2 01132 } 01133 else // pRPDG < 0 01134 { 01135 order= 1; 01136 if (pLPDG > cLPDG) nLPDG=pLPDG*1000+cLPDG*100+1; 01137 else if(pLPDG < cLPDG) nLPDG=cLPDG*1000+pLPDG*100+1; 01138 else nLPDG=cLPDG*1000+pLPDG*100+3; 01139 if (-pRPDG == R1) nRPDG=R2; 01140 else nRPDG=R1; // -pRPDG == R2 01141 } 01142 break; 01143 case 3: // ....................... cLPDG <-7 (cRPDG <0) 01144 if(pLPDG > 0) 01145 { 01146 order= 1; 01147 if (pRPDG < cRPDG) nRPDG=pRPDG*1000+cRPDG*100-1; 01148 else if(pRPDG > cRPDG) nRPDG=cRPDG*1000+pRPDG*100-1; 01149 else nRPDG=cRPDG*1000+pRPDG*100-3; 01150 if ( pLPDG == L1) nLPDG=-L2; 01151 else nLPDG=-L1; // pLPDG == L2 01152 } 01153 else // pRPDG > 0 01154 { 01155 order=-1; 01156 if (pLPDG < cRPDG) nLPDG=pLPDG*1000+cRPDG*100-1; 01157 else if(pLPDG > cRPDG) nLPDG=cRPDG*1000+pLPDG*100-1; 01158 else nLPDG=cRPDG*1000+pLPDG*100-3; 01159 if ( pRPDG == L1) nRPDG=-L2; 01160 else nRPDG=-L1; // pRPDG == L2 01161 } 01162 break; 01163 case 4: // ....................... cRPDG <-7 (cLPDG <0) 01164 if(pLPDG > 0) // pRPDG & cLPDG are anti-quarks 01165 { 01166 order=-1; 01167 if (pRPDG < cLPDG) nRPDG=pRPDG*1000+cLPDG*100-1; 01168 else if(pRPDG > cLPDG) nRPDG=cLPDG*1000+pRPDG*100-1; 01169 else nRPDG=cLPDG*1000+pRPDG*100-3; 01170 if ( pLPDG == R1) nLPDG=-R2; 01171 else nLPDG=-R1; // pLPDG == R2 01172 } 01173 else // pRPDG > 0 01174 { 01175 order= 1; 01176 if (pLPDG < cLPDG) nLPDG=pLPDG*1000+cLPDG*100-1; 01177 else if(pLPDG > cLPDG) nLPDG=cLPDG*1000+pLPDG*100-1; 01178 else nLPDG=cLPDG*1000+pLPDG*100-3; 01179 if ( pRPDG == R1) nRPDG=-R2; 01180 else nRPDG=-R1; // pRPDG == R2 01181 } 01182 break; 01183 } 01184 break; 01185 case 4: // ------------------------------------> QaQ + aDiQDiQ = QaQ (double) 01186 switch (af) 01187 { 01188 case 1: // ....................... pLPDG > 7 && pRPDG <-7 === LL/RR 01189 order= 1; 01190 if(-cLPDG == pL1) nLPDG= pL2; 01191 else nLPDG= pL1; 01192 if( cRPDG == pR1) nRPDG=-pR2; 01193 else nRPDG=-pR1; 01194 break; 01195 case 2: // ...................... pLPDG > 7 && pRPDG <-7 === LR/RL 01196 order=-1; 01197 if(-cRPDG == pL1) nLPDG= pL2; 01198 else nLPDG= pL1; 01199 if( cLPDG == pR1) nRPDG=-pR2; 01200 else nRPDG=-pR1; 01201 break; 01202 case 3: // ...................... pRPDG > 7 && pLPDG <-7 === LL/RR 01203 order= 1; 01204 if( cLPDG == pL1) nLPDG=-pL2; 01205 else nLPDG=-pL1; 01206 if(-cRPDG == pR1) nRPDG= pR2; 01207 else nRPDG= pR1; 01208 break; 01209 case 4: // ...................... pRPDG > 7 && pLPDG <-7 === LR/RL 01210 order=-1; 01211 if( cRPDG == pL1) nLPDG=-pL2; 01212 else nLPDG=-pL1; 01213 if(-cLPDG == pR1) nRPDG= pR2; 01214 else nRPDG= pR1; 01215 break; 01216 } 01217 break; 01218 case 5: // ------------------------------------> aDiQDiQ + QaQ = QaQ (double) 01219 switch (af) 01220 { 01221 case 1: // ...................... cLPDG > 7 && cRPDG <-7 === LL/RR 01222 order= 1; 01223 if(-pLPDG == L1) nLPDG= L2; 01224 else nLPDG= L1; 01225 if( pRPDG == R1) nRPDG=-R2; 01226 else nRPDG=-R1; 01227 break; 01228 case 2: // ...................... cLPDG > 7 && cRPDG <-7 === LR/RL 01229 order=-1; 01230 if(-pRPDG == L1) nRPDG= L2; 01231 else nRPDG= L1; 01232 if( pLPDG == R1) nLPDG=-R2; 01233 else nLPDG=-R1; 01234 break; 01235 case 3: // ...................... cRPDG > 7 && cLPDG <-7 === LL/RR 01236 order= 1; 01237 if( pLPDG == L1) nLPDG=-L2; 01238 else nLPDG=-L1; 01239 if(-pRPDG == R1) nRPDG= R2; 01240 else nRPDG= R1; 01241 break; 01242 case 4: // ...................... cRPDG > 7 && cLPDG <-7 === LR/RL 01243 order=-1; 01244 if( pRPDG == L1) nRPDG=-L2; 01245 else nRPDG=-L1; 01246 if(-pLPDG == R1) nLPDG= R2; 01247 else nLPDG= R1; 01248 break; 01249 } 01250 break; 01251 case 6: // ------------------------------------> QDiQ + aQaDiQ = QaQ (double) 01252 switch (af) 01253 { 01254 case 1: 01255 order=-1; 01256 if(-cRPDG == pL1) nLPDG= pL2; 01257 else nLPDG= pL1; 01258 if( pRPDG == L1) nRPDG= -L2; 01259 else nRPDG= -L1; 01260 break; 01261 case 2: 01262 order= 1; 01263 if(-cLPDG == pL1) nLPDG= pL2; 01264 else nLPDG= pL1; 01265 if( pRPDG == R1) nRPDG= -R2; 01266 else nRPDG= -R1; 01267 break; 01268 case 3: 01269 order= 1; 01270 if(-cRPDG == pR1) nRPDG= pR2; 01271 else nRPDG= pR1; 01272 if( pLPDG == L1) nLPDG= -L2; 01273 else nLPDG= -L1; 01274 break; 01275 case 4: 01276 order=-1; 01277 if(-cLPDG == pR1) nRPDG= pR2; 01278 else nRPDG= pR1; 01279 if( pLPDG == R1) nLPDG= -R2; 01280 else nLPDG= -R1; 01281 break; 01282 case 5: 01283 order=-1; 01284 if(-pRPDG == L1) nRPDG= L2; 01285 else nRPDG= L1; 01286 if( cRPDG == pL1) nLPDG=-pL2; 01287 else nLPDG=-pL1; 01288 break; 01289 case 6: 01290 order= 1; 01291 if(-pLPDG == L1) nLPDG= L2; 01292 else nLPDG= L1; 01293 if( cRPDG == pR1) nRPDG=-pR2; 01294 else nRPDG=-pR1; 01295 break; 01296 case 7: 01297 order= 1; 01298 if(-pRPDG == R1) nRPDG= R2; 01299 else nRPDG= R1; 01300 if( cLPDG == pL1) nLPDG=-pL2; 01301 else nLPDG=-pL1; 01302 break; 01303 case 8: 01304 order=-1; 01305 if(-pLPDG == R1) nLPDG= R2; 01306 else nLPDG= R1; 01307 if( cLPDG == pR1) nRPDG=-pR2; 01308 else nRPDG=-pR1; 01309 break; 01310 } 01311 break; 01312 } 01313 if(!order) G4cerr<<"-Warning-G4QInel::Constr: t="<<tf<<", a="<<af<<", cL="<<cLPDG 01314 <<", cR="<<cRPDG<<", pL="<<pLPDG<<", pR="<<pRPDG<<G4endl; 01315 else 01316 { 01317 // With theNewHypotheticalPartons the min mass must be calculated & compared 01318 G4int LT=1; 01319 if(std::abs(nLPDG) > 7) ++LT; 01320 G4int RT=1; 01321 if(std::abs(nRPDG) > 7) ++RT; 01322 G4double minM=0.; 01323 G4bool sing=true; 01324 if(cLT==2 && cRT==2) 01325 { 01326 aLPDG=0; 01327 aRPDG=0; 01328 if(cLPDG>0) 01329 { 01330 aLPDG=nLPDG/100; 01331 aRPDG=(-nRPDG)/100; 01332 } 01333 else //cRPDG>0 01334 { 01335 aRPDG=nRPDG/100; 01336 aLPDG=(-nLPDG)/100; 01337 } 01338 G4int nL1=aLPDG/10; 01339 G4int nL2=aLPDG%10; 01340 G4int nR1=aRPDG/10; 01341 G4int nR2=aRPDG%10; 01342 if(nL1!=nR1 && nL1!=nR2 && nL2!=nR1 && nL2!=nR2) // Unreducable DiQ-aDiQ 01343 { 01344 #ifdef debug 01345 G4cout<<"G4QIonIonCollis::Constr:aLPDG="<<aLPDG<<", aRPDG="<<aRPDG<<G4endl; 01346 #endif 01347 sing=false; 01348 G4QPDGCode tmp; 01349 std::pair<G4int,G4int> pB=tmp.MakeTwoBaryons(nL1, nL2, nR1, nR2); 01350 minM=G4QPDGCode(pB.first).GetMass()+G4QPDGCode(pB.second).GetMass(); 01351 } 01352 } 01353 if(sing) 01354 { 01355 std::pair<G4int,G4int> newPair = std::make_pair(nLPDG,nRPDG); 01356 G4QContent newStQC(newPair); // NewString QuarkContent 01357 #ifdef debug 01358 G4cout<<"G4QIn::Con: LPDG="<<nLPDG<<",RPDG="<<nRPDG<<",QC="<<newStQC<<G4endl; 01359 #endif 01360 G4int minPDG=newStQC.GetSPDGCode(); // PDG of the Lightest Hadron=String 01361 minM=G4QPDGCode(minPDG).GetMass() + mPi0; // Min SingleHadron=String Mass 01362 } 01363 // Compare this mass 01364 G4bool win=false; 01365 G4double pSM=0.; 01366 if(pSM2>0.) pSM=std::sqrt(pSM2); 01367 if(minC && pSM2 > maxiM2) // Up to now any positive mass is good 01368 { 01369 maxiM2=pSM2; 01370 win=true; 01371 } 01372 else if(!minC || pSM > minM) 01373 { 01374 pExcess=pSM-minM; 01375 if(minC || pExcess > excess) 01376 { 01377 minC=false; 01378 excess=pExcess; 01379 win=true; 01380 } 01381 } 01382 if(win) 01383 { 01384 sst=pst; 01385 sLPDG=nLPDG; 01386 sRPDG=nRPDG; 01387 sOrd=order; 01388 } 01389 } // End of IF(new partons are created) 01390 } // End of IF(compatible partons) 01391 //} // End of positive squared mass of the fused string 01392 } // End of the LOOP over the possible partners (with the exclusive if for itself) 01393 if(sOrd) // The best pStringCandidate was found 01394 { 01395 G4LorentzVector cL4M=cLeft->Get4Momentum(); 01396 G4LorentzVector cR4M=cRight->Get4Momentum(); 01397 G4QParton* pLeft=(*sst)->GetLeftParton(); 01398 G4QParton* pRight=(*sst)->GetRightParton(); 01399 G4LorentzVector pL4M=pLeft->Get4Momentum(); 01400 G4LorentzVector pR4M=pRight->Get4Momentum(); 01401 #ifdef debug 01402 G4cout<<"G4QIonIonCollis::Constr:cS4M="<<cS4M<<" fused/w pS4M="<<pL4M+pR4M<<G4endl; 01403 #endif 01404 if(sOrd>0) 01405 { 01406 pL4M+=cL4M; 01407 pR4M+=cR4M; 01408 } 01409 else 01410 { 01411 pL4M+=cR4M; 01412 pR4M+=cL4M; 01413 } 01414 pLeft->SetPDGCode(sLPDG); 01415 pLeft->Set4Momentum(pL4M); 01416 pRight->SetPDGCode(sRPDG); 01417 pRight->Set4Momentum(pR4M); 01418 delete (*ist); 01419 strings.erase(ist); 01420 #ifdef debug 01421 G4LorentzVector ss4M=pL4M+pR4M; 01422 G4cout<<"G4QIonIonCollision::Constr:Created,4M="<<ss4M<<",m2="<<ss4M.m2()<<G4endl; 01423 #endif 01424 if( ist != strings.begin() ) // To avoid going below begin() (for Windows) 01425 { 01426 ist--; 01427 con=false; 01428 #ifdef debug 01429 G4cout<<"G4QIonIonCollision::Construct: *** IST Decremented ***"<<G4endl; 01430 #endif 01431 } 01432 else 01433 { 01434 con=true; 01435 #ifdef debug 01436 G4cout<<"G4QIonIonCollision::Construct: *** IST Begin ***"<<G4endl; 01437 #endif 01438 break; 01439 } 01440 } // End of the IF(the best partnerString candidate was found) 01441 else 01442 { 01443 #ifdef debug 01444 G4cout<<"-Warning-G4QInel::Const: S4M="<<cS4M<<",M2="<<cSM2<<" Leave ASIS"<<G4endl; 01445 #endif 01446 ++problem; 01447 con=false; 01448 } 01449 } 01450 else con=false; 01451 } // End of loop over ist iterator 01452 #ifdef debug 01453 G4cout<<"G4QIonIonCollision::Construct: *** IST While *** , con="<<con<<G4endl; 01454 #endif 01455 } // End of "con" while 01456 #ifdef edebug 01457 // This print has meaning only if something appear between it and the StringFragmLOOP 01458 G4LorentzVector t4M=theProjNucleus.Get4Momentum()+theTargNucleus.Get4Momentum();//NucInLS 01459 G4int rC=totChg-theProjNucleus.GetZ()-theTargNucleus.GetZ(); 01460 G4int rB=totBaN-theProjNucleus.GetA()-theTargNucleus.GetA(); 01461 G4int nStr=strings.size(); 01462 G4cout<<"-EMCLS-G4QIn::Const: AfterSUPPRESION #ofS="<<nStr<<",tNuc4M(E=M)="<<sum<<G4endl; 01463 for(G4int i=0; i<nStr; i++) 01464 { 01465 G4LorentzVector strI4M=strings[i]->Get4Momentum(); 01466 t4M+=strI4M; 01467 G4int sChg=strings[i]->GetCharge(); 01468 rC-=sChg; 01469 G4int sBaN=strings[i]->GetBaryonNumber(); 01470 rB-=sBaN; 01471 G4cout<<"-EMCLS-G4QIonIonCollision::Construct: St#"<<i<<", 4M="<<strI4M<<", M=" 01472 <<strI4M.m()<<", C="<<sChg<<", B="<<sBaN<<G4endl; 01473 } 01474 G4cout<<"-EMCLS-G4QInel::Construct: r4M="<<t4M-totLS4M<<", rC="<<rC<<", rB="<<rB<<G4endl; 01475 #endif 01476 // 01477 // --- If a problem is foreseen then the DiQaDiQ strings should be reduced if possible -- 01478 // 01479 #ifdef debug 01480 G4cout<<"G4QIonIonCollision::Construct: problem="<<problem<<G4endl; 01481 #endif 01482 if(problem) 01483 { 01484 G4int nOfStr=strings.size(); 01485 #ifdef debug 01486 G4cout<<"G4QIonIonCollision::Constr:SecurityDiQaDiQReduction, #OfStr="<<nOfStr<<G4endl; 01487 #endif 01488 for (G4int astring=0; astring < nOfStr; astring++) 01489 { 01490 G4QString* curString=strings[astring]; 01491 G4QParton* cLeft=curString->GetLeftParton(); 01492 G4QParton* cRight=curString->GetRightParton(); 01493 G4int LT=cLeft->GetType(); 01494 G4int RT=cRight->GetType(); 01495 G4int sPDG=cLeft->GetPDGCode(); 01496 G4int nPDG=cRight->GetPDGCode(); 01497 if(LT==2 && RT==2) 01498 { 01499 #ifdef debug 01500 G4cout<<"G4QIonIonCollis::Constr:TrySelfReduString,L="<<sPDG<<",R="<<nPDG<<G4endl; 01501 #endif 01502 if( cLeft->ReduceDiQADiQ(cLeft, cRight) ) // DiQ-aDiQ pair was successfully reduced 01503 { 01504 sPDG=cLeft->GetPDGCode(); 01505 nPDG=cRight->GetPDGCode(); 01506 #ifdef debug 01507 G4cout<<"+G4QInel::Const:#"<<astring<<" Reduced, L="<<sPDG<<", R="<<nPDG<<G4endl; 01508 #endif 01509 } 01510 #ifdef debug 01511 else G4cout<<"--*--G4QInel::Const: #"<<astring<<" DQ-aDQ reduction Failed"<<G4endl; 01512 #endif 01513 } // End of the found DiQ/aDiQ pair 01514 else if(sPDG==3 && nPDG==-3) 01515 { 01516 sPDG= 1; 01517 nPDG=-1; 01518 cLeft->SetPDGCode(sPDG); 01519 cRight->SetPDGCode(nPDG); 01520 } 01521 else if(sPDG==-3 && nPDG==3) 01522 { 01523 sPDG=-1; 01524 nPDG= 1; 01525 cLeft->SetPDGCode(sPDG); 01526 cRight->SetPDGCode(nPDG); 01527 } 01528 } 01529 SwapPartons(); 01530 } // End of IF(problem) 01531 #ifdef edebug 01532 G4LorentzVector u4M=theProjNucleus.Get4Momentum()+theTargNucleus.Get4Momentum();//NucInLS 01533 G4int rCh=totChg-theProjNucleus.GetZ()-theTargNucleus.GetZ(); 01534 G4int rBa=totBaN-theProjNucleus.GetA()-theTargNucleus.GetA(); 01535 G4int nStri=strings.size(); 01536 G4cout<<"-EMCLS-G4QIn::Const: FinalConstruct, #ofSt="<<nStri<<",tN4M(E=M)="<<t4M<<G4endl; 01537 for(G4int i=0; i<nStri; i++) 01538 { 01539 G4LorentzVector strI4M=strings[i]->Get4Momentum(); 01540 u4M+=strI4M; 01541 G4int sChg=strings[i]->GetCharge(); 01542 rCh-=sChg; 01543 G4int sBaN=strings[i]->GetBaryonNumber(); 01544 rBa-=sBaN; 01545 G4cout<<"-EMCLS-G4QIonIonCollision::Construct: St#"<<i<<", 4M="<<strI4M<<", M=" 01546 <<strI4M.m()<<", C="<<sChg<<", B="<<sBaN<<G4endl; 01547 } 01548 G4cout<<"-EMCLS-G4QInel::Construct: r4M="<<u4M-totLS4M<<",rC="<<rCh<<",rB="<<rBa<<G4endl; 01549 #endif 01550 }
G4QIonIonCollision::~G4QIonIonCollision | ( | ) |
Definition at line 1552 of file G4QIonIonCollision.cc.
01553 { 01554 std::for_each(strings.begin(), strings.end(), DeleteQString() ); 01555 }
G4int G4QIonIonCollision::AnnihilationOrder | ( | G4int | LS, | |
G4int | MS, | |||
G4int | uP, | |||
G4int | mP, | |||
G4int | sP, | |||
G4int | nP | |||
) | [protected] |
Definition at line 3671 of file G4QIonIonCollision.cc.
References G4cerr, G4cout, and G4endl.
Referenced by Breeder().
03673 {// ^L^^ Curent ^^R^ 03674 G4int Ord=0; 03675 // Curent Partner 03676 if (LS==2 && MPS==2 ) // Fuse 2 Q-aQ strings to DiQ-aDiQ 03677 { 03678 #ifdef debug 03679 G4cout<<"G4QIonIonCollision::AnnihilationOrder:QaQ/QaQ->DiQ-aDiQ, uPDG="<<uPDG<<G4endl; 03680 #endif 03681 if ( (uPDG>0 && sPDG>0 && mPDG<0 && nPDG<0) || 03682 (uPDG<0 && sPDG<0 && mPDG>0 && nPDG>0) ) Ord= 1; // LL/RR 03683 else if( (uPDG>0 && nPDG>0 && mPDG<0 && sPDG<0) || 03684 (uPDG<0 && nPDG<0 && mPDG>0 && sPDG>0) ) Ord=-1; // LR/RL 03685 else G4cerr<<"-Warning-G4QIonIonCollision::AnnihilationOrder: Wrong 22 fusion, L=" 03686 <<uPDG<<",R="<<mPDG<<",cL="<<sPDG<<",cR="<<nPDG<<G4endl; 03687 } 03688 else if ( LS==2 && MPS==3 ) 03689 { 03690 if (uPDG > 7) // Fuse pLDiQ 03691 { 03692 #ifdef debug 03693 G4cout<<"G4QIonIonCollision::AnnihOrder:pLDiQ, sPDG="<<sPDG<<", nPDG="<<nPDG<<G4endl; 03694 #endif 03695 if ( sPDG<0 && (-sPDG==uPDG/1000 || -sPDG==(uPDG/100)%10) ) Ord= 1; // LL/RR 03696 else if( nPDG<0 && (-nPDG==uPDG/1000 || -nPDG==(uPDG/100)%10) ) Ord=-1; // LR/RL 03697 else G4cerr<<"-Warning-G4QIonIonCollision::AnnihilationOrder: Wrong pLDiQ, L="<<uPDG 03698 <<",R="<<mPDG<<",cL="<<sPDG<<",cR="<<nPDG<<G4endl; 03699 } 03700 else if (mPDG > 7) // Fuse pRDiQ 03701 { 03702 #ifdef debug 03703 G4cout<<"G4QIonIonCollision::AnnihOrder:pRDiQ, sPDG="<<sPDG<<", nPDG="<<nPDG<<G4endl; 03704 #endif 03705 if ( sPDG<0 && (-sPDG==mPDG/1000 || -sPDG==(mPDG/100)%10) ) Ord=-1; // LR/RL 03706 else if( nPDG<0 && (-nPDG==mPDG/1000 || -nPDG==(mPDG/100)%10) ) Ord= 1; // LL/RR 03707 else G4cerr<<"-Warning-G4QIonIonCollision::AnnihilationOrder: Wrong pRDiQ, L="<<uPDG 03708 <<",R="<<mPDG<<",cL="<<sPDG<<",cR="<<nPDG<<G4endl; 03709 } 03710 else if (uPDG <-7) // Fuse pLaDiQ 03711 { 03712 #ifdef debug 03713 G4cout<<"G4QIonIonCollision::AnnihOrder:pLaDiQ, sPDG="<<sPDG<<",nPDG="<<nPDG<<G4endl; 03714 #endif 03715 if ( sPDG>0 && (sPDG==(-uPDG)/1000 || sPDG==((-uPDG)/100)%10) ) Ord= 1; // LL/RR 03716 else if( nPDG>0 && (nPDG==(-uPDG)/1000 || nPDG==((-uPDG)/100)%10) ) Ord=-1; // LR/RL 03717 else G4cerr<<"-Warning-G4QIonIonCollision::AnnihilationOrder: Wrong pLaDiQ, L="<<uPDG 03718 <<", R="<<mPDG<<", cL="<<sPDG<<", cR="<<nPDG<<G4endl; 03719 } 03720 else if (mPDG <-7) // Fuse pRaDiQ 03721 { 03722 #ifdef debug 03723 G4cout<<"G4QIonIonCollision::AnnihOrder:pRaDiQ, sPDG="<<sPDG<<",nPDG="<<nPDG<<G4endl; 03724 #endif 03725 if ( sPDG>0 && (sPDG==(-mPDG)/1000 || sPDG==((-mPDG)/100)%10) ) Ord=-1; // LR/RL 03726 else if( nPDG>0 && (nPDG==(-mPDG)/1000 || nPDG==((-mPDG)/100)%10) ) Ord= 1; // LL/RR 03727 else G4cerr<<"-Warning-G4QIonIonCollision::AnnihilationOrder: Wrong pRaDiQ, L="<<uPDG 03728 <<", R="<<mPDG<<", cL="<<sPDG<<", cR="<<nPDG<<G4endl; 03729 } 03730 else if( (sPDG<0 && (-sPDG==mPDG || -sPDG==uPDG) ) || 03731 (nPDG<0 && (-nPDG==mPDG || -nPDG==uPDG) ) ) Ord= 2; // @@ Annihilation fusion 03732 #ifdef debug 03733 else G4cout<<"-Warning-G4QIonIonCollision::AnnihilatOrder:Wrong23StringFusion"<<G4endl; 03734 G4cout<<"G4QIonIonCollision::AnnihilationOrder: Ord="<<Ord<<",sPDG="<<sPDG<<",nPDG=" 03735 <<nPDG<<", uPDG="<<uPDG<<",mPDG="<<mPDG<<G4endl; 03736 #endif 03737 } 03738 else if ( LS==3 && MPS==2 ) 03739 { 03740 if (sPDG > 7) // Fuse cLDiQ 03741 { 03742 #ifdef debug 03743 G4cout<<"G4QIonIonCollision::AnnihOrder:cLDiQ, uPDG="<<uPDG<<", mPDG="<<mPDG<<G4endl; 03744 #endif 03745 if ( uPDG<0 && (-uPDG==sPDG/1000 || -uPDG==(sPDG/100)%10) ) Ord= 1; // LL/RR 03746 else if( mPDG<0 && (-mPDG==sPDG/1000 || -mPDG==(sPDG/100)%10) ) Ord=-1; // LR/RL 03747 else G4cerr<<"-Warning-G4QIonIonCollision::AnnihilationOrder: Wrong cLDiQ, L="<<uPDG 03748 <<",R="<<mPDG<<",cL="<<sPDG<<",cR="<<nPDG<<G4endl; 03749 } 03750 else if (nPDG > 7) // Fuse cRDiQ 03751 { 03752 #ifdef debug 03753 G4cout<<"G4QIonIonCollision::AnnihOrder:cRDiQ, uPDG="<<uPDG<<", mPDG="<<mPDG<<G4endl; 03754 #endif 03755 if ( uPDG<0 && (-uPDG==nPDG/1000 || -uPDG==(nPDG/100)%10) ) Ord=-1; // LR/RL 03756 else if( mPDG<0 && (-mPDG==nPDG/1000 || -mPDG==(nPDG/100)%10) ) Ord= 1; // LL/RR 03757 else G4cerr<<"-Warning-G4QIonIonCollision::AnnihilationOrder: Wrong cRDiQ, L="<<uPDG 03758 <<",R="<<mPDG<<",cL="<<sPDG<<",cR="<<nPDG<<G4endl; 03759 } 03760 else if (sPDG <-7) // Fuse cLaDiQ 03761 { 03762 #ifdef debug 03763 G4cout<<"G4QIonIonCollision::AnnihOrder:cLaDiQ, uPDG="<<uPDG<<",mPDG="<<mPDG<<G4endl; 03764 #endif 03765 if ( uPDG>0 && (uPDG==(-sPDG)/1000 || uPDG==((-sPDG)/100)%10) ) Ord= 1; // LL/RR 03766 else if( mPDG>0 && (mPDG==(-sPDG)/1000 || mPDG==((-sPDG)/100)%10) ) Ord=-1; // LR/RL 03767 else G4cerr<<"-Warning-G4QIonIonCollision::AnnihilationOrder: Wrong cLaDiQ, L="<<uPDG 03768 <<", R="<<mPDG<<", cL="<<sPDG<<", cR="<<nPDG<<G4endl; 03769 } 03770 else if (nPDG <-7) // Fuse cRaDiQ 03771 { 03772 #ifdef debug 03773 G4cout<<"G4QIonIonCollision::AnnihOrder:cRaDiQ, uPDG="<<uPDG<<",mPDG="<<mPDG<<G4endl; 03774 #endif 03775 if ( uPDG>0 && (uPDG==(-nPDG)/1000 || uPDG==((-nPDG)/100)%10) ) Ord=-1; // LR/RL 03776 else if( mPDG>0 && (mPDG==(-nPDG)/1000 || mPDG==((-nPDG)/100)%10) ) Ord= 1; // LL/RR 03777 else G4cerr<<"-Warning-G4QIonIonCollision::AnnihilationOrder: Wrong cRaDiQ, L="<<uPDG 03778 <<", R="<<mPDG<<", cL="<<sPDG<<", cR="<<nPDG<<G4endl; 03779 } 03780 else if( (uPDG<0 && (-uPDG==sPDG || -uPDG==nPDG) ) || 03781 (mPDG<0 && (-mPDG==sPDG || -mPDG==nPDG) ) ) Ord=2; // @@ Annihilation fusion 03782 #ifdef debug 03783 else G4cout<<"-Warning-G4QIonIonCollision::AnnihilatOrder:Wrong32StringFusion"<<G4endl; 03784 G4cout<<"G4QIonIonCollision::AnnihilationOrder: Ord="<<Ord<<",sPDG="<<sPDG<<",nPDG=" 03785 <<nPDG<<", uPDG="<<uPDG<<",mPDG="<<mPDG<<G4endl; 03786 #endif 03787 } 03788 else if ( (LS==2 && MPS==4) || (LS==4 && MPS==2) ) 03789 { 03790 if (uPDG > 7) // Annihilate pLDiQ 03791 { 03792 #ifdef debug 03793 G4cout<<"G4QIonIonCollision::AnnihilOrder:pLDiQ,sPDG="<<sPDG<<",nPDG="<<nPDG<<G4endl; 03794 #endif 03795 if ( sPDG<0 && (-sPDG==uPDG/1000 || -sPDG==(uPDG/100)%10) && 03796 (nPDG==(-mPDG)/1000 || nPDG==((-mPDG)/100)%10) ) Ord= 1; // LL/RR 03797 else if( nPDG<0 && (-nPDG==uPDG/1000 || -nPDG==(uPDG/100)%10) && 03798 (sPDG==(-mPDG)/1000 || sPDG==((-mPDG)/100)%10) ) Ord=-1; // LR/RL 03799 else G4cerr<<"-Warning-G4QIonIonCollision::AnnihilationOrder: Wrong 24 pLDiQ, L=" 03800 <<uPDG<<",R="<<mPDG<<",cL="<<sPDG<<",cR="<<nPDG<<G4endl; 03801 } 03802 else if (mPDG >7) // Annihilate pRDiQ 03803 { 03804 #ifdef debug 03805 G4cout<<"G4QIonIonCollision::AnnihilOrder:PRDiQ,sPDG="<<sPDG<<",nPDG="<<nPDG<<G4endl; 03806 #endif 03807 if ( sPDG<0 && (-sPDG==mPDG/1000 || -sPDG==(mPDG/100)%10) && 03808 (nPDG==(-uPDG)/1000 || nPDG==((-uPDG)/100)%10) ) Ord=-1; // LR/RL 03809 else if( nPDG<0 && (-nPDG==mPDG/1000 || -nPDG==(mPDG/100)%10) && 03810 (sPDG==(-uPDG)/1000 || sPDG==((-uPDG)/100)%10) ) Ord= 1; // LL/RR 03811 else G4cerr<<"-Warning-G4QIonIonCollision::AnnihilationOrder: Wrong 24 pLDiQ, L=" 03812 <<uPDG<<",R="<<mPDG<<",cL="<<sPDG<<",cR="<<nPDG<<G4endl; 03813 } 03814 else if (sPDG > 7) // Annihilate cLDiQ 03815 { 03816 #ifdef debug 03817 G4cout<<"G4QIonIonCollision::AnnihilOrder:cLDiQ,uPDG="<<uPDG<<",mPDG="<<mPDG<<G4endl; 03818 #endif 03819 if ( uPDG<0 && (-uPDG==sPDG/1000 || -uPDG==(sPDG/100)%10) && 03820 (mPDG==(-nPDG)/1000 || mPDG==((-nPDG)/100)%10) ) Ord= 1; // LL/RR 03821 else if( mPDG<0 && (-mPDG==sPDG/1000 || -mPDG==(sPDG/100)%10) && 03822 (uPDG==(-nPDG)/1000 || uPDG==((-nPDG)/100)%10) ) Ord=-1; // LR/RL 03823 else G4cerr<<"-Warning-G4QIonIonCollision::AnnihilationOrder: Wrong 24 cLDiQ, L=" 03824 <<uPDG<<",R="<<mPDG<<",cL="<<sPDG<<",cR="<<nPDG<<G4endl; 03825 } 03826 else if (nPDG > 7) // Annihilate cRDiQ 03827 { 03828 #ifdef debug 03829 G4cout<<"G4QIonIonCollision::AnnihilOrder:cRDiQ,uPDG="<<uPDG<<",mPDG="<<mPDG<<G4endl; 03830 #endif 03831 if ( uPDG<0 && (-uPDG==nPDG/1000 || -uPDG==(nPDG/100)%10) && 03832 (mPDG==(-sPDG)/1000 || mPDG==((-sPDG)/100)%10) ) Ord=-1; // LR/RL 03833 else if( mPDG<0 && (-mPDG==nPDG/1000 || -mPDG==(nPDG/100)%10) && 03834 (uPDG==(-sPDG)/1000 || uPDG==((-sPDG)/100)%10) ) Ord= 1; // LL/RR 03835 else G4cerr<<"-Warning-G4QIonIonCollision::AnnihilationOrder: Wrong 24 cRDiQ, L=" 03836 <<uPDG<<",R="<<mPDG<<",cL="<<sPDG<<",cR="<<nPDG<<G4endl; 03837 } 03838 #ifdef debug 03839 else G4cout<<"-Warning-G4QIonIonCollision::AnnihilOrder: Wrong 24StringFusion"<<G4endl; 03840 G4cout<<"G4QIonIonCollision::AnnihilationOrder: Ord="<<Ord<<",sPDG="<<sPDG<<",nPDG=" 03841 <<nPDG<<", uPDG="<<uPDG<<",mPDG="<<mPDG<<G4endl; 03842 #endif 03843 } 03844 else if ( LS==3 && MPS==3 ) 03845 { 03846 if (uPDG > 7) // Annihilate pLDiQ 03847 { 03848 #ifdef debug 03849 G4cout<<"G4QIonIonCollision::AnnihilOrder:pLDiQ,sPDG="<<sPDG<<",nPDG="<<nPDG<<G4endl; 03850 #endif 03851 if ( sPDG<-7 && (-nPDG==uPDG/1000 || -nPDG==(uPDG/100)%10) && 03852 (mPDG==(-sPDG)/1000 || mPDG==((-sPDG)/100)%10) ) Ord=-1; // LR/RL 03853 else if( nPDG<-7 && (-sPDG==uPDG/1000 || -sPDG==(uPDG/100)%10) && 03854 (mPDG==(-nPDG)/1000 || mPDG==((-nPDG)/100)%10) ) Ord= 1; // LL/RR 03855 else G4cerr<<"-Warning-G4QIonIonCollision::AnnihilationOrder: Wrong 33 pLDiQ, L=" 03856 <<uPDG<<",R="<<mPDG<<",cL="<<sPDG<<",cR="<<nPDG<<G4endl; 03857 } 03858 else if(mPDG > 7) // Annihilate pRDiQ 03859 { 03860 #ifdef debug 03861 G4cout<<"G4QIonIonCollision::AnnihilOrder:pRDiQ,sPDG="<<sPDG<<",nPDG="<<nPDG<<G4endl; 03862 #endif 03863 if ( sPDG<-7 && (-nPDG==mPDG/1000 || -nPDG==(mPDG/100)%10) && 03864 (uPDG==(-sPDG)/1000 || uPDG==((-sPDG)/100)%10) ) Ord= 1; // LL/RR 03865 else if( nPDG<-7 && (-sPDG==mPDG/1000 || -sPDG==(mPDG/100)%10) && 03866 (uPDG==(-nPDG)/1000 || uPDG==((-nPDG)/100)%10) ) Ord=-1; // LR/RL 03867 else G4cerr<<"-Warning-G4QIonIonCollision::AnnihilationOrder: Wrong33pRDiQ, L="<<uPDG 03868 <<",R="<<mPDG<<",cL="<<sPDG<<",cR="<<nPDG<<G4endl; 03869 } 03870 else if(sPDG > 7) // Annihilate cLDiQ 03871 { 03872 #ifdef debug 03873 G4cout<<"G4QIonIonCollision::AnnihilOrder:cLDiQ,uPDG="<<uPDG<<",mPDG="<<mPDG<<G4endl; 03874 #endif 03875 if ( uPDG<-7 && (-mPDG==sPDG/1000 || -mPDG==(sPDG/100)%10) && 03876 (nPDG==(-uPDG)/1000 || nPDG==((-uPDG)/100)%10) ) Ord=-1; // LR/RL 03877 else if( mPDG<-7 && (-uPDG==sPDG/1000 || -uPDG==(sPDG/100)%10) && 03878 (nPDG==(-mPDG)/1000 || nPDG==((-mPDG)/100)%10) ) Ord= 1; // LL/RR 03879 else G4cerr<<"-Warning-G4QIonIonCollision::AnnihilationOrder: Wrong 33 cLDiQ, L=" 03880 <<uPDG<<",R="<<mPDG<<",cL="<<sPDG<<",cR="<<nPDG<<G4endl; 03881 } 03882 else if(nPDG > 7) // Annihilate cRDiQ 03883 { 03884 #ifdef debug 03885 G4cout<<"G4QIonIonCollision::AnnihilOrder:cRDiQ,uPDG="<<uPDG<<",mPDG="<<mPDG<<G4endl; 03886 #endif 03887 if ( uPDG<-7 && (-mPDG==nPDG/1000 || -mPDG==(nPDG/100)%10) && 03888 (nPDG==(-uPDG)/1000 || nPDG==((-uPDG)/100)%10) ) Ord= 1; // LL/RR 03889 else if( mPDG<-7 && (-uPDG==nPDG/1000 || -sPDG==(nPDG/100)%10) && 03890 (sPDG==(-mPDG)/1000 || sPDG==((-mPDG)/100)%10) ) Ord=-1; // LR/RL 03891 else G4cerr<<"-Warning-G4QIonIonCollision::AnnihilationOrder: Wrong 33 cRDiQ, L=" 03892 <<uPDG<<",R="<<mPDG<<",cL="<<sPDG<<",cR="<<nPDG<<G4endl; 03893 } 03894 #ifdef debug 03895 else G4cout<<"-Warning-G4QIonIonCollision::AnnihilOrder: Wrong 33StringFusion"<<G4endl; 03896 G4cout<<"G4QIonIonCollision::AnnihilationOrder: Ord="<<Ord<<",sPDG="<<sPDG<<",nPDG=" 03897 <<nPDG<<", uPDG="<<uPDG<<",mPDG="<<mPDG<<G4endl; 03898 #endif 03899 } 03900 return Ord; 03901 }
void G4QIonIonCollision::Breeder | ( | ) | [protected] |
Definition at line 1791 of file G4QIonIonCollision.cc.
References G4QContent::AddParton(), AnnihilationOrder(), DBL_MAX, G4Quasmon::DecayQHadron(), edebug, FatalException, G4QString::FragmentString(), G4cerr, G4cout, G4endl, G4Exception(), G4QParton::Get4Momentum(), G4QString::Get4Momentum(), G4QHadron::Get4Momentum(), G4QNucleus::GetA(), G4QHadron::GetBaryonNumber(), G4QString::GetBaryonNumber(), G4QHadron::GetCharge(), G4QString::GetCharge(), G4QString::GetDirection(), G4QString::GetLeftParton(), G4QPDGCode::GetMass(), G4QNucleus::GetPDG(), G4QPDGCode::GetPDGCode(), G4QHadron::GetPDGCode(), G4QParton::GetPDGCode(), G4QString::GetQC(), G4QChipolino::GetQPDG1(), G4QChipolino::GetQPDG2(), G4QString::GetRightParton(), G4QContent::GetSPDGCode(), G4QParton::GetType(), G4QNucleus::GetZ(), LT, G4QParton::ReduceDiQADiQ(), ReducePair(), G4QHadron::Set4Momentum(), G4QParton::Set4Momentum(), G4QParton::SetPDGCode(), sqr(), and SumPartonPDG().
Referenced by Fragment().
01792 { // This is the member function, which returns the resulting vector of Hadrons & Quasmons 01793 static const G4double eps = 0.001; // Tolerance in MeV 01794 // 01795 // ------------ At this point the strings are fragmenting to hadrons in LS ------------- 01796 // 01797 #ifdef edebug 01798 G4LorentzVector totLS4M=theProjNucleus.Get4Momentum()+theTargNucleus.Get4Momentum(); //LS 01799 G4int totChg=theProjNucleus.GetZ()+theTargNucleus.GetZ(); 01800 G4int totBaN=theProjNucleus.GetA()+theTargNucleus.GetA(); 01801 G4int nStri=strings.size(); 01802 G4cout<<"-EMCLS-G4QIn::Breed: CHECKRecovery #ofS="<<nStri<<",N4M(E=M)="<<totLS4M<<G4endl; 01803 for(G4int i=0; i<nStri; i++) 01804 { 01805 G4LorentzVector strI4M=strings[i]->Get4Momentum(); 01806 totLS4M+=strI4M; 01807 G4int sChg=strings[i]->GetCharge(); 01808 totChg+=sChg; 01809 G4int sBaN=strings[i]->GetBaryonNumber(); 01810 totBaN+=sBaN; 01811 G4cout<<"-EMCLS-G4QIonIonCollision::Breeder: St#"<<i<<", 4M="<<strI4M<<", M=" 01812 <<strI4M.m()<<", C="<<sChg<<", B="<<sBaN<<G4endl; 01813 } 01814 #endif 01815 G4int nOfStr=strings.size(); 01816 #ifdef debug 01817 G4cout<<"G4QIonIonCollision::Breeder: BeforeFragmentation, #OfStr="<<nOfStr<<G4endl; 01818 #endif 01819 G4LorentzVector ft4M(0.,0.,0.,0.); 01820 G4QContent ftQC(0,0,0,0,0,0); 01821 G4bool ftBad=false; 01822 for(G4int i=0; i < nOfStr; ++i) 01823 { 01824 G4LorentzVector pS4M=strings[i]->Get4Momentum(); // String 4-momentum 01825 ft4M+=pS4M; 01826 G4QContent pSQC=strings[i]->GetQC(); // String Quark Content 01827 ftQC+=pSQC; 01828 if(pS4M.m2() < 0.) ftBad=true; 01829 #ifdef debug 01830 G4cout<<">G4QIonIonCollision::Breed:1stTest,S#"<<i<<",4M="<<pS4M<<",QC="<<pSQC<<G4endl; 01831 #endif 01832 } 01833 if(ftBad) 01834 { 01835 G4Quasmon* stringQuasmon = new G4Quasmon(ftQC, ft4M); 01836 #ifdef debug 01837 G4cout<<"->G4QIonIonCollision::Breed:*TotQ*,QC="<<ftQC<<",4M="<<ft4M<<ft4M.m()<<G4endl; 01838 #endif 01839 theQuasmons.push_back(stringQuasmon); 01840 G4LorentzVector rp4M=theProjNucleus.Get4Momentum(); // Nucleus 4-momentum in LS 01841 G4int rpPDG=theProjNucleus.GetPDG(); 01842 G4QHadron* resPNuc = new G4QHadron(rpPDG,rp4M); 01843 theResult->push_back(resPNuc); // Fill the residual projectile nucleus 01844 G4LorentzVector rt4M=theTargNucleus.Get4Momentum(); // Nucleus 4-momentum in LS 01845 G4int rtPDG=theTargNucleus.GetPDG(); 01846 G4QHadron* resTNuc = new G4QHadron(rtPDG,rt4M); 01847 theResult->push_back(resTNuc); // Fill the residual target nucleus 01848 return; 01849 } 01850 for (G4int astring=0; astring < nOfStr; astring++) 01851 { 01852 #ifdef edebug 01853 G4LorentzVector sum=theProjNucleus.Get4Momentum()+theTargNucleus.Get4Momentum(); //inLS 01854 G4int rChg=totChg-theProjNucleus.GetZ()-theTargNucleus.GetZ(); 01855 G4int rBaN=totBaN-theProjNucleus.GetA()-theTargNucleus.GetA(); 01856 G4int nOfHadr=theResult->size(); 01857 G4cout<<"-EMCLS-G4QIonIonCollision::Breed:#ofSt="<<nOfStr<<",#ofHad="<<nOfHadr<<G4endl; 01858 for(G4int i=astring; i<nOfStr; i++) 01859 { 01860 G4LorentzVector strI4M=strings[i]->Get4Momentum(); 01861 sum+=strI4M; 01862 G4int sChg=strings[i]->GetCharge(); 01863 rChg-=sChg; 01864 G4int sBaN=strings[i]->GetBaryonNumber(); 01865 rBaN-=sBaN; 01866 G4cout<<"-EMCLS-G4QI::Breed:S#"<<i<<",4M="<<strI4M<<",C="<<sChg<<",B="<<sBaN<<G4endl; 01867 } 01868 for(G4int i=0; i<nOfHadr; i++) 01869 { 01870 G4LorentzVector hI4M=(*theResult)[i]->Get4Momentum(); 01871 sum+=hI4M; 01872 G4int hChg=(*theResult)[i]->GetCharge(); 01873 rChg-=hChg; 01874 G4int hBaN=(*theResult)[i]->GetBaryonNumber(); 01875 rBaN-=hBaN; 01876 G4cout<<"-EMCLS-G4QIn::Breed: H#"<<i<<",4M="<<hI4M<<",C="<<hChg<<",B="<<hBaN<<G4endl; 01877 } 01878 G4cout<<"....-EMCLS-G4QInel::Br:r4M="<<sum-totLS4M<<",rC="<<rChg<<",rB="<<rBaN<<G4endl; 01879 #endif 01880 G4QString* curString=strings[astring]; 01881 if(!curString->GetDirection()) continue; // Historic for the dead strings: DoesNotWork 01882 #ifdef edebug 01883 G4int curStrChg = curString->GetCharge(); 01884 G4int curStrBaN = curString->GetBaryonNumber(); 01885 #endif 01886 G4LorentzVector curString4M = curString->Get4Momentum(); 01887 #ifdef debug 01888 G4cout<<"=--=>G4QIonIonCollision::Breeder: String#"<<astring<<",s4M/m="<<curString4M 01889 <<curString4M.m()<<", LPart="<<curString->GetLeftParton()->GetPDGCode() 01890 <<", RPart="<<curString->GetRightParton()->GetPDGCode()<<G4endl; 01891 #endif 01892 G4QHadronVector* theHadrons = 0; // Prototype of theStringFragmentationOUTPUT 01893 theHadrons=curString->FragmentString(true);// !! Fragmenting the String !! 01894 if (!theHadrons) // The string can not be fragmented 01895 { 01896 // First try to correct the diQ-antiDiQ strings, converting them to Q-antiQ 01897 G4QParton* cLeft=curString->GetLeftParton(); 01898 G4QParton* cRight=curString->GetRightParton(); 01899 G4int sPDG=cLeft->GetPDGCode(); 01900 G4int nPDG=cRight->GetPDGCode(); 01901 G4int LT=cLeft->GetType(); 01902 G4int RT=cRight->GetType(); 01903 G4int LS=LT+RT; 01904 if(LT==2 && RT==2) 01905 { 01906 #ifdef debug 01907 G4cout<<"G4QIonIonCollision::Breed:TryReduceString, L="<<sPDG<<",R="<<nPDG<<G4endl; 01908 #endif 01909 if( cLeft->ReduceDiQADiQ(cLeft, cRight) ) // DiQ-aDiQ pair was successfully reduced 01910 { 01911 LT=1; 01912 RT=1; 01913 LS=2; 01914 sPDG=cLeft->GetPDGCode(); 01915 nPDG=cRight->GetPDGCode(); 01916 #ifdef debug 01917 G4cout<<"G4QIonIonCollision::Breed:AfterReduction,L="<<sPDG<<",R="<<nPDG<<G4endl; 01918 #endif 01919 theHadrons=curString->FragmentString(true); 01920 cLeft=curString->GetLeftParton(); 01921 cRight=curString->GetRightParton(); 01922 #ifdef debug 01923 G4cout<<"G4QInel::Breed:L="<<cLeft->Get4Momentum()<<",R="<<cRight->Get4Momentum() 01924 <<G4endl; 01925 #endif 01926 } 01927 #ifdef debug 01928 else G4cout<<"^G4QIonIonCollision::Breed: DQ-aDQ reduction to Q-aQ Failed"<<G4endl; 01929 #endif 01930 } // End of the SelfReduction 01931 #ifdef debug 01932 G4cout<<"G4QIonIonCollision::Breed:AfterRedAttempt, theH="<<theHadrons<<", L4M=" 01933 <<cLeft->Get4Momentum()<<", R4M="<<cRight->Get4Momentum()<<G4endl; 01934 #endif 01935 unsigned next=astring+1; // The next string position 01936 if (!theHadrons) // The string can not be fragmented 01937 { 01938 G4int fusionDONE=0; // StringFusion didn't happen (1=Fuse L+L/R+R, -1=Fuse L+R/R+L) 01939 if(next < strings.size()) // TheString isn't theLastString can fuse 01940 { 01941 G4int fustr=0; // The found partner index (never can be 0) 01942 G4int swap=0; // working interger for swapping parton PDG 01943 G4double Vmin=DBL_MAX; // Prototype of the found Velocity Distance 01944 G4int dPDG=nPDG; 01945 G4int qPDG=sPDG; 01946 if(dPDG<-99 || (dPDG>0&&dPDG<7) || qPDG>99 || (qPDG<0 && qPDG>-7)) 01947 { 01948 swap=qPDG; 01949 qPDG=dPDG; 01950 dPDG=swap; 01951 } 01952 if(dPDG>99) dPDG/=100; 01953 if(qPDG<-99) qPDG=-(-qPDG)/100; 01954 #ifdef debug 01955 G4cout<<"G4QIonIonCollision::Breed:TryFuseStringS, q="<<qPDG<<", a="<<dPDG 01956 <<", n="<<next<<G4endl; 01957 #endif 01958 G4ThreeVector curV=curString4M.vect()/curString4M.e(); 01959 G4int reduce=0; // a#of reduced Q-aQ pairs 01960 G4int restr=0; // To use beyon the LOOP for printing 01961 G4int MPS=0; // PLS for the selected string 01962 for (restr=next; restr < nOfStr; restr++) 01963 { 01964 reduce=0; 01965 G4QString* reString=strings[restr]; 01966 G4QParton* Left=reString->GetLeftParton(); 01967 G4QParton* Right=reString->GetRightParton(); 01968 G4int uPDG=Left->GetPDGCode(); 01969 G4int mPDG=Right->GetPDGCode(); 01970 G4int PLT =Left->GetType(); 01971 G4int PRT =Right->GetType(); 01972 G4int aPDG=mPDG; 01973 G4int rPDG=uPDG; 01974 if(aPDG<-99 || (aPDG>0 && aPDG<7) || rPDG>99 || (rPDG<0 && rPDG>-7)) 01975 { 01976 swap=rPDG; 01977 rPDG=aPDG; 01978 aPDG=swap; 01979 } 01980 if(aPDG > 99) aPDG/=100; 01981 if(rPDG <-99) rPDG=-(-rPDG)/100; 01982 // Try to reduce two DQ-aDQ strings 01983 #ifdef debug 01984 G4cout<<"G4Qnel::Breed: TryReduce #"<<restr<<", q="<<rPDG<<",a="<<aPDG<<G4endl; 01985 #endif 01986 if(LT==2 && RT==2 && PLT==2 && PRT==2) // Have a chance for the reduction 01987 { 01988 G4int cQ1=(-qPDG)/10; 01989 G4int cQ2=(-qPDG)%10; 01990 G4int cA1=dPDG/10; 01991 G4int cA2=dPDG%10; 01992 G4int pQ1=(-rPDG)/10; 01993 G4int pQ2=(-rPDG)%10; 01994 G4int pA1=aPDG/10; 01995 G4int pA2=aPDG%10; 01996 #ifdef debug 01997 G4cout<<"G4QIonIonCollision::Breeder: cQ="<<cQ1<<","<<cQ2<<", cA="<<cA1<<"," 01998 <<cA2<<", pQ="<<pQ1<<","<<pQ2<<", pA="<<pA1<<","<<pA2<<G4endl; 01999 #endif 02000 G4bool iQA = (cA1==pQ1 || cA1==pQ2 || cA2==pQ1 || cA2==pQ2); 02001 G4bool iAQ = (cQ1==pA1 || cQ1==pA2 || cQ2==pA1 || cQ2==pA2); 02002 if(iQA) reduce++; 02003 if(iAQ) reduce++; 02004 if (reduce==2) // Two quark pairs can be reduced 02005 { 02006 if(sPDG>0 && uPDG<0) // LL/RR Reduction 02007 { 02008 std::pair<G4int,G4int> resLL=ReducePair(sPDG/100, (-uPDG)/100); 02009 G4int newCL=resLL.first; 02010 G4int newPL=resLL.second; 02011 if(!newCL || !newPL) 02012 { 02013 G4cerr<<"*G4QIonIonCollision::Breed:CL="<<newCL<<",PL="<<newPL<<G4endl; 02014 G4Exception("G4QIonIonCollision::Breeder:","72",FatalException,"2LL-"); 02015 } 02016 std::pair<G4int,G4int> resRR=ReducePair((-nPDG)/100, mPDG/100); 02017 G4int newCR=resRR.first; 02018 G4int newPR=resRR.second; 02019 if(!newCR || !newPR) 02020 { 02021 G4cerr<<"*G4QIonIonCollision::Breed:CR="<<newCR<<",PR="<<newPR<<G4endl; 02022 G4Exception("G4QIonIonCollision::Breeder:","72",FatalException,"2RR-"); 02023 } 02024 cLeft->SetPDGCode(newCL); // Reset the left quark of curString 02025 cRight->SetPDGCode(-newCR); // Reset the right quark of curString 02026 Left->SetPDGCode(-newPL); // Reset the left quark of reString 02027 Right->SetPDGCode(newPR); // Reset the right quark of reString 02028 break; // Break out of the reString internal LOOP 02029 } 02030 else if(sPDG<0 && uPDG>0) // LL/RR Reduction 02031 { 02032 std::pair<G4int,G4int> resLL=ReducePair((-sPDG)/100, uPDG/100); 02033 G4int newCL=resLL.first; 02034 G4int newPL=resLL.second; 02035 if(!newCL || !newPL) 02036 { 02037 G4cerr<<"*G4QIonIonCollision::Breed:CL="<<newCL<<",PL="<<newPL<<G4endl; 02038 G4Exception("G4QIonIonCollision::Breeder:","72",FatalException,"2LL+"); 02039 } 02040 std::pair<G4int,G4int> resRR=ReducePair(nPDG/100, (-mPDG)/100); 02041 G4int newCR=resRR.first; 02042 G4int newPR=resRR.second; 02043 if(!newCR || !newPR) 02044 { 02045 G4cerr<<"*G4QIonIonCollision::Breed:CR="<<newCR<<",PR="<<newPR<<G4endl; 02046 G4Exception("G4QIonIonCollision::Breeder:","72",FatalException,"2RR+"); 02047 } 02048 cLeft->SetPDGCode(-newCL); // Reset the left quark of curString 02049 cRight->SetPDGCode(newCR); // Reset the right quark of curString 02050 Left->SetPDGCode(newPL); // Reset the left quark of reString 02051 Right->SetPDGCode(-newPR); // Reset the right quark of reString 02052 break; // Break out of the reString internal LOOP 02053 } 02054 else if(sPDG>0 && mPDG<0) // LR Reduction 02055 { 02056 std::pair<G4int,G4int> resLL=ReducePair(sPDG/100, (-mPDG)/100); 02057 G4int newCL=resLL.first; 02058 G4int newPR=resLL.second; 02059 if(!newCL || !newPR) 02060 { 02061 G4cerr<<"*G4QIonIonCollision::Breed:CL="<<newCL<<",PR="<<newPR<<G4endl; 02062 G4Exception("G4QIonIonCollision::Breeder:","72",FatalException,"2-LR"); 02063 } 02064 std::pair<G4int,G4int> resRR=ReducePair((-nPDG)/100, uPDG/100); 02065 G4int newCR=resRR.first; 02066 G4int newPL=resRR.second; 02067 if(!newCR || !newPR) 02068 { 02069 G4cerr<<"*G4QIonIonCollision::Breed:CR="<<newCR<<",PL="<<newPL<<G4endl; 02070 G4Exception("G4QIonIonCollision::Breeder:","72",FatalException,"2-LR"); 02071 } 02072 cLeft->SetPDGCode(newCL); // Reset the left quark of curString 02073 cRight->SetPDGCode(-newCR); // Reset the right quark of curString 02074 Left->SetPDGCode(newPL); // Reset the left quark of reString 02075 Right->SetPDGCode(-newPR); // Reset the right quark of reString 02076 break; // Break out of the reString internal LOOP 02077 } 02078 else // RL Reduction 02079 { 02080 std::pair<G4int,G4int> resLL=ReducePair((-sPDG)/100, mPDG/100); 02081 G4int newCL=resLL.first; 02082 G4int newPR=resLL.second; 02083 if(!newCL || !newPR) 02084 { 02085 G4cerr<<"*G4QIonIonCollision::Breed:CL="<<newCL<<",PR="<<newPR<<G4endl; 02086 G4Exception("G4QIonIonCollision::Breeder:","72",FatalException,"2-RL"); 02087 } 02088 std::pair<G4int,G4int> resRR=ReducePair(nPDG/100, (-uPDG)/100); 02089 G4int newCR=resRR.first; 02090 G4int newPL=resRR.second; 02091 if(!newCR || !newPR) 02092 { 02093 G4cerr<<"*G4QIonIonCollision::Breed:CR="<<newCR<<",PL="<<newPL<<G4endl; 02094 G4Exception("G4QIonIonCollision::Breeder:","72",FatalException,"2-RL"); 02095 } 02096 cLeft->SetPDGCode(-newCL); // Reset the left quark of curString 02097 cRight->SetPDGCode(newCR); // Reset the right quark of curString 02098 Left->SetPDGCode(-newPL); // Reset the left quark of reString 02099 Right->SetPDGCode(newPR); // Reset the right quark of reString 02100 break; // Break out of the reString internal LOOP 02101 } 02102 } // End of IF(possible reduction) 02103 } 02104 // Check StringsCanBeFused: all QaQ+QaQ(22), single QaQ+QDiQ/aQaDtQ(23/32), 02105 // double QaQ+DiQaDiQ(24/42), QDiQ+aDiQaQ(34/43) 02106 #ifdef debug 02107 G4cout<<"G4QInel::Breed: TryFuse/w #"<<restr<<",q="<<rPDG<<",a="<<aPDG<<G4endl; 02108 #endif 02109 G4int PLS=PLT+PRT; 02110 if( (LS==2 && PLS==2) || // QaQ+QaQ always to DiQaDiQ 02111 ( ( (LS==2 && PLS==3) || (LS==3 && PLS==2) ) &&// QaQ w QDiQ/aQaDiQ(single) 02112 ( (aPDG> 7 && (-dPDG==aPDG/10 || -dPDG==aPDG%10) ) || // cAQ w DQ 02113 (dPDG> 7 && (-aPDG==dPDG/10 || -aPDG==dPDG%10) ) || // AQ w cDQ 02114 (rPDG<-7 && (qPDG==(-rPDG)/10 || qPDG==(-rPDG)%10) ) || // cQ w ADQ 02115 (qPDG<-7 && (rPDG==(-qPDG)/10 || rPDG==(-qPDG)%10) ) // Q w cADQ 02116 //|| (aPDG< 0 && -aPDG==qPDG) || (dPDG< 0 && -dPDG==rPDG) // aQ w Q 02117 ) 02118 ) || // ----------------------- 02119 ( ( LS==2 && PLS==4 &&// QaQ w DiQaDiQ (double) 02120 (aPDG> 7 && (-dPDG == aPDG/10 || -dPDG == aPDG%10) ) && 02121 (rPDG<-7 && (qPDG==(-rPDG)/10 || qPDG==(-rPDG)%10) ) 02122 ) || 02123 ( LS==4 && PLS==2 &&// DiQaDiQ w QaQ (double) 02124 (dPDG> 7 && (-aPDG == dPDG/10 || -aPDG == dPDG%10) ) && 02125 (qPDG<-7 && (rPDG==(-qPDG)/10 || rPDG==(-qPDG)%10) ) 02126 ) 02127 ) || // ----------------------- 02128 ( LS==3 && PLS==3 &&// QDiQ w aDiQaQ (double) 02129 ( (aPDG> 7 && (-dPDG == aPDG/10 || -dPDG == aPDG%10) && 02130 qPDG<-7 && (rPDG==(-qPDG)/10 || rPDG==(-qPDG)%10) 02131 ) || 02132 (dPDG> 7 && (-aPDG == dPDG/10 || -aPDG == dPDG%10) && 02133 rPDG<-7 && (dPDG==(-rPDG)/10 || dPDG==(-rPDG)%10) 02134 ) 02135 ) 02136 ) 02137 ) 02138 { 02139 G4LorentzVector reString4M = reString->Get4Momentum(); 02140 G4ThreeVector reV = reString4M.vect()/reString4M.e(); 02141 G4double dV=(curV-reV).mag2(); // Squared difference between relVelocities 02142 #ifdef debug 02143 G4cout<<"G4QIonIonCollision::Breeder: StringCand#"<<restr<<", q="<<rPDG 02144 <<", a="<<aPDG<<", L="<<uPDG<<", R="<<mPDG<<",dV="<<dV<<G4endl; 02145 #endif 02146 if(dV < Vmin) 02147 { 02148 Vmin=dV; 02149 fustr=restr; 02150 MPS=PLS; 02151 } 02152 } 02153 } 02154 if(reduce==2) // String mutual reduction happened 02155 { 02156 #ifdef debug 02157 G4cout<<"-G4QIonIonCollision::Breed:Reduced #"<<astring<<" & #"<<restr<<G4endl; 02158 #endif 02159 astring--; // String was QCreduced using another String 02160 continue; // Repeat fragmentation of the same string 02161 } 02162 if(fustr) // The partner was found -> fuse strings 02163 { 02164 #ifdef debug 02165 G4cout<<"G4QInel::Breeder: StPartner#"<<fustr<<", LT="<<LT<<",RT="<<RT<<G4endl; 02166 #endif 02167 G4QString* fuString=strings[fustr]; 02168 G4QParton* Left=fuString->GetLeftParton(); 02169 G4QParton* Right=fuString->GetRightParton(); 02170 G4int uPDG=Left->GetPDGCode(); 02171 G4int mPDG=Right->GetPDGCode(); 02172 G4int rPDG=uPDG; 02173 G4int aPDG=mPDG; 02174 if(aPDG<-99 || (aPDG>0 && aPDG<7) || rPDG>99 || (rPDG<0 && rPDG>-7)) 02175 { 02176 swap=rPDG; 02177 rPDG=aPDG; 02178 aPDG=swap; 02179 } 02180 if(aPDG > 99) aPDG/=100; 02181 if(rPDG <-99) rPDG=-(-rPDG)/100; 02182 // Check that the strings can fuse (no anti-diquarks assumed) 02183 G4LorentzVector resL4M(0.,0.,0.,0.); 02184 G4LorentzVector resR4M(0.,0.,0.,0.); 02185 G4LorentzVector L4M=cLeft->Get4Momentum(); 02186 G4LorentzVector R4M=cRight->Get4Momentum(); 02187 #ifdef debug 02188 G4cout<<"G4QIonIonCollision::Breeder:BeforeFuDir,sL="<<sPDG<<",nR="<<nPDG 02189 <<",uL="<<uPDG<<",mR="<<mPDG<<",L4M="<<L4M<<",R4M="<<R4M<<G4endl; 02190 #endif 02191 G4LorentzVector PL4M=Left->Get4Momentum(); 02192 G4LorentzVector PR4M=Right->Get4Momentum(); 02193 fusionDONE=AnnihilationOrder(LS, MPS, uPDG, mPDG, sPDG, nPDG); 02194 if (fusionDONE>0) 02195 { 02196 if(fusionDONE>1) // Anihilation algorithm 02197 { 02198 if ( (uPDG<0 || nPDG<0) && -uPDG==nPDG ) Left->SetPDGCode(sPDG); 02199 else if( (mPDG<0 || sPDG<0) && -mPDG==sPDG ) Right->SetPDGCode(nPDG); 02200 else if( (uPDG<0 || sPDG<0) && -uPDG==sPDG ) Left->SetPDGCode(nPDG); 02201 else if( (mPDG<0 || nPDG<0) && -mPDG==nPDG ) Right->SetPDGCode(sPDG); 02202 } 02203 { 02204 Left->SetPDGCode(SumPartonPDG(uPDG, sPDG)); 02205 Right->SetPDGCode(SumPartonPDG(mPDG, nPDG)); 02206 } 02207 Left->Set4Momentum(L4M+PL4M); 02208 Right->Set4Momentum(R4M+PR4M); 02209 #ifdef debug 02210 G4cout<<"G4QIonIonCollision::Breeder:LL/RR s4M="<<fuString->Get4Momentum() 02211 <<",S="<<L4M+PL4M+R4M+PR4M<<", L="<<Left->Get4Momentum()<<", R=" 02212 <<Right->Get4Momentum()<<G4endl; 02213 #endif 02214 } 02215 else if(fusionDONE<0) 02216 { 02217 Left->SetPDGCode(SumPartonPDG(uPDG, nPDG)); 02218 Left->Set4Momentum(L4M+PR4M); 02219 Right->SetPDGCode(SumPartonPDG(mPDG, sPDG)); 02220 Right->Set4Momentum(R4M+PL4M); 02221 #ifdef debug 02222 G4cout<<"G4QIonIonCollision::Breeder:LR/RL s4M="<<fuString->Get4Momentum() 02223 <<",S="<<L4M+PL4M+R4M+PR4M<<", L="<<Left->Get4Momentum()<<", R=" 02224 <<Right->Get4Momentum()<<G4endl; 02225 #endif 02226 } 02227 #ifdef debug 02228 else G4cout<<"-Warning-G4QIonIonCollision::Breeder: WrongStringFusion"<<G4endl; 02229 #endif 02230 #ifdef edebug 02231 G4cout<<"#EMC#G4QIonIonCollision::Breed:StringFused,F="<<fusionDONE<<",L="<<L4M 02232 <<",R="<<R4M<<",pL="<<PL4M<<",pR="<<PR4M<<",nL="<<Left->Get4Momentum() 02233 <<",nR="<<Right->Get4Momentum()<<",S="<<fuString->Get4Momentum()<<G4endl; 02234 #endif 02235 if(fusionDONE) 02236 { 02237 #ifdef debug 02238 G4cout<<"###G4QIonIonCollision::Breed: Str#"<<astring<<" fused/w Str#"<<fustr 02239 <<"->"<<fuString->Get4Momentum()<<fuString->Get4Momentum().m()<<G4endl; 02240 #endif 02241 continue; // @@ killing of the curString is excluded 02242 } 02243 } 02244 #ifdef debug 02245 else 02246 { 02247 02248 G4cerr<<"**G4QIonIonCollision::Breed:*NoPart*M="<<curString->Get4Momentum().m() 02249 <<", F="<<fusionDONE<<", LPDG="<<curString->GetLeftParton()->GetPDGCode() 02250 <<", RPDG="<<curString->GetRightParton()->GetPDGCode()<<G4endl; 02251 } 02252 #endif 02253 } 02254 if(!fusionDONE) // Fusion of theString failed, try hadrons 02255 { 02256 G4int nHadr=theResult->size(); // #of collected resulting hadrons upToNow 02257 #ifdef debug 02258 G4cout<<"+++4QFragmentation::Breeder:*TryHadr* M="<<curString->Get4Momentum().m() 02259 <<", nH="<<nHadr<<", LPDG="<<curString->GetLeftParton()->GetPDGCode() 02260 <<", RPDG="<<curString->GetRightParton()->GetPDGCode()<<G4endl; 02261 #endif 02262 // The only solution now is to try fusion with the secondary hadron 02263 while( (nHadr=theResult->size()) && !theHadrons) 02264 { 02265 #ifdef edebug 02266 for(G4int i=0; i<nHadr; i++) 02267 { 02268 G4LorentzVector h4M=(*theResult)[i]->Get4Momentum(); 02269 G4int hPDG=(*theResult)[i]->GetPDGCode(); 02270 G4QContent hQC=(*theResult)[i]->GetQC(); 02271 G4cout<<"-EMC-G4QInel::Breed:H#"<<i<<",4M="<<h4M<<hQC<<",PDG="<<hPDG<<G4endl; 02272 } 02273 #endif 02274 G4int fusDONE=0; // String+Hadron fusion didn't happen 02275 G4int fuhad=-1; // The found hadron index 02276 G4int newPDG=0; // PDG ofTheParton afterMerging with Hadr 02277 G4int secPDG=0; // Second PDG oParton afterMerging w/Hadr 02278 G4double maM2=-DBL_MAX; // Prototype of the max ResultingStringM2 02279 G4LorentzVector selH4M(0.,0.,0.,0.); // 4-mom of the selected hadron 02280 G4QHadron* selHP=0; // Pointer to the used hadron for erasing 02281 G4QString* cString=strings[astring]; // Must be the last string by definition 02282 G4LorentzVector cString4M = cString->Get4Momentum(); 02283 cLeft=cString->GetLeftParton(); 02284 cRight=cString->GetRightParton(); 02285 G4int sumT=cLeft->GetType()+cRight->GetType(); 02286 sPDG=cLeft->GetPDGCode(); 02287 nPDG=cRight->GetPDGCode(); 02288 G4int cMax=0; // Both or only one cuark can merge 02289 for (G4int reh=0; reh < nHadr; reh++) 02290 { 02291 G4QHadron* curHadr=(*theResult)[reh]; 02292 G4QContent curQC=curHadr->GetQC(); // Quark content of the hadron 02293 G4int partPDG1=curQC.AddParton(sPDG); 02294 G4int partPDG2=curQC.AddParton(nPDG); 02295 #ifdef debug 02296 G4cout<<"G4QIonIonCollision::Breeder: Hadron # "<<reh<<", QC="<<curQC 02297 <<", P1="<<partPDG1<<", P2="<<partPDG2<<G4endl; 02298 #endif 02299 if(partPDG1 || partPDG2) // Hadron can merge at least w/one parton 02300 { 02301 G4int cCur=1; 02302 if(sumT>3 && partPDG1 && partPDG2) cCur=2; 02303 G4LorentzVector curHadr4M = curHadr->Get4Momentum(); 02304 G4double M2=(cString4M+curHadr4M).m2();// SquaredMass of theResultingString 02305 #ifdef debug 02306 G4cout<<"G4QIonIonCollision::Breeder:*IN*Hadron#"<<reh<<",M2="<<M2<<G4endl; 02307 #endif 02308 if( (sumT<4 || cCur>=cMax) && M2 > maM2) // DeepAnnihilation for DiQ-aDiQ 02309 { 02310 maM2=M2; 02311 fuhad=reh; 02312 if(partPDG1) 02313 { 02314 fusDONE= 1; 02315 newPDG=partPDG1; 02316 if(partPDG2) 02317 { 02318 secPDG=partPDG2; 02319 cMax=cCur; 02320 } 02321 } 02322 else 02323 { 02324 fusDONE=-1; 02325 newPDG=partPDG2; 02326 } 02327 #ifdef debug 02328 G4cout<<"G4QInel::Br:*Selected*,P1="<<partPDG1<<",P2="<<partPDG2<<G4endl; 02329 #endif 02330 selH4M=curHadr4M; 02331 selHP=curHadr; 02332 } // End of IF(update selection) 02333 } // End of IF(HadronCanMergeWithTheString) 02334 } // End of the LOOP over Hadrons 02335 #ifdef debug 02336 G4cout<<"G4QIonIonCollision::Breeder: fuh="<<fuhad<<",fus="<<fusDONE<<G4endl; 02337 #endif 02338 if(fuhad>-1) // The partner was found -> fuse H&S 02339 { 02340 if (fusDONE>0) // Fuse hadron with the left parton 02341 { 02342 cLeft->SetPDGCode(newPDG); 02343 G4LorentzVector newLeft=cLeft->Get4Momentum()+selH4M; 02344 cLeft->Set4Momentum(newLeft); 02345 if(secPDG && cMax>1) 02346 { 02347 #ifdef debug 02348 G4cout<<"G4QInel::Br:TryReduce, nPDG="<<newPDG<<",sPDG="<<secPDG<<G4endl; 02349 #endif 02350 cLeft->ReduceDiQADiQ(cLeft, cRight); 02351 } 02352 #ifdef debug 02353 G4cout<<"G4QIonIonCollision::Breed: Left, s4M="<<curString->Get4Momentum() 02354 <<", L4M="<<newLeft<<", R4M="<<cRight->Get4Momentum()<<", h4M=" 02355 <<selH4M<<", newL="<<newPDG<<", oldR="<<cRight->GetPDGCode()<<G4endl; 02356 #endif 02357 } 02358 else if(fusDONE<0) // Fuse hadron with the right parton 02359 { 02360 cRight->SetPDGCode(newPDG); 02361 G4LorentzVector newRight=cRight->Get4Momentum()+selH4M; 02362 cRight->Set4Momentum(newRight); 02363 #ifdef debug 02364 G4cout<<"G4QIonIonCollision::Breed: Right, s4M="<<curString->Get4Momentum() 02365 <<", L4M="<<cLeft->Get4Momentum()<<", R4M="<<newRight<<", h4M=" 02366 <<selH4M<<", newR="<<newPDG<<", oldL="<<cLeft->GetPDGCode()<<G4endl; 02367 #endif 02368 } 02369 #ifdef debug 02370 else G4cout<<"-G4QIonIonCollision::Breed: Wrong String+HadronFusion"<<G4endl; 02371 #endif 02372 #ifdef debug 02373 if(fusDONE) G4cout<<"####G4QIonIonCollision::Breeder: String #"<<astring 02374 <<" is fused with Hadron #"<<fuhad 02375 <<", new4Mom="<<curString->Get4Momentum() 02376 <<", M2="<<curString->Get4Momentum().m2() 02377 <<", QC="<<curString->GetQC()<<G4endl; 02378 #endif 02379 } 02380 else 02381 { 02382 #ifdef debug 02383 G4cerr<<"**G4QIonIonCollision::Breed:*NoH*,M="<<curString->Get4Momentum().m() 02384 <<", LPDG="<<curString->GetLeftParton()->GetPDGCode() 02385 <<", RPDG="<<curString->GetRightParton()->GetPDGCode()<<G4endl; 02386 // @@ Temporary exception for the future solution 02387 //G4Exception("G4QIonIonCollision::Breed:","72",FatalException,"SHNotFused"); 02388 #endif 02389 break; // Breake the While LOOP 02390 } // End of the namespace where both Fusion and reduction have failed 02391 // The fused hadron must be excluded from theResult 02392 #ifdef debug 02393 G4cout<<"G4QIonIonCollision::Breed: before HR, nH="<<theResult->size()<<G4endl; 02394 G4int icon=0; // Loop counter 02395 #endif 02396 G4QHadronVector::iterator ih; 02397 #ifdef debug 02398 G4bool found=false; 02399 #endif 02400 for(ih = theResult->begin(); ih != theResult->end(); ih++) 02401 { 02402 #ifdef debug 02403 G4cout<<"G4QInelast::Breeder:#"<<icon<<", i="<<(*ih)<<", sH="<<selHP<<G4endl; 02404 icon++; 02405 #endif 02406 if((*ih)==selHP) 02407 { 02408 #ifdef debug 02409 G4cout<<"G4QInel::Breed: *HadronFound*, PDG="<<selHP->GetPDGCode()<<G4endl; 02410 #endif 02411 G4LorentzVector p4M=selHP->Get4Momentum(); 02412 curString4M+=p4M; 02413 #ifdef edebug 02414 G4int Chg=selHP->GetCharge(); 02415 G4int BaN=selHP->GetBaryonNumber(); 02416 curStrChg+=Chg; 02417 curStrBaN+=BaN; 02418 G4cout<<"-EMC->>>G4QIonIonCollision::Breed: S+=H, 4M="<<curString4M<<", M=" 02419 <<curString4M.m()<<", Charge="<<curStrChg<<", B="<<curStrBaN<<G4endl; 02420 #endif 02421 delete selHP; // delete the Hadron 02422 theResult->erase(ih); // erase the Hadron from theResult 02423 #ifdef debug 02424 found=true; 02425 #endif 02426 break; // beak the LOOP over hadrons 02427 } 02428 } // End of the LOOP over hadrons 02429 #ifdef debug 02430 if(!found) G4cout<<"*G4QIonIonCollision::Breed:nH="<<theResult->size()<<G4endl; 02431 #endif 02432 // New attempt of the string decay 02433 theHadrons=curString->FragmentString(true); 02434 #ifdef debug 02435 G4cout<<"G4QInel::Breeder: tH="<<theHadrons<<",nH="<<theResult->size()<<G4endl; 02436 #endif 02437 } // End of the while LOOP over the fusion with hadrons 02438 #ifdef debug 02439 G4cout<<"*G4QIonIonCollision::Breed: *CanTryToDecay?* nH="<<theHadrons<<", next=" 02440 <<next<<" =? nS="<<strings.size()<<", nR="<<theResult->size()<<G4endl; 02441 #endif 02442 if(!theHadrons && next == strings.size() && !(theResult->size()))// TryToDecay 02443 { 02444 G4QContent miQC=curString->GetQC(); // QContent of the Lightest Hadron 02445 G4int miPDG=miQC.GetSPDGCode(); // PDG of the Lightest Hadron 02446 if(miPDG == 10) // ==> Decay Hadron-Chipolino 02447 { 02448 G4QChipolino QCh(miQC); // define theTotalResidual as aChipolino 02449 G4QPDGCode h1QPDG=QCh.GetQPDG1(); // QPDG of the first hadron 02450 G4double h1M =h1QPDG.GetMass();// Mass of the first hadron 02451 G4QPDGCode h2QPDG=QCh.GetQPDG2(); // QPDG of the second hadron 02452 G4double h2M =h2QPDG.GetMass();// Mass of the second hadron 02453 G4double ttM =curString4M.m(); // Real Mass of the Chipolino 02454 if(h1M+h2M<ttM+eps) // Two particles decay of Chipolino 02455 { 02456 G4LorentzVector h14M(0.,0.,0.,h1M); 02457 G4LorentzVector h24M(0.,0.,0.,h2M); 02458 if(std::fabs(ttM-h1M-h2M)<=eps) 02459 { 02460 G4double part1=h1M/(h1M+h2M); 02461 h14M=part1*curString4M; 02462 h24M=curString4M-h14M; 02463 } 02464 else 02465 { 02466 if(!G4QHadron(curString4M).DecayIn2(h14M,h24M)) 02467 { 02468 G4cerr<<"***G4QIonIonCollision::Breeder: tM="<<ttM<<"->h1="<<h1QPDG 02469 <<"("<<h1M<<")+h2="<<h1QPDG<<"("<<h2M<<")="<<h1M+h2M<<G4endl; 02470 G4Exception("G4QIonIonCollision::Breeder:","72",FatalException,"QDec"); 02471 } 02472 } 02473 G4QHadron* h1H = new G4QHadron(h1QPDG.GetPDGCode(),h14M); 02474 theResult->push_back(h1H); // (delete equivalent) 02475 #ifdef debug 02476 G4LorentzVector f4M=h1H->Get4Momentum(); 02477 G4int fPD=h1H->GetPDGCode(); 02478 G4int fCg=h1H->GetCharge(); 02479 G4int fBN=h1H->GetBaryonNumber(); 02480 G4cout<<"-EMC->>G4QIonIonCollision::Breed:String=HadrChiPro1's filled,f4M=" 02481 <<f4M<<", fPDG="<<fPD<<", fCg="<<fCg<<", fBN="<<fBN<<G4endl; 02482 #endif 02483 G4QHadron* h2H = new G4QHadron(h2QPDG.GetPDGCode(),h24M); 02484 theResult->push_back(h2H); // (delete equivalent) 02485 #ifdef debug 02486 G4LorentzVector s4M=h2H->Get4Momentum(); 02487 G4int sPD=h2H->GetPDGCode(); 02488 G4int sCg=h2H->GetCharge(); 02489 G4int sBN=h2H->GetBaryonNumber(); 02490 G4cout<<"-EMC->>G4QIonIonCollision::Breed:String=HadrChiPro2's filled,s4M=" 02491 <<s4M<<", sPDG="<<sPD<<", sCg="<<sCg<<", sBN="<<sBN<<G4endl; 02492 #endif 02493 #ifdef edebug 02494 G4cout<<"-EMC-..Chi..G4QIonIonCollision::Breeder: DecayCHECK, Ch4M=" 02495 <<curString4M<<", d4M="<<curString4M-h14M-h24M<<G4endl; 02496 #endif 02497 break; // Go out of the main StringDecay LOOP 02498 } 02499 else 02500 { 02501 G4Quasmon* stringQuasmon = new G4Quasmon(miQC, curString4M);// String->Quas 02502 theQuasmons.push_back(stringQuasmon); 02503 break; // Go out of the main StringDecay LOOP 02504 } 02505 } 02506 else if(miPDG) // Decay Hadron as a Resonans 02507 { 02508 if (miPDG>0 && miPDG%10 < 3) miPDG+=2; // Convert to Resonans 02509 else if(miPDG<0 && (-miPDG)%10< 3) miPDG-=2; // Convert to antiResonans 02510 G4Quasmon Quasm; 02511 G4QHadron* sHad = new G4QHadron(miPDG,curString4M); 02512 G4QHadronVector* tmpQHadVec=Quasm.DecayQHadron(sHad); // It deletes sHad 02513 G4int tmpN=tmpQHadVec->size(); 02514 #ifdef debug 02515 G4cout<<"G4QIonIonCollision::Breeder: Decay the Last, Res#H="<<tmpN<<G4endl; 02516 #endif 02517 if(tmpN>1) 02518 { 02519 for(G4int aH=0; aH < tmpN; aH++) 02520 { 02521 theResult->push_back((*tmpQHadVec)[aH]);//TheDecayProdOfHadronDecIsFilled 02522 #ifdef debug 02523 G4QHadron* prodH =(*tmpQHadVec)[aH]; 02524 G4LorentzVector p4M=prodH->Get4Momentum(); 02525 G4int PDG=prodH->GetPDGCode(); 02526 G4int Chg=prodH->GetCharge(); 02527 G4int BaN=prodH->GetBaryonNumber(); 02528 G4cout<<"-EMC->>G4QIonIonCollis::Breed:String=Hadr,H#"<<aH<<" filled,4M=" 02529 <<p4M<<", PDG="<<PDG<<", Chg="<<Chg<<", BaN="<<BaN<<G4endl; 02530 #endif 02531 } 02532 } 02533 else 02534 { 02535 G4Quasmon* stringQuasmon = new G4Quasmon(miQC, curString4M);// String->Quas 02536 #ifdef debug 02537 G4cout<<"G4QIonIonCollision::Breeder:==> to Quasm="<<miQC<<curString4M 02538 <<", pNuc="<<theProjNucleus<<theProjNucleus.Get4Momentum()<<", tNuc=" 02539 <<theTargNucleus<<theTargNucleus.Get4Momentum()<<", NString=" 02540 <<strings.size()<<", nR="<<theResult->size()<<", nQ=" 02541 <<theQuasmons.size()<<G4endl; 02542 #endif 02543 theQuasmons.push_back(stringQuasmon); 02544 delete sHad; 02545 tmpQHadVec->clear(); 02546 delete tmpQHadVec; // WhoCallsDecayQHadron is responsible for clear&delete 02547 break; // Go out of the main StringDecay LOOP 02548 } 02549 tmpQHadVec->clear(); 02550 delete tmpQHadVec; // Who calls DecayQHadron is responsible for clear&delete 02551 break; // Go out of the main String Decay LOOP 02552 } 02553 } // End of the DecayOfTheLast 02554 } // End of IF(String-Hadron fusion) 02555 } // End of IF(NO_Hadrons) for String-String and String-Hadron fusion 02556 // The last hope is to CORREC the string, using other strings (ForwardInLOOP) 02557 #ifdef debug 02558 G4cout<<"G4QIonIonCollision::Breeder: theH="<<theHadrons<<"?=0, next="<<next<<G4endl; 02559 #endif 02560 if(!theHadrons && next < strings.size()) // ForwardInLOOP strings exist 02561 { 02562 // @@ string can be not convertable to one hadron (2,0.0,0,2,0) - To be improved 02563 G4QContent miQC=curString->GetQC(); // QContent of the Lightest Hadron 02564 G4int miPDG=miQC.GetSPDGCode();// PDG of the Lightest Hadron 02565 #ifdef debug 02566 G4cout<<">>>G4QIonIonCollision::Breeder: SQC="<<miQC<<", miSPDG="<<miPDG<<G4endl; 02567 #endif 02568 G4double miM=0.; // Prototype of the Mass of the Cur LightestHadron 02569 if(miPDG!=10) miM=G4QPDGCode(miPDG).GetMass(); // Mass of theCurLightestOneHadron 02570 else 02571 { 02572 G4QChipolino QCh(miQC); // define the TotalString as a Chipolino 02573 miM=QCh.GetQPDG1().GetMass()+QCh.GetQPDG2().GetMass();//MinMass of theStringChipo 02574 } 02575 G4double cM2=curString4M.m2(); // Actual squared mass of the Cur String 02576 #ifdef debug 02577 G4cout<<">>>G4QIonIonCollision::Breeder: minMass="<<miM<<", realM2="<<cM2<<G4endl; 02578 #endif 02579 G4double cM=0.; 02580 if(cM2>0.) 02581 { 02582 cM=std::sqrt(cM2); 02583 if(std::fabs(cM-miM) < eps) // Convert to hadron(2 hadrons) w/o calculations 02584 { 02585 if(miPDG!=10) 02586 { 02587 G4QHadron* sHad = new G4QHadron(miPDG,curString4M); 02588 theResult->push_back(sHad);// Fill the curString as a hadron 02589 #ifdef debug 02590 G4cout<<">>>G4QIonIonCollision::Breeder:S->H="<<miPDG<<curString4M<<G4endl; 02591 #endif 02592 } 02593 else 02594 { 02595 G4QChipolino QCh(miQC); // define TotalResidual as a Chipolino 02596 G4QPDGCode h1QPDG=QCh.GetQPDG1(); // QPDG of the first hadron 02597 G4double h1M =h1QPDG.GetMass(); // Mass of the first hadron 02598 G4QPDGCode h2QPDG=QCh.GetQPDG2(); // QPDG of the second hadron 02599 G4double h2M =h2QPDG.GetMass(); // Mass of the second hadron 02600 G4double pt1 =h1M/(h1M+h2M); // 4-mom part of the first hadron 02601 G4LorentzVector h14M=pt1*curString4M; // 4-mom of the first hadron 02602 G4LorentzVector h24M=curString4M-h14M;// 4-mom of the second hadron 02603 G4QHadron* h1H = new G4QHadron(h1QPDG.GetPDGCode(),h14M); 02604 theResult->push_back(h1H); // (delete equivalent) 02605 #ifdef debug 02606 G4LorentzVector f4M=h1H->Get4Momentum(); 02607 G4int fPD=h1H->GetPDGCode(); 02608 G4int fCg=h1H->GetCharge(); 02609 G4int fBN=h1H->GetBaryonNumber(); 02610 G4cout<<"-EMC->>G4QIonIonCollision::Breed:Str=2HadrAR Prod-F is filled, f4M=" 02611 <<f4M<<", fPDG="<<fPD<<", fCg="<<fCg<<", fBN="<<fBN<<G4endl; 02612 #endif 02613 G4QHadron* h2H = new G4QHadron(h2QPDG.GetPDGCode(),h24M); 02614 theResult->push_back(h2H); // (delete equivalent) 02615 #ifdef debug 02616 G4LorentzVector s4M=h2H->Get4Momentum(); 02617 G4int sPD=h2H->GetPDGCode(); 02618 G4int sCg=h2H->GetCharge(); 02619 G4int sBN=h2H->GetBaryonNumber(); 02620 G4cout<<"-EMC->>G4QIonIonCollision::Breed:Str=2HadrAR Prod-S is filled, s4M=" 02621 <<s4M<<", sPDG="<<sPD<<", sCg="<<sCg<<", sBN="<<sBN<<G4endl; 02622 #endif 02623 } 02624 continue; // Continue the LOOP over the curStrings 02625 } 02626 else // Try to recover (+/-) to min Mass 02627 { 02628 G4ThreeVector cP=curString4M.vect(); // Momentum of the curString 02629 G4double cE=curString4M.e(); // Energy of the curString 02630 G4ThreeVector curV=cP/cE; // curRelativeVelocity 02631 G4double miM2=miM*miM; 02632 G4int restr=0; // To use beyon the LOOP for printing 02633 G4int fustr=0; // Selected String-partner (0 = NotFound) 02634 G4double selX=0.; // Selected value of x 02635 G4double maD=-DBL_MAX; // Maximum Free Mass 02636 G4double Vmin=DBL_MAX; // min Velocity Distance 02637 G4LorentzVector s4M(0.,0.,0.,0.); // Selected 4-mom of the hadron 02638 #ifdef debug 02639 G4cout<<"G4QIonIonCollision::Breeder: TryRecover, cV="<<curV<<G4endl; 02640 #endif 02641 nOfStr=strings.size(); 02642 for(restr=next; restr < nOfStr; ++restr) if(restr != astring) 02643 { 02644 G4QString* pString=strings[restr]; 02645 G4LorentzVector p4M=pString->Get4Momentum(); 02646 G4ThreeVector pP=p4M.vect(); // Momentum of the partnerString 02647 G4double pE=p4M.e(); // Energy of the partnerString 02648 G4double D2=cE*pE-cP.dot(pP); 02649 G4double pM2=p4M.m2(); 02650 G4double dM4=pM2*(miM2-cM2); 02651 G4double D=D2*D2+dM4; 02652 G4double x=-1.; // Bad preexpectation 02653 if(D >= 0. && pM2>.01) x=(std::sqrt(D)-D2)/pM2; // what we should get from p 02654 #ifdef debug 02655 else G4cout<<"G4QIonIonCollis::Breed:D="<<D<<",D2="<<D2<<",dM="<<dM4<<G4endl; 02656 G4cout<<"G4QIonIonCollision::Breed: pM2="<<pM2<<",D2="<<D2<<",x="<<x<<G4endl; 02657 #endif 02658 if(x > 0. && x < 1.) // We are getting x part of p4M 02659 { 02660 G4QContent pQC=pString->GetQC(); // Quark Content of The Partner 02661 G4int pPDG=pQC.GetSPDGCode();// PDG of The Lightest Hadron for the Partner 02662 G4double pM=0.; // Mass of the LightestHadron 02663 if(pPDG==10) 02664 { 02665 G4QChipolino QCh(pQC); // define the TotalString as a Chipolino 02666 pM=QCh.GetQPDG1().GetMass()+QCh.GetQPDG2().GetMass();// Mass of Chipolino 02667 } 02668 else pM=G4QPDGCode(pPDG).GetMass();// Mass of theLightestHadron for Partner 02669 G4double rM=std::sqrt(pM2); // Real mass of the string-partner 02670 G4double delta=(1.-x)*rM-pM;// @@ Minimum CM disterbance measure 02671 if(delta > 0. && delta > maD) 02672 { 02673 maD=delta; 02674 #ifdef debug 02675 G4cout<<"G4QIonIonCollision::Breed: Subtr,S#"<<restr<<",d="<<maD<<G4endl; 02676 #endif 02677 fustr=restr; 02678 selX=x; 02679 s4M=p4M; 02680 } 02681 } 02682 else if(x <= 0.) // We are adding to p4M, so use RelVelocity 02683 { 02684 G4ThreeVector pV=pP/pE; // curRelativeVelocity 02685 G4double dV=(curV-pV).mag2();// SquaredDifferenceBetweenRelVelocities 02686 if(dV < Vmin) 02687 { 02688 #ifdef debug 02689 G4cout<<"G4QIonIonCollision::Breed: FreeAdd,S#"<<restr<<",x="<<x<<G4endl; 02690 #endif 02691 Vmin=dV; 02692 fustr=restr; 02693 selX=x; 02694 s4M=p4M; 02695 } 02696 } 02697 #ifdef debug 02698 G4cout<<"G4QIonIonCollision::Breed:EndOfLOOP r="<<restr<<"<"<<nOfStr<<G4endl; 02699 #endif 02700 } // End of the LOOP over string-partners for Correction 02701 #ifdef debug 02702 G4cout<<"G4QIonIonCollision::Breeder: AfterLOOP fustr="<<fustr<<G4endl; 02703 #endif 02704 if(fustr) 02705 { 02706 #ifdef edebug 02707 G4LorentzVector sum4M=s4M+curString4M; 02708 G4cout<<"G4QIonIonCollision::Breeder: Found Sum4M="<<sum4M<<G4endl; 02709 #endif 02710 G4QString* pString=strings[fustr]; 02711 curString4M+=selX*s4M; 02712 if(std::abs(miPDG)%10 > 2) // Decay String-Hadron-Resonance 02713 { 02714 G4Quasmon Quasm; 02715 G4QHadron* sHad = new G4QHadron(miPDG,curString4M); 02716 G4QHadronVector* tmpQHadVec=Quasm.DecayQHadron(sHad); // It deletes sHad 02717 #ifdef debug 02718 G4cout<<"G4QIonIonCollision::Breed:DecStH,nH="<<tmpQHadVec->size()<<G4endl; 02719 #endif 02720 for(unsigned aH=0; aH < tmpQHadVec->size(); aH++) 02721 { 02722 theResult->push_back((*tmpQHadVec)[aH]);//TheDecayProdOfHadron is filled 02723 #ifdef debug 02724 G4QHadron* prodH =(*tmpQHadVec)[aH]; 02725 G4LorentzVector p4M=prodH->Get4Momentum(); 02726 G4int PDG=prodH->GetPDGCode(); 02727 G4int Chg=prodH->GetCharge(); 02728 G4int BaN=prodH->GetBaryonNumber(); 02729 G4cout<<"-EMC->>G4QIonIonCollision::Breed:St=Had,pH#"<<aH<<" filled, 4M=" 02730 <<p4M<<", PDG="<<PDG<<", Chg="<<Chg<<", BaN="<<BaN<<G4endl; 02731 #endif 02732 } 02733 tmpQHadVec->clear(); 02734 delete tmpQHadVec; // Who calls DecayQHadron is responsibleRorClear&Delete 02735 } 02736 else if(miPDG == 10) // ==> Decay Hadron-Chipolino 02737 { 02738 G4QChipolino QCh(miQC); // define theTotalResid as aChipolino 02739 G4QPDGCode h1QPDG=QCh.GetQPDG1(); // QPDG of the first hadron 02740 G4double h1M =h1QPDG.GetMass();// Mass of the first hadron 02741 G4QPDGCode h2QPDG=QCh.GetQPDG2(); // QPDG of the second hadron 02742 G4double h2M =h2QPDG.GetMass();// Mass of the second hadron 02743 G4double ttM =curString4M.m(); // Real Mass of the Chipolino 02744 if(h1M+h2M<ttM+eps) // Two particles decay of Chipolino 02745 { 02746 G4LorentzVector h14M(0.,0.,0.,h1M); 02747 G4LorentzVector h24M(0.,0.,0.,h2M); 02748 if(std::fabs(ttM-h1M-h2M)<=eps) 02749 { 02750 G4double part1=h1M/(h1M+h2M); 02751 h14M=part1*curString4M; 02752 h24M=curString4M-h14M; 02753 } 02754 else 02755 { 02756 if(!G4QHadron(curString4M).DecayIn2(h14M,h24M)) 02757 { 02758 G4cerr<<"***G4QIonIonCollision::Breeder: tM="<<ttM<<"->h1="<<h1QPDG 02759 <<"("<<h1M<<")+h2="<<h1QPDG<<"("<<h2M<<")="<<h1M+h2M<<G4endl; 02760 G4Exception("G4QIonIonCollision::Breeder:","72",FatalException,"CD"); 02761 } 02762 } 02763 G4QHadron* h1H = new G4QHadron(h1QPDG.GetPDGCode(),h14M); 02764 theResult->push_back(h1H); // (delete equivalent) 02765 #ifdef debug 02766 G4LorentzVector f4M=h1H->Get4Momentum(); 02767 G4int fPD=h1H->GetPDGCode(); 02768 G4int fCg=h1H->GetCharge(); 02769 G4int fBN=h1H->GetBaryonNumber(); 02770 G4cout<<"-EMC->>>G4QIonIonCollision::Breed:Str=Hadr Prod-F's filled,f4M=" 02771 <<f4M<<", fPDG="<<fPD<<", fCg="<<fCg<<", fBN="<<fBN<<G4endl; 02772 #endif 02773 G4QHadron* h2H = new G4QHadron(h2QPDG.GetPDGCode(),h24M); 02774 theResult->push_back(h2H); // (delete equivalent) 02775 #ifdef debug 02776 G4LorentzVector s4M=h2H->Get4Momentum(); 02777 G4int sPD=h2H->GetPDGCode(); 02778 G4int sCg=h2H->GetCharge(); 02779 G4int sBN=h2H->GetBaryonNumber(); 02780 G4cout<<"-EMC->>>G4QIonIonCollision::Breed:Str=Hadr Prod-S's filled,s4M=" 02781 <<s4M<<", sPDG="<<sPD<<", sCg="<<sCg<<", sBN="<<sBN<<G4endl; 02782 #endif 02783 #ifdef edebug 02784 G4cout<<"-EMC-Chipo.G4QIonIonCollision::Breed:DecCHECK,c4M="<<curString4M 02785 <<", ChQC="<<miQC<<", d4M="<<curString4M-h14M-h24M<<G4endl; 02786 #endif 02787 } 02788 else 02789 { 02790 G4cerr<<"***G4QInel::Breeder: tM="<<ttM<<miQC<<"->h1="<<h1QPDG<<"(" <<h1M 02791 <<")+h2="<<h1QPDG<<"("<<h2M<<") = "<<h1M+h2M<<G4endl; 02792 G4Exception("G4QIonIonCollision::Breeder:","72",FatalException,"ChiDec"); 02793 } 02794 } 02795 else 02796 { 02797 G4QHadron* sHad = new G4QHadron(miPDG,curString4M); 02798 theResult->push_back(sHad); // The original string-hadron is filled 02799 #ifdef debug 02800 G4cout<<"-EMC->>>G4QIonIonCollision::Breeder:Str=Hadr Filled, 4M=" 02801 <<curString4M<<", PDG="<<miPDG<<G4endl; 02802 #endif 02803 } 02804 G4double corF=1-selX; 02805 G4QParton* Left=pString->GetLeftParton(); 02806 G4QParton* Right=pString->GetRightParton(); 02807 Left->Set4Momentum(corF*Left->Get4Momentum()); 02808 Right->Set4Momentum(corF*Right->Get4Momentum()); 02809 #ifdef edebug 02810 G4cout<<"-EMC-...Cor...G4QIonIonCollision::Breeder:CorCHECK Sum="<<sum4M 02811 <<" =? "<<curString4M+pString->Get4Momentum()<<", M="<<miM<<" ?= " 02812 <<curString4M.m()<<G4endl; 02813 #endif 02814 #ifdef debug 02815 G4cout<<">>>G4QIonIonCollision::Breeder:*Corrected* String->Hadr="<<miPDG 02816 <<curString4M<<" by String #"<<fustr<<G4endl; 02817 #endif 02818 continue; // Continue the LOOP over the curStrings 02819 } // End of Found combination for the String on string correction 02820 } // End of the Try-to-recover String+String correction algorithm 02821 } // End of IF(CM2>0.) 02822 } // End of IF(Can try to correct by String-String) 02823 #ifdef debug 02824 else G4cerr<<"***G4QIonIonCollision::Breed:**No SSCorrection**, next="<<next<<G4endl; 02825 #endif 02826 // ------------ At this point we can reduce the 3/-3 meson to 1/-1 meson ------------ 02827 G4QParton* lpcS=curString->GetLeftParton(); 02828 G4QParton* rpcS=curString->GetRightParton(); 02829 G4int lPDGcS=lpcS->GetPDGCode(); 02830 G4int rPDGcS=rpcS->GetPDGCode(); 02831 if (lPDGcS==3 && rPDGcS==-3) 02832 { 02833 lpcS->SetPDGCode( 1); 02834 rpcS->SetPDGCode(-1); 02835 } 02836 else if(lPDGcS==-3 && rPDGcS==3) 02837 { 02838 lpcS->SetPDGCode(-1); 02839 rpcS->SetPDGCode( 1); 02840 } 02841 // -------- Now the only hope is Correction, using the already prodused Hadrons ----- 02842 G4int nofRH=theResult->size(); // #of resulting Hadrons 02843 #ifdef debug 02844 G4cout<<"G4QIonIonCollision::Breeder: theH="<<theHadrons<<", #OfH="<<nofRH<<G4endl; 02845 #endif 02846 if(!theHadrons && nofRH) // Hadrons are existing for SH Correction 02847 { 02848 #ifdef debug 02849 G4cerr<<"!G4QIonIonCollision::Breeder:CanTrySHCor, nH="<<theResult->size()<<G4endl; 02850 #endif 02851 // @@ string can be not convertable to one hadron (2,0.0,0,2,0) - To be improved 02852 G4QContent miQC=curString->GetQC(); // QContent of the Lightest Hadron 02853 G4int miPDG=miQC.GetSPDGCode(); // PDG of the Lightest Hadron 02854 G4double miM=0.; // Prototype ofMass of theCurLightestHadron 02855 if(miPDG==10) // Mass of the Cur Lightest ChipolinoHadron 02856 { 02857 G4QChipolino QCh(miQC); // define the TotalString as a Chipolino 02858 miM=QCh.GetQPDG1().GetMass()+QCh.GetQPDG2().GetMass();//MinMass of theStringChipo 02859 } 02860 else miM=G4QPDGCode(miPDG).GetMass(); // Mass of the Cur Lightest Hadron 02861 G4double spM=0.; // Mass of the selected Hadron-Partner 02862 G4ThreeVector cP=curString4M.vect(); // Momentum of the curString 02863 G4double cE=curString4M.e(); // Energy of the curString 02864 G4ThreeVector curV=cP/cE; // curRelativeVelocity 02865 G4int reha=0; // Hadron # to use beyon the LOOP 02866 G4int fuha=0; // Selected Hadron-partner (0 = NotFound) 02867 G4double dMmin=DBL_MAX; // min Excess of the mass 02868 G4LorentzVector s4M(0.,0.,0.,0.); // Selected 4-mom of the Hadron+String 02869 G4double sM=0.; // Selected Mass of the Hadron+String 02870 for (reha=next; reha < nofRH; reha++) // LOOP over already collected hadrons 02871 { 02872 G4QHadron* pHadron=(*theResult)[reha];// Pointer to the current Partner-Hadron 02873 G4LorentzVector p4M=pHadron->Get4Momentum(); 02874 G4double pM=p4M.m(); // Mass of the Partner-Hadron 02875 G4LorentzVector t4M=p4M+curString4M; // Total momentum of the compound 02876 G4double tM2=t4M.m2(); // Squared total mass of the compound 02877 if(tM2 >= sqr(pM+miM+eps)) // Condition of possible correction 02878 { 02879 G4double tM=std::sqrt(tM2); // Mass of the Hadron+String compound 02880 G4double dM=tM-pM-miM; // Excess of the compound mass 02881 if(dM < dMmin) 02882 { 02883 dMmin=dM; 02884 fuha=reha; 02885 spM=pM; 02886 s4M=t4M; 02887 sM=tM; 02888 } 02889 } 02890 } // End of the LOOP over string-partners for Correction 02891 if(fuha) // The hadron-partner was found 02892 { 02893 G4QHadron* pHadron=(*theResult)[fuha];// Necessary for update 02894 G4LorentzVector mi4M(0.,0.,0.,miM); // Prototype of the new String=Hadron 02895 if(miM+spM<sM+eps) // Decay into CorrectedString+theSameHadron 02896 { 02897 G4LorentzVector sp4M(0.,0.,0.,spM); 02898 if(std::fabs(sM-miM-spM)<=eps) 02899 { 02900 G4double part1=miM/(miM+spM); 02901 mi4M=part1*s4M; 02902 sp4M=s4M-mi4M; 02903 } 02904 else 02905 { 02906 if(!G4QHadron(s4M).DecayIn2(mi4M,sp4M)) 02907 { 02908 G4cerr<<"***G4QIonIonCollision::Breeder: *SH*, tM="<<sM<<"->h1=("<<miPDG 02909 <<")"<<miM<<" + h2="<<spM<<" = "<<miM+spM<<G4endl; 02910 G4Exception("G4QIonIonCollision::Breeder:","72",FatalException,"SHChiDec"); 02911 } 02912 } 02913 pHadron->Set4Momentum(sp4M); 02914 #ifdef debug 02915 G4cout<<"-EMC->...>G4QIonIonCollision::Breed: H# "<<fuha<<" is updated, new4M=" 02916 <<sp4M<<G4endl; 02917 #endif 02918 } 02919 else 02920 { 02921 G4cerr<<"***G4QInel::Breeder: HS Failed, tM="<<sM<<"->h1M=("<<miPDG<<")"<<miM 02922 <<"+h2M="<<spM<<" = "<<miM+spM<<G4endl; 02923 G4Exception("G4QIonIonCollision::Breeder:","72",FatalException,"HSChipoliDec"); 02924 } 02925 if(std::abs(miPDG)%10 > 2) // Decay Hadron-Resonans 02926 { 02927 G4Quasmon Quasm; 02928 G4QHadron* sHad = new G4QHadron(miPDG,mi4M); 02929 G4QHadronVector* tmpQHadVec=Quasm.DecayQHadron(sHad); // It deletes sHad 02930 #ifdef debug 02931 G4cout<<"G4QInelast::Breeder: *HS* DecStrHad, nH="<<tmpQHadVec->size()<<G4endl; 02932 #endif 02933 for(unsigned aH=0; aH < tmpQHadVec->size(); aH++) 02934 { 02935 theResult->push_back((*tmpQHadVec)[aH]);// TheDecayProductOfTheHadronIsFilled 02936 #ifdef debug 02937 G4QHadron* prodH =(*tmpQHadVec)[aH]; 02938 G4LorentzVector p4M=prodH->Get4Momentum(); 02939 G4int PDG=prodH->GetPDGCode(); 02940 G4int Chg=prodH->GetCharge(); 02941 G4int BaN=prodH->GetBaryonNumber(); 02942 G4cout<<"-EMC->>>G4QIonIonCollision::Breed:Str+Hadr PrH#"<<aH<<" filled, 4M=" 02943 <<p4M<<", PDG="<<PDG<<", Chg="<<Chg<<", BaN="<<BaN<<G4endl; 02944 #endif 02945 } 02946 tmpQHadVec->clear(); 02947 delete tmpQHadVec; // Who calls DecayQHadron is responsible for clear & delete 02948 } 02949 else if(miPDG == 10) // ==> Decay Hadron-Chipolino 02950 { 02951 G4QChipolino QCh(miQC); // define the TotalResidual as a Chipolino 02952 G4QPDGCode h1QPDG=QCh.GetQPDG1(); // QPDG of the first hadron 02953 G4double h1M =h1QPDG.GetMass();// Mass of the first hadron 02954 G4QPDGCode h2QPDG=QCh.GetQPDG2(); // QPDG of the second hadron 02955 G4double h2M =h2QPDG.GetMass();// Mass of the second hadron 02956 G4double ttM =curString4M.m(); // Real Mass of the Chipolino 02957 if(h1M+h2M<miM+eps) // Two particles decay of Chipolino 02958 { 02959 G4LorentzVector h14M(0.,0.,0.,h1M); 02960 G4LorentzVector h24M(0.,0.,0.,h2M); 02961 if(std::fabs(ttM-h1M-h2M)<=eps) 02962 { 02963 G4double part1=h1M/(h1M+h2M); 02964 h14M=part1*mi4M; 02965 h24M=mi4M-h14M; 02966 } 02967 else 02968 { 02969 if(!G4QHadron(mi4M).DecayIn2(h14M,h24M)) 02970 { 02971 G4cerr<<"***G4QIonIonCollision::Breeder: HS tM="<<ttM<<"->h1="<<h1QPDG 02972 <<"("<<h1M<<")+h2="<<h1QPDG<<"("<<h2M<<")="<<h1M+h2M<<G4endl; 02973 G4Exception("G4QIonIonCollision::Breeder:","72",FatalException,"ChiDec"); 02974 } 02975 } 02976 G4QHadron* h1H = new G4QHadron(h1QPDG.GetPDGCode(),h14M); 02977 theResult->push_back(h1H); // (delete equivalent) 02978 #ifdef debug 02979 G4LorentzVector f4M=h1H->Get4Momentum(); 02980 G4int fPD=h1H->GetPDGCode(); 02981 G4int fCg=h1H->GetCharge(); 02982 G4int fBN=h1H->GetBaryonNumber(); 02983 G4cout<<"-EMC->>G4QIonIonCollision::Breed: CorStrHadr Prod-1 is filled, f4M=" 02984 <<f4M<<", fPDG="<<fPD<<", fCg="<<fCg<<", fBN="<<fBN<<G4endl; 02985 #endif 02986 G4QHadron* h2H = new G4QHadron(h2QPDG.GetPDGCode(),h24M); 02987 theResult->push_back(h2H); // (delete equivalent) 02988 #ifdef debug 02989 G4LorentzVector n4M=h2H->Get4Momentum(); 02990 G4int nPD=h2H->GetPDGCode(); 02991 G4int nCg=h2H->GetCharge(); 02992 G4int nBN=h2H->GetBaryonNumber(); 02993 G4cout<<"-EMC->>>G4QIonIonCollision::Breed:CorStrHadr Prod-2 is filled, n4M=" 02994 <<n4M<<", nPDG="<<nPD<<", nCg="<<nCg<<", nBN="<<nBN<<G4endl; 02995 #endif 02996 #ifdef edebug 02997 G4cout<<"-EMC-...HS-Chipo...G4QIonIonCollision::Breeder:DecCHECK, Ch4M=" 02998 <<mi4M<<", ChQC="<<miQC<<", d4M="<<mi4M-h14M-h24M<<G4endl; 02999 #endif 03000 } 03001 } 03002 else 03003 { 03004 G4QHadron* sHad = new G4QHadron(miPDG, mi4M); 03005 theResult->push_back(sHad); // The original string=hadron is filled 03006 #ifdef debug 03007 G4cout<<">>..>>G4QIonIonCollision::Breeder: CorStr=Hadr is Filled, 4M=" 03008 <<curString4M<<", StPDG="<<miPDG<<G4endl; 03009 #endif 03010 } 03011 #ifdef edebug 03012 G4cout<<"-EMC-...Cor...G4QIonIonCollision::Breeder:StHadCor CHECK Sum="<<s4M 03013 <<" =? "<<mi4M+pHadron->Get4Momentum()<<G4endl; 03014 #endif 03015 #ifdef debug 03016 G4cout<<">>>G4QIonIonCollision::Breeder:*Corrected* String+Hadr="<<miPDG 03017 <<mi4M<<" by Hadron #"<<reha<<G4endl; 03018 #endif 03019 continue; // Continue the LOOP over the curStrings 03020 } 03021 else 03022 { 03023 #ifdef debug 03024 G4cout<<"-EMC->>>G4QIonIonCollision::Breeder: Str+Hadr Failed, 4M="<<curString4M 03025 <<", PDG="<<miPDG<<G4endl; 03026 #endif 03027 } 03028 // @@@ convert string to Quasmon with curString4M 03029 G4QContent curStringQC=curString->GetQC(); 03030 G4Quasmon* stringQuasmon = new G4Quasmon(curStringQC, curString4M); 03031 theQuasmons.push_back(stringQuasmon); 03032 continue; // Continue the LOOP over the curStrings 03033 } // End of IF(Can try the String-Hadron correction 03034 } // End of IF(NO_Hadrons) = Problem solution namespace 03035 G4Quasmon tmpQ; // @@ an issue of Q to decay resonances 03036 G4int nHfin=0; 03037 if(theHadrons) nHfin=theHadrons->size(); 03038 else // !! Sum Up all strings and convert them in a Quasmon (Exception for development) 03039 { 03040 G4LorentzVector ss4M(0.,0.,0.,0.); 03041 G4QContent ssQC(0,0,0,0,0,0); 03042 G4int tnSt=strings.size(); 03043 for(G4int i=astring; i < tnSt; ++i) 03044 { 03045 G4LorentzVector pS4M=strings[i]->Get4Momentum(); // String 4-momentum 03046 ss4M+=pS4M; 03047 G4QContent pSQC=strings[i]->GetQC(); // String Quark Content 03048 ssQC+=pSQC; 03049 #ifdef debug 03050 G4cout<<"=--=>G4QIonIonCollision::Breed:S#"<<i<<",4M="<<pS4M<<",QC="<<pSQC<<G4endl; 03051 #endif 03052 } 03053 #ifdef debug 03054 G4cout<<"==>G4QIonIonCollision::Breed:AllStrings are summed up in a Quasmon"<<G4endl; 03055 #endif 03056 G4Quasmon* stringQuasmon = new G4Quasmon(ssQC, ss4M); 03057 theQuasmons.push_back(stringQuasmon); 03058 break; // break the LOOP over Strings 03059 } 03060 #ifdef debug 03061 G4cout<<"G4QIonIonCollision::Breeder: Trying to decay hadrons #ofHRes="<<nHfin<<G4endl; 03062 #endif 03063 for(G4int aTrack=0; aTrack<nHfin; aTrack++) 03064 { 03065 G4QHadron* curHadron=(*theHadrons)[aTrack]; 03066 G4int hPDG=curHadron->GetPDGCode(); 03067 #ifdef edebug 03068 G4LorentzVector curH4M=curHadron->Get4Momentum(); 03069 G4int curHCh=curHadron->GetCharge(); 03070 G4int curHBN=curHadron->GetBaryonNumber(); 03071 #endif 03072 #ifdef debug 03073 G4cout<<">>...>>G4QIonIonCollision::Breeder:S#"<<astring<<",H#"<<aTrack<<",PDG=" 03074 <<hPDG<<",4M="<<curHadron->Get4Momentum()<<G4endl; 03075 #endif 03076 if(std::abs(hPDG)%10 > 2) 03077 { 03078 G4QHadronVector* tmpQHadVec=tmpQ.DecayQHadron(curHadron); // It deletes curHadron 03079 #ifdef debug 03080 G4cout<<"G4QIonIonCollision::Breed:-DECAY'S DONE-,nH="<<tmpQHadVec->size()<<G4endl; 03081 #endif 03082 //G4int tmpS=tmpQHadVec->size(); // "The elegant method" (tested) is commented 03083 //theResult->resize(tmpS+theResult->size()); // Resize theQHadrons length 03084 //copy(tmpQHadVec->begin(), tmpQHadVec->end(), theResult->end()-tmpS); 03085 for(unsigned aH=0; aH < tmpQHadVec->size(); aH++) 03086 { 03087 theResult->push_back((*tmpQHadVec)[aH]);// TheDecayProduct of TheHadron is filled 03088 #ifdef edebug 03089 G4QHadron* prodH =(*tmpQHadVec)[aH]; 03090 G4LorentzVector p4M=prodH->Get4Momentum(); 03091 G4int PDG=prodH->GetPDGCode(); 03092 G4int Chg=prodH->GetCharge(); 03093 G4int BaN=prodH->GetBaryonNumber(); 03094 curH4M-=p4M; 03095 curString4M-=p4M; 03096 curStrChg-=Chg; 03097 curStrBaN-=BaN; 03098 curHCh-=Chg; 03099 curHBN-=BaN; 03100 G4cout<<"-EMC->.>>G4QIonIonCollision::Breed:Str*Filled, 4M="<<p4M<<", PDG="<<PDG 03101 <<", Chg="<<Chg<<", BaN="<<BaN<<G4endl; 03102 #endif 03103 } 03104 #ifdef edebug 03105 G4cout<<"-EMC-.G4QIn::Br:Dec,r4M="<<curH4M<<",rC="<<curHCh<<",rB="<<curHBN<<G4endl; 03106 #endif 03107 tmpQHadVec->clear(); 03108 delete tmpQHadVec; // Who calls DecayQHadron is responsible for clear & delete 03109 } 03110 else // Chipolino is not checked here 03111 { 03112 theResult->push_back(curHadron); // The original hadron is filled 03113 #ifdef edebug 03114 curString4M-=curH4M; 03115 G4int curCh=curHadron->GetCharge(); 03116 G4int curBN=curHadron->GetBaryonNumber(); 03117 curStrChg-=curCh; 03118 curStrBaN-=curBN; 03119 G4cout<<"-EMC->>..>>G4QIonIonCollision::Breeder: curH filled 4M="<<curH4M<<",PDG=" 03120 <<curHadron->GetPDGCode()<<", Chg="<<curCh<<", BaN="<<curBN<<G4endl; 03121 #endif 03122 } 03123 } 03124 // clean up (the issues are filled to theResult) 03125 if(theHadrons) delete theHadrons; 03126 #ifdef edebug 03127 G4cout<<"-EMC-.......G4QIonIonCollision::Breeder: StringDecay CHECK, r4M="<<curString4M 03128 <<", rChg="<<curStrChg<<", rBaN="<<curStrBaN<<G4endl; 03129 #endif 03130 } // End of the main LOOP over decaying strings 03131 G4LorentzVector rp4M=theProjNucleus.Get4Momentum(); // projNucleus 4-momentum in LS 03132 G4int rpPDG=theProjNucleus.GetPDG(); 03133 G4QHadron* resPNuc = new G4QHadron(rpPDG,rp4M); 03134 theResult->push_back(resPNuc); // Fill the residual nucleus 03135 G4LorentzVector rt4M=theTargNucleus.Get4Momentum(); // Nucleus 4-momentum in LS 03136 G4int rtPDG=theTargNucleus.GetPDG(); 03137 G4QHadron* resTNuc = new G4QHadron(rtPDG,rt4M); 03138 theResult->push_back(resTNuc); // Fill the residual nucleus 03139 #ifdef edebug 03140 G4LorentzVector s4M(0.,0.,0.,0.); // Sum of the Result in LS 03141 G4int rCh=totChg; 03142 G4int rBN=totBaN; 03143 G4int nHadr=theResult->size(); 03144 G4int nQuasm=theQuasmons.size(); 03145 G4cout<<"-EMCLS-G4QInel::Breeder: #ofHadr="<<nHadr<<", #OfQ="<<nQuasm<<", rPA="<<rp4M.m() 03146 <<"="<<G4QNucleus(rpPDG).GetGSMass()<<", rTA="<<rt4M.m()<<"=" 03147 <<G4QNucleus(rtPDG).GetGSMass()<<G4endl; 03148 for(G4int i=0; i<nHadr; i++) 03149 { 03150 G4LorentzVector hI4M=(*theResult)[i]->Get4Momentum(); 03151 s4M+=hI4M; 03152 G4int hChg=(*theResult)[i]->GetCharge(); 03153 rCh-=hChg; 03154 G4int hBaN=(*theResult)[i]->GetBaryonNumber(); 03155 rBN-=hBaN; 03156 G4cout<<"-EMCLS-G4QIonIonCollision::Breeder: Hadron#"<<i<<", 4M="<<hI4M<<", PDG=" 03157 <<(*theResult)[i]->GetPDGCode()<<", C="<<hChg<<", B="<<hBaN<<G4endl; 03158 } 03159 for(G4int i=0; i<nQuasm; i++) 03160 { 03161 G4LorentzVector hI4M=theQuasmons[i]->Get4Momentum(); 03162 s4M+=hI4M; 03163 G4int hChg=theQuasmons[i]->GetCharge(); 03164 rCh-=hChg; 03165 G4int hBaN=theQuasmons[i]->GetBaryonNumber(); 03166 rBN-=hBaN; 03167 G4cout<<"-EMCLS-G4QIonIonCollision::Breeder: Quasmon#"<<i<<", 4M="<<hI4M<<", C="<<hChg 03168 <<", B="<<hBaN<<G4endl; 03169 } 03170 G4cout<<"-EMCLS-G4QInel::Breed: LS r4M="<<s4M-totLS4M<<",rC="<<rCh<<",rB="<<rBN<<G4endl; 03171 #endif 03172 // Now we need to coolect particles for creation of a Quasmon @@ improve !! 03173 G4int nRes=theResult->size(); 03174 if(nRes>2) 03175 { 03176 G4QHadronVector::iterator ih; 03177 G4QHadronVector::iterator nih; 03178 G4QHadronVector::iterator mih; 03179 G4QHadronVector::iterator lst=theResult->end(); 03180 lst--; 03181 G4double minMesEn=DBL_MAX; 03182 G4double minBarEn=DBL_MAX; 03183 G4bool nfound=false; 03184 G4bool mfound=false; 03185 for(ih = theResult->begin(); ih < theResult->end(); ++ih) if(ih != lst) 03186 { 03187 G4LorentzVector h4M=(*ih)->Get4Momentum(); 03188 G4int hPDG=(*ih)->GetPDGCode(); 03189 #ifdef debug 03190 G4cout<<"%->G4QIonIonCollision::Breeder: TRY hPDG="<<hPDG<<", h4M="<<h4M<<G4endl; 03191 #endif 03192 if(hPDG>1111 && hPDG<3333) 03193 { 03194 G4double bE=h4M.e()-(*ih)->GetMass(); 03195 if(bE < minBarEn) 03196 { 03197 minBarEn=bE; 03198 nih=ih; 03199 nfound=true; 03200 } 03201 } 03202 else if(hPDG>-1111) 03203 { 03204 G4double mE=h4M.e(); 03205 if(mE < minMesEn) 03206 { 03207 minMesEn=mE; 03208 mih=ih; 03209 mfound=true; 03210 } 03211 } 03212 } 03213 if(nfound && mfound && minMesEn+minBarEn < maxEn) 03214 { 03215 G4QHadron* resNuc = (*theResult)[nRes-1]; // ResidualNucleus is theLastH 03216 theResult->pop_back(); // Fill the QHad-nucleus later 03217 G4LorentzVector q4M=(*nih)->Get4Momentum()+(*mih)->Get4Momentum(); 03218 G4QContent qQC=(*nih)->GetQC()+(*mih)->GetQC(); 03219 maxEn -= minBarEn+minMesEn; // Reduce the energy limit 03220 #ifdef debug 03221 G4cout<<"%->G4QIonIonCollision::Breeder:Exclude,4M="<<(*nih)->Get4Momentum() 03222 <<", & 4m="<<(*mih)->Get4Momentum()<<G4endl; 03223 #endif 03224 delete (*nih); 03225 delete (*mih); 03226 if(nih > mih) // Exclude nucleon first 03227 { 03228 theResult->erase(nih); // erase Baryon from theResult 03229 theResult->erase(mih); // erase Meson from theResult 03230 } 03231 else // Exclude meson first 03232 { 03233 theResult->erase(mih); // erase Baryon from theResult 03234 theResult->erase(nih); // erase Meson from theResult 03235 } 03236 #ifdef debug 03237 G4cout<<"%->G4QI::Breed: BeforeLOOP, dE="<<maxEn<<", nR="<<theResult->size()<<G4endl; 03238 #endif 03239 if(maxEn > 0.) // More hadrons to absorbe 03240 { 03241 for(ih = theResult->begin(); ih < theResult->end(); ih++) 03242 { 03243 G4LorentzVector h4M=(*ih)->Get4Momentum(); 03244 G4int hPDG=(*ih)->GetPDGCode(); 03245 G4double dE=0.; 03246 if (hPDG> 1111 && hPDG<3333) dE=h4M.e()-(*ih)->GetMass(); // Baryons 03247 else if(hPDG>-1111) dE=h4M.e(); // Mesons 03248 else continue; // Don't consider anti-baryons 03249 if(dE < maxEn) 03250 { 03251 maxEn-=dE; 03252 q4M+=h4M; 03253 qQC+=(*ih)->GetQC(); 03254 #ifdef debug 03255 G4cout<<"%->G4QIonIonCollision::Breed:Exclude,4M="<<h4M<<",dE="<<maxEn<<G4endl; 03256 #endif 03257 delete (*ih); 03258 theResult->erase(ih); // erase Hadron from theResult 03259 --ih; 03260 } 03261 } 03262 } 03263 G4Quasmon* softQuasmon = new G4Quasmon(qQC, q4M); // SoftPart -> Quasmon 03264 #ifdef debug 03265 G4cout<<"%->G4QIonIonCollision::Breed:QuasmonIsFilled,4M="<<q4M<<",QC="<<qQC<<G4endl; 03266 #endif 03267 theQuasmons.push_back(softQuasmon); 03268 theResult->push_back(resNuc); 03269 } 03270 #ifdef edebug 03271 G4LorentzVector f4M(0.,0.,0.,0.); // Sum of the Result in LS 03272 G4int fCh=totChg; 03273 G4int fBN=totBaN; 03274 G4int nHd=theResult->size(); 03275 G4int nQm=theQuasmons.size(); 03276 G4cout<<"-EMCLS-G4QIonIonCollision::Breeder:#ofHadr="<<nHd<<", #OfQ="<<nQm 03277 <<",rPA="<<rt4M.m()<<"="<<G4QNucleus(rtPDG).GetGSMass() 03278 <<",rTA="<<rt4M.m()<<"="<<G4QNucleus(rtPDG).GetGSMass()<<G4endl; 03279 for(G4int i=0; i<nHd; i++) 03280 { 03281 G4LorentzVector hI4M=(*theResult)[i]->Get4Momentum(); 03282 f4M+=hI4M; 03283 G4int hChg=(*theResult)[i]->GetCharge(); 03284 fCh-=hChg; 03285 G4int hBaN=(*theResult)[i]->GetBaryonNumber(); 03286 fBN-=hBaN; 03287 G4cout<<"-EMCLS-G4QIonIonCollision::Breeder: Hadron#"<<i<<", 4M="<<hI4M<<", PDG=" 03288 <<(*theResult)[i]->GetPDGCode()<<", C="<<hChg<<", B="<<hBaN<<G4endl; 03289 } 03290 for(G4int i=0; i<nQm; i++) 03291 { 03292 G4LorentzVector hI4M=theQuasmons[i]->Get4Momentum(); 03293 f4M+=hI4M; 03294 G4int hChg=theQuasmons[i]->GetCharge(); 03295 fCh-=hChg; 03296 G4int hBaN=theQuasmons[i]->GetBaryonNumber(); 03297 fBN-=hBaN; 03298 G4cout<<"-EMCLS-G4QIonIonCollision::Breeder: Quasmon#"<<i<<", 4M="<<hI4M<<", C=" 03299 <<hChg<<", B="<<hBaN<<G4endl; 03300 } 03301 G4cout<<"-EMCLS-G4QInel::Breed:, r4M="<<f4M-totLS4M<<", rC="<<fCh<<",rB="<<fBN<<G4endl; 03302 #endif 03303 } // End of the soft Quasmon Creation 03304 return; 03305 } // End of Breeder
Definition at line 3551 of file G4QIonIonCollision.cc.
References FatalException, G4cerr, G4cout, G4endl, G4Exception(), G4UniformRand, and sqr().
Referenced by ExciteDiffParticipants(), and ExciteSingDiffParticipants().
03552 { 03553 // choose an x between Xmin and Xmax with P(x) ~ 1/x @@ M.K. -> 1/sqrt(x) 03554 //G4double range=Xmax-Xmin; 03555 if(Xmax == Xmin) return Xmin; 03556 if( Xmin < 0. || Xmax < Xmin) 03557 { 03558 G4cerr<<"***G4QIonIonCollision::ChooseX: Xmin="<<Xmin<<", Xmax="<<Xmax<< G4endl; 03559 G4Exception("G4QIonIonCollision::ChooseX:","72",FatalException,"Bad X or X-Range"); 03560 } 03561 //G4double x; 03562 //do {x=Xmin+G4UniformRand()*range;} while ( Xmin/x < G4UniformRand() ); 03563 G4double sxi=std::sqrt(Xmin); 03564 G4double x=sqr(sxi+G4UniformRand()*(std::sqrt(Xmax)-sxi)); 03565 #ifdef debug 03566 G4cout<<"G4QIonIonCollision::ChooseX: DiffractiveX="<<x<<G4endl; 03567 #endif 03568 return x; 03569 } // End of ChooseX
G4bool G4QIonIonCollision::ExciteDiffParticipants | ( | G4QHadron * | aPartner, | |
G4QHadron * | bPartner | |||
) | const [protected] |
Definition at line 3308 of file G4QIonIonCollision.cc.
References ChooseX(), G4cout, G4endl, GaussianPt(), G4QHadron::Get4Momentum(), G4QHadron::GetMass(), G4QHadron::Set4Momentum(), and sqr().
Referenced by G4QIonIonCollision().
03310 { 03311 G4LorentzVector Pprojectile=projectile->Get4Momentum(); 03312 G4double Mprojectile=projectile->GetMass(); 03313 G4double Mprojectile2=Mprojectile*Mprojectile; 03314 G4LorentzVector Ptarget=target->Get4Momentum(); 03315 G4double Mtarget=target->GetMass(); 03316 G4double Mtarget2=Mtarget*Mtarget; 03317 #ifdef debug 03318 G4cout<<"G4QInel::ExciteDiffPartici: Ep="<<Pprojectile.e()<<", Et="<<Ptarget.e()<<G4endl; 03319 #endif 03320 // Transform momenta to cms and then rotate parallel to z axis; 03321 G4LorentzVector Psum=Pprojectile+Ptarget; 03322 G4LorentzRotation toCms(-Psum.boostVector()); // Boost Rotation to CMS 03323 G4LorentzVector Ptmp=toCms*Pprojectile; 03324 if(Ptmp.pz()<=0.) // "String" moving backwards in CMS, abort collision !! ? M.K. 03325 { 03326 #ifdef debug 03327 G4cout<<"G4QIonIonCollision::ExciteDiffParticipants:*1* abort Collision!! *1*"<<G4endl; 03328 #endif 03329 return false; 03330 } 03331 toCms.rotateZ(-Ptmp.phi()); 03332 toCms.rotateY(-Ptmp.theta()); 03333 #ifdef debug 03334 G4cout<<"G4QIonIonCollision::ExciteDiffParticipantts:BeforBoost Pproj="<<Pprojectile 03335 <<",Ptarg="<<Ptarget<<G4endl; 03336 #endif 03337 G4LorentzRotation toLab(toCms.inverse()); // Boost Rotation to LabSys (LS) 03338 Pprojectile.transform(toCms); 03339 Ptarget.transform(toCms); 03340 #ifdef debug 03341 G4cout<< "G4QInelast::ExciteDiffParticipantts: AfterBoost Pproj="<<Pprojectile<<",Ptarg=" 03342 <<Ptarget<<", cms4M="<<Pprojectile+Ptarget<<G4endl; 03343 G4cout<<"G4QIonIonCollis::ExciteDiffParticipants:ProjX+="<<Pprojectile.plus()<<",ProjX-=" 03344 <<Pprojectile.minus()<<", tX+="<< Ptarget.plus()<<",tX-="<<Ptarget.minus()<<G4endl; 03345 #endif 03346 G4LorentzVector Qmomentum(0.,0.,0.,0.); 03347 G4int whilecount=0; 03348 #ifdef debug 03349 G4cout<<"G4QIonIonCollision::ExciteDiffParticipants: Before DO"<<G4endl; 03350 #endif 03351 do 03352 { 03353 // Generate pt 03354 G4double maxPtSquare=sqr(Ptarget.pz()); 03355 #ifdef debug 03356 G4cout<<"G4QIonIonCollision::ExciteDiffParticipants: maxPtSq="<<maxPtSquare<<G4endl; 03357 if(whilecount++>=500 && whilecount%100==0) // @@ M.K. Hardwired limits 03358 G4cout<<"G4QIonIonCollision::ExciteDiffParticipantts: can loop, loopCount="<<whilecount 03359 <<", maxPtSquare="<<maxPtSquare<<G4endl; 03360 #endif 03361 if(whilecount>1000) // @@ M.K. Hardwired limits 03362 { 03363 #ifdef debug 03364 G4cout<<"G4QIonIonCollision::ExciteDiffParticipants: *2* abort Loop!! *2*"<<G4endl; 03365 #endif 03366 return false; // Ignore this interaction 03367 } 03368 Qmomentum=G4LorentzVector(GaussianPt(widthOfPtSquare,maxPtSquare),0); 03369 #ifdef debug 03370 G4cout<<"G4QIonIonCollision::ExciteDiffParticipants: generated Pt="<<Qmomentum 03371 <<", ProjPt="<<Pprojectile+Qmomentum<<", TargPt="<<Ptarget-Qmomentum<<G4endl; 03372 #endif 03373 // Momentum transfer 03374 G4double Xmin=0.; 03375 G4double Xmax=1.; 03376 G4double Xplus =ChooseX(Xmin,Xmax); 03377 G4double Xminus=ChooseX(Xmin,Xmax); 03378 #ifdef debug 03379 G4cout<<"G4QIonIonCollision::ExciteDiffParticipant:X+="<<Xplus<<",X-="<<Xminus<<G4endl; 03380 #endif 03381 G4double pt2=Qmomentum.vect().mag2(); 03382 G4double Qplus =-pt2/Xminus/Ptarget.minus(); 03383 G4double Qminus= pt2/Xplus /Pprojectile.plus(); 03384 Qmomentum.setPz((Qplus-Qminus)/2); 03385 Qmomentum.setE( (Qplus+Qminus)/2); 03386 #ifdef debug 03387 G4cout<<"G4QInelast::ExciteDiffParticip: Qplus="<<Qplus<<", Qminus="<<Qminus<<", pt2=" 03388 <<pt2<<", Qmomentum="<<Qmomentum<<", ProjM="<<(Pprojectile+Qmomentum).mag() 03389 <<", TargM="<<(Ptarget-Qmomentum).mag()<<G4endl; 03390 #endif 03391 } while((Pprojectile+Qmomentum).mag2()<=Mprojectile2 || 03392 (Ptarget-Qmomentum).mag2()<=Mtarget2); 03393 Pprojectile += Qmomentum; 03394 Ptarget -= Qmomentum; 03395 #ifdef debug 03396 G4cout<<"G4QInelast::ExciteDiffParticipan: Proj(Q)="<<Pprojectile<<", Targ(Q)="<<Ptarget 03397 <<", Proj(back)="<<toLab*Pprojectile<<", Targ(bac)="<< toLab*Ptarget << G4endl; 03398 #endif 03399 // Transform back and update SplitableHadron Participant. 03400 Pprojectile.transform(toLab); 03401 Ptarget.transform(toLab); 03402 #ifdef debug 03403 G4cout<< "G4QIonIonCollision::ExciteDiffParticipants:TargetMass="<<Ptarget.mag()<<G4endl; 03404 #endif 03405 target->Set4Momentum(Ptarget); 03406 #ifdef debug 03407 G4cout<<"G4QIonIonCollision::ExciteDiffParticipant:ProjMass="<<Pprojectile.mag()<<G4endl; 03408 #endif 03409 projectile->Set4Momentum(Pprojectile); 03410 return true; 03411 } // End of ExciteDiffParticipants
G4bool G4QIonIonCollision::ExciteSingDiffParticipants | ( | G4QHadron * | aPartner, | |
G4QHadron * | bPartner | |||
) | const [protected] |
Definition at line 3415 of file G4QIonIonCollision.cc.
References ChooseX(), G4cout, G4endl, G4UniformRand, GaussianPt(), G4QHadron::Get4Momentum(), G4QHadron::GetMass(), G4QHadron::GetMass2(), G4QHadron::Set4Momentum(), and sqr().
Referenced by G4QIonIonCollision().
03417 { 03418 G4LorentzVector Pprojectile=projectile->Get4Momentum(); 03419 G4double Mprojectile=projectile->GetMass(); 03420 G4double Mprojectile2=Mprojectile*Mprojectile; 03421 G4LorentzVector Ptarget=target->Get4Momentum(); 03422 G4double Mtarget=target->GetMass(); 03423 G4double Mtarget2=Mtarget*Mtarget; 03424 #ifdef debug 03425 G4cout<<"G4QInel::ExSingDiffPartici: Ep="<<Pprojectile.e()<<", Et="<<Ptarget.e()<<G4endl; 03426 #endif 03427 G4bool KeepProjectile= G4UniformRand() > 0.5; 03428 // Reset minMass of the non diffractive particle to its value, (minus for rounding...) 03429 if(KeepProjectile ) 03430 { 03431 #ifdef debug 03432 G4cout<<"-1/2-G4QIonIonCollision::ExSingDiffParticipants: Projectile is fixed"<<G4endl; 03433 #endif 03434 Mprojectile2 = projectile->GetMass2()*(1.-perCent); // Isn't it too big reduction? M.K. 03435 } 03436 else 03437 { 03438 #ifdef debug 03439 G4cout<<"---1/2---G4QIonIonCollision::ExSingDiffParticipants: Target is fixed"<<G4endl; 03440 #endif 03441 Mtarget2 = target->GetMass2()*(1.-perCent); // @@ Isn't it too big reduction? M.K. 03442 } 03443 // @@ From this point it repeats the Diffractional excitation (? Use flag ?) 03444 // Transform momenta to cms and then rotate parallel to z axis; 03445 G4LorentzVector Psum=Pprojectile+Ptarget; 03446 G4LorentzRotation toCms(-Psum.boostVector()); // Boost Rotation to CMS 03447 G4LorentzVector Ptmp=toCms*Pprojectile; 03448 if(Ptmp.pz()<=0.) // "String" moving backwards in CMS, abort collision !! ? M.K. 03449 { 03450 #ifdef debug 03451 G4cout<<"G4QIonIonCollision::ExciteSingDiffParticipants: *1*abortCollision*1*"<<G4endl; 03452 #endif 03453 return false; 03454 } 03455 toCms.rotateZ(-Ptmp.phi()); 03456 toCms.rotateY(-Ptmp.theta()); 03457 #ifdef debug 03458 G4cout<<"G4QInel::ExciteSingDiffParticipantts: Be4Boost Pproj="<<Pprojectile<<", Ptarg=" 03459 <<Ptarget<<G4endl; 03460 #endif 03461 G4LorentzRotation toLab(toCms.inverse()); // Boost Rotation to LabSys (LS) 03462 Pprojectile.transform(toCms); 03463 Ptarget.transform(toCms); 03464 #ifdef debug 03465 G4cout<< "G4QInelast::ExciteDiffParticipantts: AfterBoost Pproj="<<Pprojectile<<"Ptarg=" 03466 <<Ptarget<<", cms4M="<<Pprojectile+Ptarget<<G4endl; 03467 03468 G4cout<<"G4QInelast::ExciteDiffParticipantts: ProjX+="<<Pprojectile.plus()<<", ProjX-=" 03469 <<Pprojectile.minus()<<", TargX+="<< Ptarget.plus()<<", TargX-="<<Ptarget.minus() 03470 <<G4endl; 03471 #endif 03472 G4LorentzVector Qmomentum(0.,0.,0.,0.); 03473 G4int whilecount=0; 03474 do 03475 { 03476 // Generate pt 03477 G4double maxPtSquare=sqr(Ptarget.pz()); 03478 if(whilecount++>=500 && whilecount%100==0) // @@ M.K. Hardwired limits 03479 #ifdef debug 03480 G4cout<<"G4QIonIonCollision::ExciteSingDiffParticipantts: can loop, loopCount=" 03481 <<whilecount<<", maxPtSquare="<<maxPtSquare<<G4endl; 03482 #endif 03483 if(whilecount>1000) // @@ M.K. Hardwired limits 03484 { 03485 #ifdef debug 03486 G4cout<<"G4QIonIonCollision::ExciteSingDiffParticipants:*2* abortLoop!! *2*"<<G4endl; 03487 #endif 03488 return false; // Ignore this interaction 03489 } 03490 Qmomentum=G4LorentzVector(GaussianPt(widthOfPtSquare,maxPtSquare),0); 03491 #ifdef debug 03492 G4cout<<"G4QInelast::ExciteSingDiffParticipants: generated Pt="<<Qmomentum 03493 <<", ProjPt="<<Pprojectile+Qmomentum<<", TargPt="<<Ptarget-Qmomentum<<G4endl; 03494 #endif 03495 // Momentum transfer 03496 G4double Xmin=0.; 03497 G4double Xmax=1.; 03498 G4double Xplus =ChooseX(Xmin,Xmax); 03499 G4double Xminus=ChooseX(Xmin,Xmax); 03500 #ifdef debug 03501 G4cout<<"G4QInel::ExciteSingDiffPartici: X-plus="<<Xplus<<", X-minus="<<Xminus<<G4endl; 03502 #endif 03503 G4double pt2=G4ThreeVector(Qmomentum.vect()).mag2(); 03504 G4double Qplus =-pt2/Xminus/Ptarget.minus(); 03505 G4double Qminus= pt2/Xplus /Pprojectile.plus(); 03506 if (KeepProjectile) 03507 Qminus=(projectile->GetMass2()+pt2)/(Pprojectile.plus()+Qplus) - Pprojectile.minus(); 03508 else Qplus=Ptarget.plus() - (target->GetMass2()+pt2)/(Ptarget.minus()-Qminus); 03509 Qmomentum.setPz((Qplus-Qminus)/2); 03510 Qmomentum.setE( (Qplus+Qminus)/2); 03511 #ifdef debug 03512 G4cout<<"G4QInel::ExciteDiffParticip: Qplus="<<Qplus<<", Qminus="<<Qminus<<", pt2=" 03513 <<pt2<<", Qmomentum="<<Qmomentum<<", ProjM="<<(Pprojectile+Qmomentum).mag() 03514 <<", TargM="<<(Ptarget-Qmomentum).mag()<<G4endl; 03515 #endif 03516 // while is different from the Double Diffractive Excitation (@@ !) 03517 //} while((Pprojectile+Qmomentum).mag2()<= Mprojectile2 || 03518 // (Ptarget-Qmomentum).mag2()<=Mtarget2); 03519 } while((Ptarget-Qmomentum).mag2()<=Mtarget2 || 03520 (Pprojectile+Qmomentum).mag2()<=Mprojectile2 || 03521 (Ptarget-Qmomentum).e() < 0. || (Pprojectile+Qmomentum).e() < 0.); 03522 Pprojectile += Qmomentum; 03523 Ptarget -= Qmomentum; 03524 #ifdef debug 03525 G4cout<<"G4QIonIonCollision::ExciteSingDiffParticipan: Proj(Q)="<<Pprojectile<<"(E=" 03526 <<Pprojectile.e()<<"), Targ(Q)="<<Ptarget<<"(E="<<Ptarget.e() 03527 <<"), Proj(back)="<<toLab*Pprojectile<<", Targ(bac)="<< toLab*Ptarget << G4endl; 03528 #endif 03529 // Transform back and update SplitableHadron Participant. 03530 Pprojectile.transform(toLab); 03531 Ptarget.transform(toLab); 03532 #ifdef debug 03533 G4cout<< "G4QIonIonCollision::ExciteSingleDiffParticipants: TgM="<<Ptarget.mag()<<G4endl; 03534 #endif 03535 target->Set4Momentum(Ptarget); 03536 #ifdef debug 03537 G4cout<<"G4QIonIonCollision::ExciteSingleParticipants:ProjM="<<Pprojectile.mag()<<G4endl; 03538 #endif 03539 projectile->Set4Momentum(Pprojectile); 03540 return true; 03541 } // End of ExciteSingleDiffParticipants
G4QHadronVector * G4QIonIonCollision::Fragment | ( | ) |
Definition at line 1557 of file G4QIonIonCollision.cc.
References G4QEnvironment::AddQuasmon(), Breeder(), FatalException, G4QEnvironment::Fragment(), G4cerr, G4cout, G4endl, G4Exception(), G4QHadron::Get4Momentum(), G4QNucleus::GetA(), G4QPDGCode::GetMass(), G4QNucleus::GetPDG(), G4QHadron::GetPDGCode(), G4QHadron::GetQC(), G4QNucleus::GetQCZNS(), G4QCHIPSWorld::GetQParticle(), G4QChipolino::GetQPDG1(), G4QChipolino::GetQPDG2(), G4QContent::GetSPDGCode(), G4QNucleus::GetZ(), G4QParticle::MinMassOfFragm(), G4Quasmon::Set4Momentum(), and G4Quasmon::SetQC().
Referenced by G4QInelastic::PostStepDoIt().
01558 { // This is the member function fragmenting Strings & Quasmons (in nuclear matter) 01559 #ifdef debug 01560 G4cout<<"*******>G4QIonIonCollision::Fragment: ***Called***, Res="<<theResult<<G4endl; 01561 #endif 01562 G4int striNum=strings.size(); // Find out if there're strings 01563 G4int hadrNum=theResult->size(); // Find out if there're hadrons 01564 #ifdef edebug 01565 G4int nQm=theQuasmons.size(); 01566 G4LorentzVector totLS4M=theProjNucleus.Get4Momentum()+theTargNucleus.Get4Momentum(); //LS 01567 G4int totChg=theProjNucleus.GetZ()+theTargNucleus.GetZ(); 01568 G4int totBaN=theTargNucleus.GetA()+theTargNucleus.GetA(); 01569 G4cout<<"-EMCLS-G4QIn::Fragment: CHECKRecovery, #ofS="<<striNum<<", Nuc4M(E=M)="<<totLS4M 01570 <<",#Q="<<nQm<<",#H="<<hadrNum<<G4endl; 01571 for(G4int i=0; i < striNum; i++) 01572 { 01573 G4LorentzVector strI4M=strings[i]->Get4Momentum(); 01574 totLS4M+=strI4M; 01575 G4int sChg=strings[i]->GetCharge(); 01576 totChg+=sChg; 01577 G4int sBaN=strings[i]->GetBaryonNumber(); 01578 totBaN+=sBaN; 01579 G4cout<<"-EMCLS-G4QIonIonCollision::Fragment:String#"<<i<<", 4M="<<strI4M<<", M=" 01580 <<strI4M.m()<<", C="<<sChg<<", B="<<sBaN<<G4endl; 01581 } 01582 for(G4int i=0; i < nQm; i++) 01583 { 01584 G4LorentzVector hI4M=theQuasmons[i]->Get4Momentum(); 01585 totLS4M+=hI4M; 01586 G4int hChg=theQuasmons[i]->GetCharge(); 01587 totChg+=hChg; 01588 G4int hBaN=theQuasmons[i]->GetBaryonNumber(); 01589 totBaN+=hBaN; 01590 G4cout<<"-EMCLS-G4QIonIonCollision::Fragment: Quasmon#"<<i<<", 4M="<<hI4M<<", C="<<hChg 01591 <<", B="<<hBaN<<G4endl; 01592 } 01593 for(G4int i=0; i < hadrNum; i++) 01594 { 01595 G4LorentzVector hI4M=(*theResult)[i]->Get4Momentum(); 01596 totLS4M+=hI4M; 01597 G4int hChg=(*theResult)[i]->GetCharge(); 01598 totChg+=hChg; 01599 G4int hBaN=(*theResult)[i]->GetBaryonNumber(); 01600 totBaN+=hBaN; 01601 G4cout<<"-EMCLS-G4QIn::Fragment:H#"<<i<<",4M="<<hI4M<<",C="<<hChg<<",B="<<hBaN<<G4endl; 01602 } 01603 #endif 01604 #ifdef debug 01605 G4cout<<"***>G4QIonIonCollision::Fragm: #OfStr="<<striNum<<", #OfRes="<<hadrNum<<G4endl; 01606 #endif 01607 if(!striNum && hadrNum) // Quasi-elastic or decoupled p 01608 { 01609 #ifdef debug 01610 G4cout<<"***>G4QIonIonCollision::Fragm:**Quasi-Elastic**, #OfResult="<<hadrNum<<G4endl; 01611 #endif 01612 return theResult; 01613 } 01614 else if(striNum) Breeder(); // Strings fragmentation 01615 else // No strings, make HadrNucleus 01616 { 01617 if(hadrNum) 01618 { 01619 for(G4int ih=0; ih<hadrNum; ih++) delete (*theResult)[ih]; 01620 theResult->clear(); 01621 } 01622 G4LorentzVector rp4M=theProjNucleus.Get4Momentum(); // Nucleus 4-momentum in LS 01623 G4int rpPDG=theProjNucleus.GetPDG(); // Nuclear PDG 01624 G4QHadron* resPNuc = new G4QHadron(rpPDG,rp4M); // Nucleus -> Hadron 01625 theResult->push_back(resPNuc); // Fill the residual nucleus 01626 G4LorentzVector rt4M=theTargNucleus.Get4Momentum(); // Nucleus 4-momentum in LS 01627 G4int rtPDG=theTargNucleus.GetPDG(); // Nuclear PDG 01628 G4QHadron* resTNuc = new G4QHadron(rtPDG,rt4M); // Nucleus -> Hadron 01629 theResult->push_back(resTNuc); // Fill the residual nucleus 01630 } 01631 G4int nQuas=theQuasmons.size(); // Size of the Quasmon OUTPUT 01632 G4int theRS=theResult->size(); // Size of Hadron Output by now 01633 #ifdef debug 01634 G4cout<<"***>G4QIonIonCollision::Fragm:beforeEnv, #OfQ="<<nQuas<<",#OfR="<<theRS<<G4endl; 01635 #endif 01636 if(nQuas && theRS) 01637 { 01638 01639 G4QHadron* resNuc = (*theResult)[theRS-1]; // Pointer to Residual Nucleus 01640 G4LorentzVector resNuc4M = resNuc->Get4Momentum(); // 4-Momentum of the Nucleuz 01641 G4int resNucPDG= resNuc->GetPDGCode(); // PDG Code of the Nucleus 01642 G4QNucleus theEnv(resNucPDG); // NucleusHadron->NucleusAtRest 01643 delete resNuc; // Delete resNucleus as aHadron 01644 theResult->pop_back(); // Exclude the nucleus from HV 01645 --theRS; // Reduce the OUTPUT by theNucl 01646 #ifdef pdebug 01647 G4cout<<"G4QIonIonCollision::Fragm:#OfRemainingHadron="<<theRS<<", A="<<theEnv<<G4endl; 01648 #endif 01649 // Now we need to be sure that the compound nucleus is heavier than the Ground State 01650 for(G4int j=theRS-1; j>-2; --j) // try to reach M_compound>M_GS 01651 { 01652 G4LorentzVector qsum4M=resNuc4M; // Proto compound 4-momentum 01653 G4QContent qsumQC=theEnv.GetQCZNS(); // Proto compound Quark Content 01654 #ifdef pdebug 01655 G4cout<<"G4QIonIonCollis::Fragm: rN4M"<<qsum4M<<qsum4M.m()<<",rNQC="<<qsumQC<<G4endl; 01656 #endif 01657 G4Quasmon* firstQ=0; // Prototype of theFirstQuasmon 01658 G4LorentzVector first4M; // Proto of the FirstQuasmon 4M 01659 G4QContent firstQC; // Proto of the FirstQuasmon QC 01660 for(G4int i=0; i<nQuas; ++i) // LOOP over Quasmons 01661 { 01662 G4Quasmon* curQuasm=theQuasmons[i]; // current Quasmon 01663 G4LorentzVector cur4M=curQuasm->Get4Momentum(); // 4-Mom of the Quasmon 01664 G4QContent curQC=curQuasm->GetQC(); // Quark Content of the Quasmon 01665 qsum4M+=cur4M; // Add quasmon's 4-momentum 01666 qsumQC+=curQC; // Add quasmon's Quark Content 01667 #ifdef pdebug 01668 G4cout<<"G4QIonIonCollis::Fragm: Q#"<<i<<", QQC="<<curQC<<", sQC="<<qsumQC<<G4endl; 01669 #endif 01670 if(!i) // Remember 1-st for correction 01671 { 01672 firstQ =curQuasm; 01673 first4M=cur4M; 01674 firstQC=curQC; 01675 } 01676 } 01677 G4int miPDG=qsumQC.GetSPDGCode(); // PDG of minM of hadron/fragm. 01678 G4double gsM=0.; // Proto minM of had/frag forQC 01679 if(miPDG == 10) 01680 { 01681 G4QChipolino QCh(qsumQC); // define TotNuc as a Chipolino 01682 gsM=QCh.GetQPDG1().GetMass()+QCh.GetQPDG2().GetMass(); // Sum of Hadron Masses 01683 //gsM=theWorld->GetQParticle(QCh.GetQPDG1())->MinMassOfFragm() + 01684 // theWorld->GetQParticle(QCh.GetQPDG2())->MinMassOfFragm(); 01685 } 01686 // @@ it is not clear, why it does not work ?! 01687 //else if(miPDG>80000000) // Compound Nucleus 01688 //{ 01689 // G4QNucleus rtN(qsumQC); // CreatePseudoNucl for totComp 01690 // gsM=rtN.GetGSMass(); // MinMass of residQ+(Env-ParC) 01691 //} 01692 else if(miPDG < 80000000 && std::abs(miPDG)%10 > 2) 01693 gsM=theWorld->GetQParticle(G4QPDGCode(miPDG))->MinMassOfFragm(); 01694 else gsM=G4QPDGCode(miPDG).GetMass(); // minM of hadron/fragm. for QC 01695 G4double reM=qsum4M.m(); // real mass of the compound 01696 #ifdef pdebug 01697 G4cout<<"G4QIonIonCollision::Fragm: PDG="<<miPDG<<", rM="<<reM<<",GSM="<<gsM<<G4endl; 01698 #endif 01699 if(reM > gsM) break; // CHIPS can be called 01700 if(j > -1) // Can try to add hadrons to Q0 01701 { 01702 G4QHadron* cH = (*theResult)[j]; // Pointer to the last Hadron 01703 G4LorentzVector h4M = cH->Get4Momentum(); // 4-Momentum of the Hadron 01704 G4QContent hQC = cH->GetQC(); // QC of the Hadron 01705 firstQ->Set4Momentum(first4M+h4M); // Update the Quasmon's 4-Mom 01706 firstQ->SetQC(firstQC+hQC); // Update the Quasmon's QCont 01707 delete cH; // Delete the Hadron 01708 theResult->pop_back(); // Exclude the hadron from HV 01709 #ifdef pdebug 01710 G4cout<<"G4QInel::Fragm: H#"<<j<<", hQC="<<hQC<<",hPDG="<<cH->GetPDGCode()<<G4endl; 01711 #endif 01712 } 01713 else 01714 { 01715 G4cerr<<"***G4QIonIonCollision::Fra:PDG="<<miPDG<<",M="<<reM<<",GSM="<<gsM<<G4endl; 01716 G4Exception("G4QIonIonCollision::Fragment:","27",FatalException,"Can'tRecoverGSM"); 01717 } 01718 } 01719 G4double nucE=resNuc4M.e(); // Total energy of the nuclEnv 01720 if(nucE<1.E-12) nucE=0.; // Computer accuracy safety 01721 G4ThreeVector nucVel(0.,0.,0.); // Proto of the NucleusVelocity 01722 G4QHadronVector* output=0; // NucleusFragmentation Hadrons 01723 G4QEnvironment* pan= new G4QEnvironment(theEnv); // ---> DELETED --->----------+ 01724 #ifdef pdebug 01725 G4cout<<"G4QIonIonCollision::Fragm: nucE="<<nucE<<", nQ="<<nQuas<<G4endl; // | 01726 #endif 01727 if(nucE) nucVel=resNuc4M.vect()/nucE; // The NucleusVelocity | 01728 for(G4int i=0; i<nQuas; ++i) // LOOP over Quasmons | 01729 { // | 01730 G4Quasmon* curQuasm=theQuasmons[i]; // current Quasmon | 01731 if(nucE) curQuasm->Boost(-nucVel); // Boost it to CMS of Nucleus | 01732 pan->AddQuasmon(curQuasm); // Fill the predefined Quasmon| 01733 #ifdef pdebug 01734 G4LorentzVector cQ4M=curQuasm->Get4Momentum(); // Just for printing | 01735 G4cout<<"G4QIonIonCollision::Fragment: Quasmon# "<<i<<" added, 4M="<<cQ4M<<G4endl;//| 01736 #endif 01737 } // | 01738 try // | 01739 { // | 01740 delete output; // | 01741 output = pan->Fragment();// DESTROYED after theHadrons are transferred to theResult | 01742 } // | | 01743 catch (G4QException& error) // | | 01744 { // | | 01745 G4cerr<<"***G4QIonIonCollision::Fragment: G4QE Exception is catched"<<G4endl; // | | 01746 // G4Exception("G4QIonIonCollision::Fragment:","27",FatalException,"CHIPSCrash");// | | 01747 G4Exception("G4QIonIonCollision::Fragment()", "HAD_CHPS_0027", 01748 FatalException, "CHIPSCrash"); 01749 } // | | 01750 delete pan; // Delete the Nuclear Environment <-----<--+-+ 01751 if(output) // Output exists | 01752 { // | 01753 G4int nOut=output->size(); // #ofHadrons in the Nuclear Fragmentation | 01754 for(G4int j=0; j<nOut; j++) // LOOP over Hadrons transferring to LS | 01755 { // | 01756 G4QHadron* curHadron=(*output)[j]; // Hadron from the nucleus fragmentation | 01757 if(nucE) curHadron->Boost(nucVel); // Boost it back to Laboratory System | 01758 theResult->push_back(curHadron); // Transfer it to the result | 01759 } // | 01760 delete output; // Delete the OUTPUT <-----<-----<-----<---+ 01761 } 01762 } 01763 else if(!striNum) G4cout<<"-Warning-G4QIonIonCollision::Fragment:NothingWasDone"<<G4endl; 01764 #ifdef debug 01765 G4cout<<"=--=>G4QIonIonCollision::Fragment: Final #OfResult="<<theResult->size()<<G4endl; 01766 #endif 01767 G4int nQ =theQuasmons.size(); 01768 if(nQ) theQuasmons.clear(); // @@ Not necesary ? 01769 #ifdef edebug 01770 G4LorentzVector f4M(0.,0.,0.,0.); // Sum of the Result in LS 01771 G4int fCh=totChg; 01772 G4int fBN=totBaN; 01773 G4int nHd=theResult->size(); 01774 G4cout<<"-EMCLS-G4QIonIonCollision::Fragment: #ofHadr="<<nHd<<",#OfQuasm="<<nQ<<G4endl; 01775 for(G4int i=0; i<nHd; i++) 01776 { 01777 G4LorentzVector hI4M=(*theResult)[i]->Get4Momentum(); 01778 f4M+=hI4M; 01779 G4int hChg=(*theResult)[i]->GetCharge(); 01780 fCh-=hChg; 01781 G4int hBaN=(*theResult)[i]->GetBaryonNumber(); 01782 fBN-=hBaN; 01783 G4cout<<"-EMCLS-G4QIonIonCollision::Fragment: Hadron#"<<i<<", 4M="<<hI4M<<", PDG=" 01784 <<(*theResult)[i]->GetPDGCode()<<", C="<<hChg<<", B="<<hBaN<<G4endl; 01785 } 01786 G4cout<<"-EMCLS-G4QInel::Fragm:, r4M="<<f4M-totLS4M<<", rC="<<fCh<<",rB="<<fBN<<G4endl; 01787 #endif 01788 return theResult; 01789 } // End of fragmentation
G4ThreeVector G4QIonIonCollision::GaussianPt | ( | G4double | widthSquare, | |
G4double | maxPtSquare | |||
) | const [protected] |
Definition at line 3572 of file G4QIonIonCollision.cc.
References G4cout, G4endl, G4UniformRand, and sqr().
Referenced by ExciteDiffParticipants(), and ExciteSingDiffParticipants().
03573 { 03574 #ifdef debug 03575 G4cout<<"G4QIonIonCollision::GaussianPt: w^2="<<widthSq<<",maxPt2="<<maxPtSquare<<G4endl; 03576 #endif 03577 G4double pt2=0.; 03578 G4double rm=maxPtSquare/widthSq; // Negative 03579 if(rm>-.01) pt2=widthSq*(std::sqrt(1.-G4UniformRand()*(1.-sqr(1.+rm)))-1.); 03580 else pt2=widthSq*std::log(1.-G4UniformRand()*(1.-std::exp(rm))); 03581 pt2=std::sqrt(pt2); // It is not pt2, but pt now 03582 G4double phi=G4UniformRand()*twopi; 03583 return G4ThreeVector(pt2*std::cos(phi),pt2*std::sin(phi),0.); 03584 } // End of GaussianPt
G4bool G4QIonIonCollision::IsSingleDiffractive | ( | ) | [inline, protected] |
Definition at line 76 of file G4QIonIonCollision.hh.
References G4UniformRand.
Referenced by G4QIonIonCollision().
00076 {G4bool r=false; if(G4UniformRand()<1.) r=true; return r;}
Definition at line 3651 of file G4QIonIonCollision.cc.
References G4cout, and G4endl.
Referenced by Breeder().
03652 { 03653 #ifdef debug 03654 G4cout<<"G4QIonIonCollision::ReducePair: **Called** P1="<<P1<<", P2="<<P2<<G4endl; 03655 #endif 03656 G4int P11=P1/10; 03657 G4int P12=P1%10; 03658 G4int P21=P2/10; 03659 G4int P22=P2%10; 03660 if (P11==P21) return std::make_pair(P12,P22); 03661 else if(P11==P22) return std::make_pair(P12,P21); 03662 else if(P12==P21) return std::make_pair(P11,P22); 03663 else if(P12==P22) return std::make_pair(P11,P21); 03664 //#ifdef debug 03665 G4cout<<"-Warning-G4QIonIonCollision::ReducePair:**Failed**,P1="<<P1<<",P2="<<P2<<G4endl; 03666 //#endif 03667 return std::make_pair(0,0); // Reduction failed 03668 }
void G4QIonIonCollision::SetParameters | ( | G4int | nC, | |
G4double | strTens, | |||
G4double | tubeDens, | |||
G4double | SigPt | |||
) | [static] |
Definition at line 3543 of file G4QIonIonCollision.cc.
03544 { 03545 nCutMax = nC; // max number of pomeron cuts 03546 stringTension = stT; // string tension for absorbed energy 03547 tubeDensity = tbD; // Flux Tube Density of nuclear nucleons 03548 widthOfPtSquare = -2*SigPt*SigPt; // width^2 of pt for string excitation 03549 }
Definition at line 3586 of file G4QIonIonCollision.cc.
References FatalException, G4cerr, G4endl, and G4Exception().
Referenced by Breeder().
03587 { 03588 if (PDG1 < 7 && PDG1 > 0 && PDG2 < 7 && PDG2 > 0) // Sum up two Q in DiQ (S=0) 03589 { 03590 if(PDG1 > PDG2) return PDG1*1000+PDG2*100+1; 03591 else return PDG2*1000+PDG1*100+1; 03592 } 03593 else if (PDG1 >-7 && PDG1 < 0 && PDG2 >-7 && PDG2 < 0) // Sum up two AQ in ADiQ (S=0) 03594 { 03595 if(-PDG1 > -PDG2) return PDG1*1000+PDG2*100-1; 03596 else return PDG2*1000+PDG1*100-1; 03597 } 03598 else if (PDG1 <-99 && PDG2 < 7 && PDG2 > 0) // Sum up Q and ADiQ in AQ 03599 { 03600 G4int PDG=-PDG1/100; 03601 if(PDG2==PDG/10) return -PDG%10; 03602 if(PDG2==PDG%10) return -PDG/10; 03603 else 03604 { 03605 G4cerr<<"***4QFragmentation::SumPartonPDG: PDG1="<<PDG1<<", PDG2="<<PDG2<<G4endl; 03606 G4Exception("G4QIonIonCollision::SumPartonPDG:","72",FatalException,"Q&ADiQnoMatch"); 03607 } 03608 } 03609 else if (PDG2 <-99 && PDG1 < 7 && PDG1 > 0) // Sum up ADiQ and Q in AQ 03610 { 03611 G4int PDG=-PDG2/100; 03612 if(PDG1==PDG/10) return -PDG%10; 03613 if(PDG1==PDG%10) return -PDG/10; 03614 else 03615 { 03616 G4cerr<<"***4QFragmentation::SumPartonPDG: PDG1="<<PDG1<<", PDG2="<<PDG2<<G4endl; 03617 G4Exception("G4QIonIonCollision::SumPartonPDG:","72",FatalException,"ADiQ&QnoMatch"); 03618 } 03619 } 03620 else if (PDG1 > 99 && PDG2 >-7 && PDG2 < 0) // Sum up DiQ and AQ in Q 03621 { 03622 G4int PDG=PDG1/100; 03623 if(PDG2==-PDG/10) return PDG%10; 03624 if(PDG2==-PDG%10) return PDG/10; 03625 else 03626 { 03627 G4cerr<<"***4QFragmentation::SumPartonPDG: PDG1="<<PDG1<<", PDG2="<<PDG2<<G4endl; 03628 G4Exception("G4QIonIonCollision::SumPartonPDG:","72",FatalException,"DiQ&AQnoMatch"); 03629 } 03630 } 03631 else if (PDG2 > 99 && PDG1 >-7 && PDG1 < 0) // Sum up AQ and DiQ in Q 03632 { 03633 G4int PDG=PDG2/100; 03634 if(PDG1==-PDG/10) return PDG%10; 03635 if(PDG1==-PDG%10) return PDG/10; 03636 else 03637 { 03638 G4cerr<<"***G4QIonIonCollision::SumPartonPDG: PDG1="<<PDG1<<", PDG2="<<PDG2<<G4endl; 03639 G4Exception("G4QIonIonCollision::SumPartonPDG:","72",FatalException,"AQ&DiQnoMatch"); 03640 } 03641 } 03642 else 03643 { 03644 G4cerr<<"***G4QIonIonCollision::SumPartonPDG: PDG1="<<PDG1<<", PDG2="<<PDG2<<G4endl; 03645 G4Exception("G4QIonIonCollision::SumPartonPDG:","72",FatalException,"Can'tSumUpParts"); 03646 } 03647 return 0; 03648 } // End of SumPartonPDG
void G4QIonIonCollision::SwapPartons | ( | ) | [protected] |
Definition at line 3903 of file G4QIonIonCollision.cc.
References DBL_MAX, G4cout, G4endl, G4QParton::Get4Momentum(), G4QParton::GetPDGCode(), and G4QParton::GetType().
Referenced by G4QIonIonCollision().
03904 { 03905 static const G4double baryM=800.; // Mass excess for baryons 03906 G4QStringVector::iterator ist; 03907 for(ist = strings.begin(); ist < strings.end(); ist++) 03908 { 03909 G4LorentzVector cS4M=(*ist)->Get4Momentum(); 03910 G4double cSM2=cS4M.m2(); // Squared mass of the String 03911 if(cSM2<0.1) // Parton Swapping is needed 03912 { 03913 G4QParton* cLeft=(*ist)->GetLeftParton(); // Current String Left Parton 03914 G4QParton* cRight=(*ist)->GetRightParton(); // Current String Right Parton 03915 G4int cLPDG=cLeft->GetPDGCode(); 03916 G4int cRPDG=cRight->GetPDGCode(); 03917 G4LorentzVector cL4M=cLeft->Get4Momentum(); 03918 G4LorentzVector cR4M=cRight->Get4Momentum(); 03919 G4int cLT=cLeft->GetType(); 03920 G4int cRT=cRight->GetType(); 03921 G4QStringVector::iterator sst; // Selected partner string 03922 G4QStringVector::iterator pst; // LOOP iterator 03923 G4double maxM=-DBL_MAX; // Swapping providing the max mass 03924 G4int sDir=0; // Selected direction of swapping 03925 for(pst = strings.begin(); pst < strings.end(); pst++) if(pst != ist) 03926 { 03927 G4QParton* pLeft=(*pst)->GetLeftParton(); // Partner String Left Parton 03928 G4QParton* pRight=(*pst)->GetRightParton(); // Partner String Right Parton 03929 G4int pLPDG=pLeft->GetPDGCode(); 03930 G4int pRPDG=pRight->GetPDGCode(); 03931 G4LorentzVector pL4M=pLeft->Get4Momentum(); 03932 G4LorentzVector pR4M=pRight->Get4Momentum(); 03933 G4int pLT=pLeft->GetType(); 03934 G4int pRT=pRight->GetType(); 03935 G4double LM=0.; 03936 G4double RM=0.; 03937 if(((cLPDG<-7 || (cLPDG>0 && cLPDG< 7) ) && (pLPDG<-7 || (pLPDG>0 && pLPDG< 7) ))|| 03938 ((cLPDG> 7 || (cLPDG<0 && cLPDG>-7) ) && (pLPDG> 7 || (pLPDG<0 && cLPDG>-7) ))) 03939 { 03940 G4double pLM2=(cL4M+pR4M).m2(); // new partner M2 03941 G4double cLM2=(cR4M+pL4M).m2(); // new partner M2 03942 if(pLM2>0. && cLM2>0.) 03943 { 03944 G4double pLM=std::sqrt(pLM2); 03945 if(cLT+pRT==3) pLM-=baryM; 03946 G4double cLM=std::sqrt(cLM2); 03947 if(cRT+pLT==3) cLM-=baryM; 03948 LM=std::min(pLM2,cLM2); 03949 } 03950 } 03951 if(((cRPDG<-7 || (cRPDG>0 && cRPDG< 7) ) && (pRPDG<-7 || (pRPDG>0 && pRPDG< 7) ))|| 03952 ((cRPDG> 7 || (cRPDG<0 && cRPDG>-7) ) && (pRPDG> 7 || (pRPDG<0 && cRPDG>-7) )) ) 03953 { 03954 G4double pRM2=(cR4M+pL4M).m2(); // new partner M2 03955 G4double cRM2=(cL4M+pR4M).m2(); // new partner M2 03956 if(pRM2>0. && cRM2>0.) 03957 { 03958 G4double pRM=std::sqrt(pRM2); 03959 if(cRT+pLT==3) pRM-=baryM; 03960 G4double cRM=std::sqrt(cRM2); 03961 if(cLT+pRT==3) cRM-=baryM; 03962 RM=std::min(pRM,cRM); 03963 } 03964 } 03965 G4int dir=0; 03966 G4double sM=0.; 03967 if( LM && LM > RM ) 03968 { 03969 dir= 1; 03970 sM=LM; 03971 } 03972 else if(RM) 03973 { 03974 dir=-1; 03975 sM=RM; 03976 } 03977 if(sM > maxM) 03978 { 03979 sst=pst; 03980 maxM=sM; 03981 sDir=dir; 03982 } 03983 } 03984 if(sDir) 03985 { 03986 G4QParton* pLeft=(*sst)->GetLeftParton(); // Partner String Left Parton 03987 G4QParton* pRight=(*sst)->GetRightParton(); // Partner String Right Parton 03988 G4QParton* swap=pLeft; // Prototype remember the partner Left 03989 if(sDir>0) // swap left partons 03990 { 03991 (*sst)->SetLeftParton(cLeft); 03992 (*ist)->SetLeftParton(swap); 03993 } 03994 else 03995 { 03996 swap=pRight; 03997 (*sst)->SetRightParton(cRight); 03998 (*ist)->SetRightParton(swap); 03999 } 04000 } 04001 #ifdef debug 04002 else G4cout<<"***G4QIonIonCollision::SwapPartons:**Failed**,cLPDG="<<cLPDG<<",cRPDG=" 04003 <<cRPDG<<",-->cM2="<<cSM2<<G4endl; 04004 #endif 04005 04006 } 04007 } 04008 }