Geant4 9.6.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
Public Member Functions | Static Public Member Functions | Protected Member Functions | List of all members
G4QIonIonCollision Class Reference

#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 ()
 

Detailed Description

Definition at line 59 of file G4QIonIonCollision.hh.

Constructor & Destructor Documentation

◆ G4QIonIonCollision()

G4QIonIonCollision::G4QIonIonCollision ( G4QNucleus &  pNucleus,
const G4QNucleus &  tNucleus 
)

Definition at line 62 of file G4QIonIonCollision.cc.

63{
64 static const G4double mProt= G4QPDGCode(2212).GetMass(); // Mass of proton
65 static const G4double mPi0= G4QPDGCode(111).GetMass(); // Mass of Pi0
66 theWorld= G4QCHIPSWorld::Get(); // Get a pointer to the CHIPS World
67 theResult = new G4QHadronVector; // Define theResultingHadronVector
68 G4bool stringsInitted=false; // Strings are initiated
69 G4LorentzRotation toZ; // Lorentz Transformation to the projectileSystem
70 G4LorentzVector proj4M=pNucleus.Get4Momentum(); // Projectile nucleus 4-momentum in LS
71 //G4double projM=pNucleus.GetGSMass(); // Ground state mass of the projectile nucleus
72 G4double targM=tNucleus.GetGSMass(); // Ground state mass of the target nucleus
73#ifdef edebug
74 G4cout<<"G4QIn::Constr:*Called*,pA="<<pNucleus<<proj4M<<",tA="<<tNucleus<<targM<<G4endl;
75#endif
76 G4int pZ=pNucleus.GetZ();
77 G4int pN=pNucleus.GetN();
78 G4int pA=pZ+pN;
79 G4int pPDG=90000000+pZ*1000+pN;
80 G4int tZ=tNucleus.GetZ();
81 G4int tN=tNucleus.GetN();
82 G4int tA=tZ+tN;
83 G4int tPDG=90000000+tZ*1000+tN;
84 toZ.rotateZ(-proj4M.phi());
85 toZ.rotateY(-proj4M.theta());
86 G4LorentzVector zProj4M=toZ*proj4M; // Proj 4-momentum in LS rotated to Z axis
87 pNucleus.Set4Momentum(zProj4M); // Now the projectile nucleus moves along Z axis
88#ifdef edebug
89 G4int totChg=pZ+tZ; // Charge of the projectile+target for the CHECK
90 G4int totBaN=pA+tA; // Baryon Number of Proj+Targ for CHECK
91 G4LorentzVector tgLS4M(0.,0.,0.,targM); // Target 4-momentum in LS
92 G4LorentzVector totLS4M=proj4M+tgLS4M; // Total 4-momentum in LS
93 G4LorentzVector totZLS4M=zProj4M+tgLS4M;// Total 4-momentum in LS with momentum along Z
94 G4cout<<"-EMC-G4QIonIonCollision::Constr: tLS4M="<<totLS4M<<",tZLS4M="<<totZLS4M<<G4endl;
95 // === From nere all consideration is made in the rotated LS frame (proj is along Z) ===
96#endif
97 G4LorentzRotation toLab(toZ.inverse()); // Lorentz Transfornation "ZLS"->LS (at the end)
98 G4QInteractionVector theInteractions; // A vector of interactions (taken from the Body)
99 G4QPartonPairVector thePartonPairs; // The parton pairs (taken from the Body)
100 G4int ModelMode=SOFT; // The current model type (taken from the Body)
101 theTargNucleus=G4QNucleus(tZ,tN); // Create theTargNucleus to Move From LS to CM
102 theTargNucleus.InitByPDG(tPDG); // Reinit the Nucleus for the new Attempt?
103 theTargNucleus.Init3D(); // 3D-initialisation(nucleons)ofTheTGNucleusClone
104#ifdef debug
105 G4cout<<"G4QIonIonCollision::Constr:TargNuclMom="<<theTargNucleus.Get4Momentum()<<G4endl;
106#endif
107 theProjNucleus=G4QNucleus(pZ,pN); // Create theProjNucleus to Move From LS to CM
108 theProjNucleus.InitByPDG(pPDG); // Reinit the Nucleus for the new Attempt?
109 theProjNucleus.Init3D(); // 3D-initialisation(nucleons)ofThePrNucleusClone
110#ifdef debug
111 G4cout<<"G4QIonIonCollision::Constr:ProjNuclMom="<<theProjNucleus.Get4Momentum()<<G4endl;
112#endif
113#ifdef edebug
114 G4LorentzVector sumP1=theProjNucleus.GetNucleons4Momentum();// Sum ofPrNucleons4M inRotLS
115 G4LorentzVector sumT1=theTargNucleus.GetNucleons4Momentum();// Sum ofTgNucleons4M inRotLS
116 G4cout<<"-EMC-G4QIonIonCollision::Construct: pNuc4M="<<sumP1<<", tNuc4M="<<sumT1<<G4endl;
117#endif
118 G4ThreeVector theCMVelocity(0.,0.,0.); // Prototype of the "CM" velocity
119 G4ThreeVector theALVelocity(0.,0.,0.); // Prototype of "CMAntiLab" velocity
120 // Very important! Here the projectile 4M is recalculated in the ZLS (previously in LS)
121 proj4M = pNucleus.Get4Momentum(); // 4-mom of theProjectileHadron inZLS
122 G4double pz_per_projectile = proj4M.pz()/pA; // Momentum per nucleon
123 G4double e_per_projectile = proj4M.e()/pA+targM/tA;// Add MassOfTargetNucleon
124 G4double vz = pz_per_projectile/e_per_projectile; // Speed of the "oneNuclenCM"
125#ifdef debug
126 G4cout<<"G4QIonIonCollision::Construct: (КщеЯ)Projectile4M="<<proj4M<<", Vz="<<vz
127 <<", pA="<<pA<<", pE="<<e_per_projectile<<G4endl;
128#endif
129 theCMVelocity.setZ(vz); // CM (w/r to one nucleon) velosity
130 theProjNucleus.DoLorentzBoost(-theCMVelocity); // BoostProjNucleus to "CM"
131 theTargNucleus.DoLorentzBoost(-theCMVelocity); // BoostProjNucleus to "CM"
132 G4LorentzVector prN4M=theProjNucleus.Get4Momentum();
133 G4double avz=prN4M.pz()/prN4M.e(); // Speed of AntiLabSys in CMS
134 theALVelocity.setZ(avz); // CM (w/r to one nucleon) velosity
135#ifdef edebug
136 G4LorentzVector sumP2=theProjNucleus.GetNucleons4Momentum();// Sum of Nucleons4M in RotCM
137 G4LorentzVector sumT2=theTargNucleus.GetNucleons4Momentum();// Sum of Nucleons4M in RotCM
138 G4cout<<"-EMC-G4QIonIonCollision::Construct: Boosted, vCM="<<vz<<", vAL="<<avz<<", tN4M="
139 <<sumT2<<", pN4M="<<sumP2<<G4endl;
140#endif
141 G4int maxCuts = 7; // @@ Can be reduced for low energies
142 //
143 // >>----------->> Find collisions meeting collision conditions and the NN interaction XS
144 //
145 G4double outerRadius = theProjNucleus.GetOuterRadius()+theTargNucleus.GetOuterRadius();
146 G4QProbability theProbability(2212); // *** proj/targ nucleons ***
147 // Clean up all previous interactions and reset the counters
148#ifdef debug
149 G4cout<<"G4QIonIonCollision::Construct: theIntSize="<<theInteractions.size()<<G4endl;
150#endif
151 G4int attCnt=-1;
152 G4int maxAtt=27;
153 // From here a loop over nucleons of the projectile Nucleus (@@ Z-sorting isn't implem!!)
154 //
155 G4QHadron* pNucleon; // Proto of the Projectile Nucleon pointer
156 G4QHadron* tNucleon; // Proto of the Target Nucleon pointer
157 G4QNucleus curProjNucleus;
158 G4QNucleus curTargNucleus;
159 while(!theInteractions.size() && ++attCnt < maxAtt)// TillAtLeastOneInteractionIsCreated
160 {
161 // std::for_each(theInteractions.begin(),theInteractions.end(), DeleteQInteraction());
162 // Here we need to clean up the ProjNucleon and the TargNucleon in the Interactions !
163 G4int nint=theInteractions.size(); // For the 1st attempt should be zero
164 for(G4int ni=0; ni < nint; ++ni)
165 {
166 G4QInteraction* curInt = theInteractions[ni];
167 delete curInt->GetProjectile();
168 delete curInt->GetTarget();
169 delete curInt;
170 }
171 theInteractions.clear();
172 // Now we need to make a copy of the projectile and the target nuclei with 3D info
173 if(attCnt) // This is not theFirstAttempt->Clean
174 {
175 curProjNucleus.DeleteNucleons();
176 curTargNucleus.DeleteNucleons();
177 }
178 curProjNucleus = theProjNucleus;
179 curTargNucleus = theTargNucleus;
180 // choose random impact parameter
181 std::pair<G4double, G4double> theImpactParameter;
182 theImpactParameter = curProjNucleus.ChooseImpactXandY(outerRadius);
183 G4double impactX = theImpactParameter.first;
184 G4double impactY = theImpactParameter.second;
185#ifdef debug
186 G4cout<<"G4QIonIonCollision::Con:At#="<<attCnt<<",X="<<impactX<<",Y="<<impactY<<G4endl;
187#endif
188 curProjNucleus.StartLoop(); // Prepare Loop ovder nucleons
189 G4int pnCnt=0;
190 G4int pnA=curProjNucleus.GetA();
191#ifdef debu
192 G4cout<<"G4QIonIonCollision::Constr: Before the WHILE over pNucleons,pA="<<pnA<<G4endl;
193#endif
194 while((pNucleon=curProjNucleus.GetNextNucleon()) && pnCnt < pnA) // @@ can be for LOOP?
195 {
196 ++pnCnt;
197 G4QInteractionVector curInteractions; // A temporary vector of interactions
198 G4LorentzVector pNuc4M=pNucleon->Get4Momentum();
199 G4ThreeVector pNucR=pNucleon->GetPosition(); // Position of the pNucleon WRTo pNucl
200 G4double pImpX = impactX+pNucR.x();
201 G4double pImpY = impactY+pNucR.y();
202 G4double prEn=proj4M.e(); // For mesons
203 G4int proB=pNucleus.GetBaryonNumber();
204 if (proB>0) prEn-=pNucleus.GetMass(); // For baryons
205 else if(proB<0) prEn+=mProt; // For anti-baryons
206#ifdef debug
207 G4double impactUsed = 0.;
208 G4cout<<"G4QIonIonCollision::Construct: estimated energy, prEn="<<prEn<<G4endl;
209#endif
210 G4int totalCuts = 0;
211 // @@ This is a fake (random) s calculation @@ can be done inside the TARG-while
212 G4int tnA=curTargNucleus.GetA();
213 G4double pImp=std::sqrt(pImpX*pImpX+pImpY*pImpY);
214 G4double eflen=curTargNucleus.GetThickness(pImp); // EffectiveLength
215 maxEn=eflen*stringTension; // maxAbsorbedEnergy in IndUnits=MeV
216 maxNuc=static_cast<int>(eflen*tubeDensity+0.5);// max #0f involved nuclear nucleons
217#ifdef debug
218 G4cout<<"G4QIonIonCollision::Con:pE="<<prEn<<" < mE="<<maxEn<<",mN="<<maxNuc<<G4endl;
219#endif
220 if(prEn < maxEn) // Create DIN interaction and go out
221 {
222 G4QHadron* aProjectile = new G4QHadron(*pNucleon);// Copy selected PrNuc for String
223 curTargNucleus.StartLoop(); // Initialize newSelection forNucleons
224 tNucleon=curTargNucleus.GetNextNucleon(); // Select a nucleon
225 G4QHadron* aTarget = new G4QHadron(*tNucleon);// Copy selected nucleon for String
226 G4QInteraction* anInteraction = new G4QInteraction(aProjectile);
227 anInteraction->SetTarget(aTarget);
228 anInteraction->SetNumberOfDINRCollisions(1); // Consider this interaction as DINR
229 curInteractions.push_back(anInteraction); //--> now the Interaction is not empty
230 curProjNucleus.DoLorentzBoost(-theALVelocity);// Boost theResPrNucleus toRotAntiLab
231 curProjNucleus.SubtractNucleon(pNucleon); // Pointer to the used ProjNucleon
232 curProjNucleus.DoLorentzBoost(theALVelocity);// Boost theResProjNucleus back to CM
233 curTargNucleus.DoLorentzBoost(theCMVelocity);// Boost theResTgNucleus to RotatedLS
234 curTargNucleus.SubtractNucleon(tNucleon); // Pointer to the used TargNucleon
235 curTargNucleus.DoLorentzBoost(-theCMVelocity);// Boost theResTargNucleus back to CM
236#ifdef debug
237 G4cout<<"G4QIonIonCollision::Construct: DINR interaction is created"<<G4endl;
238#endif
239 break; // Break the WHILE of the pNucleons
240 }
241 // LOOP over nuclei of the target nucleus to select collisions
242 curTargNucleus.StartLoop(); // To get the same nucleon
243 G4int tnCnt = 0;
244#ifdef debu
245 G4cout<<"G4QIonIonCollision::Constr: Before WHILE over tNucleons, tA="<<tnA<<G4endl;
246#endif
247 while((tNucleon=curTargNucleus.GetNextNucleon()) && tnCnt<tnA && totalCuts<maxCuts)
248 {
249 ++tnCnt;
250 G4LorentzVector tNuc4M=tNucleon->Get4Momentum();
251#ifdef debug
252 G4cout<<"G4QIonIonCollision::Constr: OuterR="<<outerRadius<<", mC="<<maxCuts
253 <<", pA="<<curProjNucleus<<", tA="<<curTargNucleus<<G4endl;
254#endif
255 // Check the reaction threshold
256 G4double s_value = (tNuc4M + pNuc4M).mag2(); // Squared CM Energy of compound
257 G4double ThresholdMass = pNucleon->GetMass() + tNucleon->GetMass();
258#ifdef debug
259 G4cout<<"G4QInel::Constr: s="<<s_value<<", ThreshM="<<sqr(ThresholdMass)<<G4endl;
260#endif
261 ModelMode = SOFT; // NOT-Diffractive hadronization
262 if (s_value < 0.) // At ThP=0 is impossible(virtNucl)
263 {
264 G4cerr<<"*G4QInelast::Constr:s="<<s_value<<",pN4M="<<pNuc4M<<",tN4M="<<tNuc4M<<G4endl;
265 G4Exception("G4QIonIonCollision::Construct:","72",FatalException,"LowEn(NegS)");
266 }
267 if (s_value < sqr(ThresholdMass)) // --> Only diffractive interaction
268 {
269#ifdef debug
270 G4cout<<"G4QIonIonCollision::Constr: ***OnlyDiffraction***, ThM="<<ThresholdMass
271 <<">sqrt(s)="<<std::sqrt(s_value)<<" -> only Diffraction is possible"<<G4endl;
272#endif
273 ModelMode = DIFFRACTIVE;
274 }
275 G4ThreeVector tNucR=tNucleon->GetPosition(); // Position of the tNucleon WRTo tNucl
276 G4double dImpX = pImpX-tNucR.x();
277 G4double dImpY = pImpY-tNucR.y();
278 G4double Distance2=dImpX*dImpX+dImpY*dImpY;
279#ifdef sdebug
280 G4cout<<"G4QIonIonCollision::Construct: s="<<s_value<<", D2="<<Distance2<<G4endl;
281#endif
282 // Needs to be moved to Probability class @@
283 if(s_value<=10000.)
284 {
285 G4cout<<"-Warning-G4QIonIonCollision::Construct: s < .01 GeV^2, p4M="
286 <<pNucleon->Get4Momentum()<<",t4M="<<tNucleon->Get4Momentum()<<G4endl;
287 continue; // skip the rest of the targetNucleons
288 }
289 G4double Probability = theProbability.GetPomInelProbability(s_value, Distance2);// P_INEL
290 // test for inelastic collision
291#ifdef sdebug
292 G4cout<<"G4QIonIonCollision::Construct: Probubility="<<Probability<<G4endl;
293#endif
294 G4double rndNumber = G4UniformRand(); // For the printing purpose
295 // ModelMode = DIFFRACTIVE;
296#ifdef sdebug
297 G4cout<<"G4QIonIonCollision::Construct: NLOOP prob="<<Probability<<", rndm="
298 <<rndNumber<<", d="<<std::sqrt(Distance2)<<G4endl;
299#endif
300 if (Probability > rndNumber) // Inelastic (diffractive or soft) interaction (JOB)
301 {
302 G4QHadron* aProjectile = new G4QHadron(*pNucleon);// Copy selected pNuc to String
303 G4QHadron* aTarget = new G4QHadron(*tNucleon);// Copy selected tNucleon to String
304#ifdef edebug
305 G4cout<<"--->EMC-->G4QIonIonCollision::Construct: TargNucleon is filled, 4M/PDG="
306 <<aTarget->Get4Momentum()<<aTarget->GetPDGCode()<<G4endl;
307#endif
308 // Now the energy of the nucleons must be updated in CMS
309 curProjNucleus.DoLorentzBoost(-theALVelocity);// Boost theResNucleus toRotatedLS
310 curProjNucleus.SubtractNucleon(pNucleon); // Pointer to the used nucleon
311 curProjNucleus.DoLorentzBoost(theALVelocity); // Boost theResNucleus back to CM
312 curTargNucleus.DoLorentzBoost(theCMVelocity); // Boost theResNucleus toRotatedLS
313 curTargNucleus.SubtractNucleon(tNucleon); // Pointer to the used nucleon
314 curTargNucleus.DoLorentzBoost(-theCMVelocity);// Boost theResNucleus back to CM
315 if((theProbability.GetPomDiffProbability(s_value,Distance2)/Probability >
316 G4UniformRand() && ModelMode==SOFT ) || ModelMode==DIFFRACTIVE)
317 {
318 // ------------->> diffractive interaction @@ IsSingleDiffractive called once
319 if(IsSingleDiffractive()) ExciteSingDiffParticipants(aProjectile, aTarget);
320 else ExciteDiffParticipants(aProjectile, aTarget);
321 G4QInteraction* anInteraction = new G4QInteraction(aProjectile);
322 anInteraction->SetTarget(aTarget);
323 anInteraction->SetNumberOfDiffractiveCollisions(1); // Why not increment? M.K.?
324 curInteractions.push_back(anInteraction);//--> now theInteractions not empty
325 // @@ Why not breake the NLOOP, if only one diffractive can happend?
326 totalCuts++; // UpdateOfNucleons in't necessary
327#ifdef debug
328 G4cout<<"G4QIonIonCollision::Constr:NLOOP DiffInteract,tC="<<totalCuts<<G4endl;
329#endif
330 }
331 else
332 {
333 // -------------->> nondiffractive = soft interaction
334 // sample nCut+1 (cut Pomerons) pairs of strings can be produced
335 G4int nCut; // Result in a chosen number of cuts
336 G4double* running = new G4double[nCutMax];// @@ This limits the max cuts
337 for(nCut = 0; nCut < nCutMax; nCut++) // Calculates multiCut probabilities
338 {
339 running[nCut]= theProbability.GetCutPomProbability(s_value, Distance2, nCut+1);
340 if(nCut) running[nCut] += running[nCut-1];// Sum up with the previous one
341 }
342 G4double random = running[nCutMax-1]*G4UniformRand();
343 for(nCut = 0; nCut < nCutMax; nCut++) if(running[nCut] > random) break; // tNuc
344 delete [] running;
345#ifdef debug
346 G4cout<<"G4QIonIonCollision::Construct: NLOOP-Soft Chosen nCut="<<nCut<<G4endl;
347#endif
348 // @@ If nCut>0 interaction with a few nucleons is possible
349 // @@ nCut is found with big efforts and now nCut=0 ?? M.K. ?? !!
350 //nCut = 0; // @@ in original code ?? @@
351 aTarget->IncrementCollisionCount(nCut+1); // @@ What about multyNucleon target?
352 aProjectile->IncrementCollisionCount(nCut+1);
353 G4QInteraction* anInteraction = new G4QInteraction(aProjectile);
354 anInteraction->SetTarget(aTarget);
355 anInteraction->SetNumberOfSoftCollisions(nCut+1);
356 curInteractions.push_back(anInteraction); // Now curInteractions are not empty
357 totalCuts += nCut+1;
358#ifdef debug
359 G4cout<<"G4QIonIonCollision::Constr:NLOOP SoftInteract,tC="<<totalCuts<<G4endl;
360 impactUsed=Distance2;
361#endif
362 }
363 }// Probability selection
364 } // End of While over target nucleons
365 G4int nCurInt=curInteractions.size();
366 for(G4int ni=0; ni < nCurInt; ++ni) theInteractions.push_back(curInteractions[ni]);
367 curInteractions.clear(); // Probably, not necessary...
368 } // End of WHILE over projectile nucleons
369 } // End of WHILE over attempts to find at least one nucleus-nucleus interaction
370 theProjNucleus.DeleteNucleons();
371 theTargNucleus.DeleteNucleons();
372 theProjNucleus = curProjNucleus;
373 theTargNucleus = curTargNucleus;
374 curProjNucleus.DeleteNucleons();
375 curTargNucleus.DeleteNucleons();
376 G4int nInt=theInteractions.size();
377#ifdef debug
378 G4cout<<"G4QIonIonCollision::Constr: #ofInteractions = "<<nInt<<", #OfDINR = "
379 <<theInteractions[0]->GetNumberOfDINRCollisions()<<G4endl;
380#endif
381 if(!nInt || (nInt==1 && theInteractions[0]->GetNumberOfDINRCollisions()==1)) // QFreeInel
382 {
383 G4QHadron* aTarget;
384 G4QHadron* aProjectile;
385 if(nInt) // Take Targ/Proj from the Interaction
386 {
387 aTarget=theInteractions[0]->GetTarget();
388 aProjectile=theInteractions[0]->GetProjectile();
389 delete theInteractions[0];
390 theInteractions.clear();
391 }
392 else // Create a new target nucleon (?)
393 {
394 theProjNucleus.StartLoop(); // To get the same nucleon from projec
395 pNucleon=theProjNucleus.GetNextNucleon(); // Get theNucleon to create projectNuc
396 aProjectile = new G4QHadron(*pNucleon); // Copy selected pNucleon for String
397 theProjNucleus.DoLorentzBoost(-theALVelocity); // Boost theResProjNucleus toRotatedLS
398 theProjNucleus.SubtractNucleon(pNucleon); // Pointer to SelProjNucleon to delete
399 theProjNucleus.DoLorentzBoost(theALVelocity); // Boost theResProNucleus back to CMS
400 theTargNucleus.StartLoop(); // To get the same nucleon from target
401 tNucleon=theTargNucleus.GetNextNucleon(); // Get theNucleon to create targetNucl
402 aTarget = new G4QHadron(*tNucleon); // Copy selected tNucleon for String
403 theTargNucleus.DoLorentzBoost(theCMVelocity); // Boost theResTargNucleus toRotatedLS
404 theTargNucleus.SubtractNucleon(tNucleon); // Pointer to SelTargNucleon to delete
405 theTargNucleus.DoLorentzBoost(-theCMVelocity); // Boost theResTargNucleus back to CMS
406 }
407 G4QContent QQC=aTarget->GetQC()+aProjectile->GetQC(); // QContent of the compound
408 G4LorentzVector Q4M=aTarget->Get4Momentum()+aProjectile->Get4Momentum(); // 4-mom of Q
409 delete aTarget;
410 delete aProjectile;
411 // @@ 4-Mom should be converted to LS for the TargQuasmon and to AL for the ProjQuasmon
412 Q4M.boost(theCMVelocity);
413 Q4M=toLab*Q4M;
414 G4Quasmon* stringQuasmon = new G4Quasmon(QQC, Q4M);
415 theQuasmons.push_back(stringQuasmon);
416 theTargNucleus.DoLorentzBoost(theCMVelocity); // BoostTheResidualNucleus toRotatedLS
417 theTargNucleus.DoLorentzRotation(toLab); // Recove Z-direction in LS ("LS"->LS)
418 return;
419 }
420 //
421 // ------------------ now build the parton pairs for the strings ------------------
422 //
423 for(G4int i=0; i<nInt; i++)
424 {
425 theInteractions[i]->SplitHadrons();
426#ifdef edebug
427 G4QHadron* projH=theInteractions[i]->GetProjectile(); // Projectile of theInteraction
428 G4QHadron* targH=theInteractions[i]->GetTarget(); // Target of the Interaction
429 G4LorentzVector pSP(0.,0.,0.,0.); // Sum of parton's 4mom's for proj
430 G4LorentzVector tSP(0.,0.,0.,0.); // Sum of parton's 4mom's for proj
431 std::list<G4QParton*> projCP=projH->GetColor(); // Pointers to proj Color-partons
432 std::list<G4QParton*> projAC=projH->GetAntiColor();// PointersTo projAntiColorPartons
433 std::list<G4QParton*> targCP=targH->GetColor(); // Pointers to targ Color-partons
434 std::list<G4QParton*> targAC=targH->GetAntiColor();// PointersTo targAntiColorPartons
435 std::list<G4QParton*>::iterator picp = projCP.begin();
436 std::list<G4QParton*>::iterator pecp = projCP.end();
437 std::list<G4QParton*>::iterator piac = projAC.begin();
438 std::list<G4QParton*>::iterator peac = projAC.end();
439 std::list<G4QParton*>::iterator ticp = targCP.begin();
440 std::list<G4QParton*>::iterator tecp = targCP.end();
441 std::list<G4QParton*>::iterator tiac = targAC.begin();
442 std::list<G4QParton*>::iterator teac = targAC.end();
443 for(; picp!=pecp&& piac!=peac&& ticp!=tecp&& tiac!=teac; ++picp,++piac,++ticp,++tiac)
444 {
445 pSP+=(*picp)->Get4Momentum();
446 pSP+=(*piac)->Get4Momentum();
447 tSP+=(*ticp)->Get4Momentum();
448 tSP+=(*tiac)->Get4Momentum();
449 }
450 G4cout<<"-EMC-G4QIonIonCollision::Construct: Interaction#"<<i<<",dP4M="
451 <<projH->Get4Momentum()-pSP<<",dT4M="<<targH->Get4Momentum()-tSP<<G4endl;
452#endif
453 }
454 //
455 // >>........>> make soft collisions (ordering is vital)
456 //
457 G4QInteractionVector::iterator it;
458#ifdef debug
459 G4cout<<"G4QIonIonCollision::Constr: Creation ofSoftCollisionPartonPair STARTS"<<G4endl;
460#endif
461 G4bool rep=true;
462 while(rep && theInteractions.size())
463 {
464 for(it = theInteractions.begin(); it != theInteractions.end(); ++it)
465 {
466 G4QInteraction* anIniteraction = *it;
467 G4QPartonPair* aPair=0;
468 G4int nSoftCollisions = anIniteraction->GetNumberOfSoftCollisions();
469#ifdef debug
470 G4cout<<"G4QIonIonCollision::Construct: #0f SoftCollisions ="<<nSoftCollisions<<G4endl;
471#endif
472 if (nSoftCollisions)
473 {
474 G4QHadron* pProjectile = anIniteraction->GetProjectile();
475 G4QHadron* pTarget = anIniteraction->GetTarget();
476 for (G4int j = 0; j < nSoftCollisions; j++)
477 {
478 aPair = new G4QPartonPair(pTarget->GetNextParton(),
479 pProjectile->GetNextAntiParton(),
480 G4QPartonPair::SOFT, G4QPartonPair::TARGET);
481 thePartonPairs.push_back(aPair); // A target pair (Why TAGRET?)
482 aPair = new G4QPartonPair(pProjectile->GetNextParton(),
483 pTarget->GetNextAntiParton(),
484 G4QPartonPair::SOFT, G4QPartonPair::PROJECTILE);
485 thePartonPairs.push_back(aPair); // A projectile pair (Why Projectile?)
486#ifdef debug
487 G4cout<<"--->G4QIonIonCollision::Constr: SOFT, 2 parton pairs are filled"<<G4endl;
488#endif
489 }
490 delete *it;
491 it=theInteractions.erase(it); // SoftInteractions're converted&erased
492 if( it != theInteractions.begin() )// To avoid going below begin() (for Windows)
493 {
494 it--;
495 rep=false;
496 }
497 else
498 {
499 rep=true;
500 break;
501 }
502 }
503 else rep=false;
504 }
505 }
506#ifdef debug
507 G4cout<<"G4QIonIonCollision::Constr: -> Parton pairs for SOFT strings are made"<<G4endl;
508#endif
509 //
510 // >>.............>> make the rest as the diffractive interactions
511 //
512 for(unsigned i = 0; i < theInteractions.size(); i++) // Interactions are reduced bySoft
513 {
514 // The double or single diffraction is defined by the presonce of proj/targ partons
515 G4QInteraction* anIniteraction = theInteractions[i];
516 G4QPartonPair* aPartonPair;
517#ifdef debug
518 G4cout<<"G4QIonIonCollision::Constr: CreationOfDiffractivePartonPairs, i="<<i<<G4endl;
519#endif
520 // the projectile diffraction parton pair is created first
521 G4QHadron* aProjectile = anIniteraction->GetProjectile();
522 G4QParton* aParton = aProjectile->GetNextParton();
523 if (aParton)
524 {
525 aPartonPair = new G4QPartonPair(aParton, aProjectile->GetNextAntiParton(),
526 G4QPartonPair::DIFFRACTIVE,
527 G4QPartonPair::PROJECTILE);
528 thePartonPairs.push_back(aPartonPair);
529#ifdef debug
530 G4cout<<"G4QIonIonCollision::Constr: proj Diffractive PartonPair is filled"<<G4endl;
531#endif
532 }
533 // then the target diffraction parton pair is created
534 G4QHadron* aTarget = anIniteraction->GetTarget();
535 aParton = aTarget->GetNextParton();
536 if (aParton)
537 {
538 aPartonPair = new G4QPartonPair(aParton, aTarget->GetNextAntiParton(),
539 G4QPartonPair::DIFFRACTIVE, G4QPartonPair::TARGET);
540 thePartonPairs.push_back(aPartonPair);
541#ifdef debug
542 G4cout<<"G4QIonIonCollision::Constr: target Diffractive PartonPair's filled"<<G4endl;
543#endif
544 }
545 }
546#ifdef debug
547 G4cout<<"G4QIonIonCollision::Construct: DiffractivePartonPairs are created"<<G4endl;
548#endif
549 //
550 // >>.......>> clean-up Interactions, if necessary
551 //
552 std::for_each(theInteractions.begin(),theInteractions.end(), DeleteQInteraction());
553 theInteractions.clear();
554#ifdef debug
555 G4cout<<"G4QIonIonCollision::Construct: Temporary objects are cleaned up"<<G4endl;
556#endif
557 // This function prepares theBoost for transformation of secondaries to LS (-ProjRot!)
558 theProjNucleus.DoLorentzBoost(theCMVelocity);// Boost theResidualProjNucleus to RotatedLS
559 theTargNucleus.DoLorentzBoost(theCMVelocity);// Boost theResidualTargNucleus to RotatedLS
560 // @@ Nucleus isn't completely in LS, it needs the toZ (-ProjRot) rotation to consE/M
561#ifdef debug
562 G4cout<<">>........>>G4QIonIonCollision::Construct: >>..>> Strings are created "<<G4endl;
563#endif
564 G4QPartonPair* aPair;
565 G4QString* aString=0;
566 while(thePartonPairs.size()) // @@ At present noDifference in stringBuild (? M.K.)
567 {
568 aPair = thePartonPairs.back(); // Get the parton pair
569 thePartonPairs.pop_back(); // Clean up thePartonPairPointer in the Vector
570#ifdef debug
571 G4cout<<"G4QIonIonCollision::Construct:StringType="<<aPair->GetCollisionType()<<G4endl;
572#endif
573 aString= new G4QString(aPair);
574#ifdef debug
575 G4cout<<"G4QIonIonCollision::Construct:NewString4M="<<aString->Get4Momentum()<<G4endl;
576#endif
577 aString->Boost(theCMVelocity); // ! Strings are moved to ZLS when pushed !
578 strings.push_back(aString);
579 stringsInitted=true;
580 delete aPair;
581 } // End of the String Creation LOOP
582#ifdef edebug
583 G4LorentzVector sum=theProjNucleus.Get4Momentum()+theTargNucleus.Get4Momentum(); // in LS
584 G4int rChg=totChg-theProjNucleus.GetZ()-theTargNucleus.GetZ();
585 G4int rBaN=totBaN-theProjNucleus.GetA()-theTargNucleus.GetA();
586 G4int nStrings=strings.size();
587 G4cout<<"-EMC-G4QIonIonCollision::Constr:#ofString="<<nStrings<<",tNuc4M="<<sum<<G4endl;
588 for(G4int i=0; i<nStrings; i++)
589 {
590 G4QString* prString=strings[i];
591 G4LorentzVector strI4M=prString->Get4Momentum();
592 sum+=strI4M;
593 G4int sChg=prString->GetCharge();
594 G4int sBaN=prString->GetBaryonNumber();
595 G4int LPDG=prString->GetLeftParton()->GetPDGCode();
596 G4int RPDG=prString->GetRightParton()->GetPDGCode();
597 G4QContent LQC =prString->GetLeftParton()->GetQC();
598 G4QContent RQC =prString->GetRightParton()->GetQC();
599 rChg-=sChg;
600 rBaN-=sBaN;
601 G4cout<<"-EMC-G4QIonIonCollision::Construct: String#"<<i<<", 4M="<<strI4M<<", LPDG="
602 <<LPDG<<LQC<<",RPDG="<<RPDG<<RQC<<", Ch="<<sChg<<", BN="<<sBaN<<G4endl;
603 }
604 G4cout<<"-EMC-G4QInel::Constr: r4M="<<sum-totZLS4M<<", rC="<<rChg<<", rB="<<rBaN<<G4endl;
605#endif
606 if(!stringsInitted)
607 {
608 G4cerr<<"*****G4QIonIonCollision::Construct:**** No strings are created *****"<<G4endl;
609 G4Exception("G4QIonIonCollision::Construct:","72",FatalException,"No Strings created");
610 }
611#ifdef debug
612 G4cout<<"G4QIonIonCollision::Constr:BeforeRotation, #0fStrings="<<strings.size()<<G4endl;
613#endif
614 //
615 // ---------------- At this point the strings must be already created in "LS" -----------
616 //
617 for(unsigned astring=0; astring < strings.size(); astring++)
618 strings[astring]->LorentzRotate(toLab); // Recove Z-direction in LS ("LS"->LS)
619 theProjNucleus.DoLorentzRotation(toLab); // Recove Z-dir in LS ("LS"->LS) for ProjNucleus
620 theTargNucleus.DoLorentzRotation(toLab); // Recove Z-dir in LS ("LS"->LS) for TargNucleus
621 // Now everything is in LS system
622#ifdef edebug
623 G4LorentzVector sm=theProjNucleus.Get4Momentum()+theTargNucleus.Get4Momentum();// NucInLS
624 G4int rCg=totChg-theProjNucleus.GetZ()-theTargNucleus.GetZ();
625 G4int rBC=totBaN-theProjNucleus.GetA()-theTargNucleus.GetA();
626 G4int nStrs=strings.size();
627 G4cout<<"-EMCLS-G4QIonIonCollision::Constr:#ofS="<<nStrings<<",Nuc4M(E=M)="<<sum<<G4endl;
628 for(G4int i=0; i<nStrs; i++)
629 {
630 G4LorentzVector strI4M=strings[i]->Get4Momentum();
631 sm+=strI4M;
632 G4int sChg=strings[i]->GetCharge();
633 rCg-=sChg;
634 G4int sBaN=strings[i]->GetBaryonNumber();
635 rBC-=sBaN;
636 G4cout<<"-EMCLS-G4QInel::Construct:String#"<<i<<",4M="<<strI4M<<strI4M.m()<<",Charge="
637 <<sChg<<",BaryN="<<sBaN<<G4endl;
638 }
639 G4cout<<"-EMCLS-...G4QInel::Constr: r4M="<<sm-totLS4M<<", rC="<<rCg<<",rB="<<rBC<<G4endl;
640#endif
641 //
642 // --- Strings are created, but we should try to get rid of negative mass strings -----
643 //
644 SwapPartons();
645 //
646 // --- Strings are created, but we should get rid of too light strings (Mmin+MPi0) -----
647 //
648 G4int problem=0; // 0="no problem", incremented by ASIS
649 G4QStringVector::iterator ist;
650 G4bool con=true;
651 while(con && strings.size())
652 {
653 for(ist = strings.begin(); ist < strings.end(); ++ist)
654 {
655 G4bool bad=true;
656 G4LorentzVector cS4M=(*ist)->Get4Momentum();
657 G4double cSM2=cS4M.m2(); // Squared mass of the String
658 G4QParton* cLeft=(*ist)->GetLeftParton();
659 G4QParton* cRight=(*ist)->GetRightParton();
660 G4int cLT=cLeft->GetType();
661 G4int cRT=cRight->GetType();
662 G4int cLPDG=cLeft->GetPDGCode();
663 G4int cRPDG=cRight->GetPDGCode();
664 G4int aLPDG=0;
665 G4int aRPDG=0;
666 if (cLPDG > 7) aLPDG= cLPDG/100;
667 else if(cLPDG <-7) aLPDG=(-cLPDG)/100;
668 if (cRPDG > 7) aRPDG= cRPDG/100;
669 else if(cRPDG <-7) aRPDG=(-cRPDG)/100;
670 G4int L1=0;
671 G4int L2=0;
672 if(aLPDG)
673 {
674 L1=aLPDG/10;
675 L2=aLPDG%10;
676 }
677 G4int R1=0;
678 G4int R2=0;
679 if(aRPDG)
680 {
681 R1=aRPDG/10;
682 R2=aRPDG%10;
683 }
684 G4double cSM=cSM2;
685 if(cSM2>0.) cSM=std::sqrt(cSM2);
686#ifdef debug
687 G4cout<<"G4QIonIonCollision::Construct: Needs Fusion? cLPDG="<<cLPDG<<",cRPDG="<<cRPDG
688 <<",cM(cM2 if neg)="<<cSM<<G4endl;
689#endif
690 if(cSM>0.) // Mass can be calculated
691 {
692 G4bool single=true;
693 G4double miM=0.; // Proto of the Min HadronString Mass
694 if(cLT==2 && cRT==2)
695 {
696 if(L1!=R1 && L1!=R2 && L2!=R1 && L2!=R2) // Unreducable DiQ-aDiQ
697 {
698 single=false;
699 G4QPDGCode tmp;
700 std::pair<G4int,G4int> paB=tmp.MakeTwoBaryons(L1, L2, R1, R2);
701 miM=G4QPDGCode(paB.first).GetMass()+G4QPDGCode(paB.second).GetMass();
702 }
703 }
704 if(single) miM=G4QPDGCode((*ist)->GetQC().GetSPDGCode()).GetMass() + mPi0;//MinHSMass
705 //if(single) miM=G4QPDGCode((*ist)->GetQC().GetSPDGCode()).GetMass();//MinHSMass
706#ifdef debug
707 G4cout<<"G4QInel::Const:*IsItGood? realM="<<std::sqrt(cSM2)<<" > GSM="<<miM<<G4endl;
708#endif
709 if(std::sqrt(cSM2) > miM) bad=false; // String is OK
710 }
711 if(bad) // String should be merged with others
712 {
713#ifdef debug
714 G4cout<<"G4QInel::Const:TryFuse,L1="<<L1<<",L2="<<L2<<",R1="<<R1<<",R2="<<R2<<G4endl;
715#endif
716 G4int cST=cLT+cRT;
717 G4double excess=-DBL_MAX; // The value to be maximized excess M
718 G4double maxiM2=-DBL_MAX; // The value to be maximized M2
719 G4QStringVector::iterator sst; // Selected partner string
720 G4QStringVector::iterator pst;
721 G4int sLPDG=0; // selectedLeft (like inStringPartner)
722 G4int sRPDG=0; // selectedRight(like inStringPartner)
723 G4int sOrd=0; // selected Order of PartonFusion
724 G4bool minC=true; // for the case when M2<0
725 if(cSM2>0.) minC=false; // If M2>0 already don'tSearchFor M2>0
726 for(pst = strings.begin(); pst < strings.end(); pst++) if(pst != ist)
727 {
728 G4LorentzVector pS4M=(*pst)->Get4Momentum()+cS4M; // Summed 4-momentum
729 G4int nLPDG=0; // new Left (like in theStringPartner)
730 G4int nRPDG=0; // new Right(like in theStringPartner)
731 G4double pExcess=-DBL_MAX; // Prototype of the excess
732 G4double pSM2=pS4M.m2(); // Squared mass of the Fused Strings
733#ifdef debug
734 G4cout<<"->G4QIonIonCollision::Construct: sum4M="<<pS4M<<", M2="<<pSM2<<G4endl;
735#endif
736 //if(pSM2>0.) // The partner can be a candidate
737 //{
738 G4QParton* pLeft=(*pst)->GetLeftParton();
739 G4QParton* pRight=(*pst)->GetRightParton();
740 G4int pLT=pLeft->GetType();
741 G4int pRT=pRight->GetType();
742 G4int pLPDG=pLeft->GetPDGCode();
743 G4int pRPDG=pRight->GetPDGCode();
744 G4int LPDG=0;
745 G4int RPDG=0;
746 if (pLPDG > 7) LPDG= pLPDG/100;
747 else if(pLPDG <-7) LPDG=(-pLPDG)/100;
748 if (pRPDG > 7) RPDG= pRPDG/100;
749 else if(pRPDG <-7) RPDG=(-pRPDG)/100;
750 G4int pL1=0;
751 G4int pL2=0;
752 if(LPDG)
753 {
754 pL1=LPDG/10;
755 pL2=LPDG%10;
756 }
757 G4int pR1=0;
758 G4int pR2=0;
759 if(RPDG)
760 {
761 pR1=RPDG/10;
762 pR2=RPDG%10;
763 }
764 G4int pST=pLT+pRT;
765#ifdef debug
766 G4cout<<"G4QInel::Construct: Partner/w pLPDG="<<pLPDG<<", pRPDG="<<pRPDG<<", pM2="
767 <<pSM2<<G4endl;
768#endif
769 // Different fromCompactAlrorithm ofStringFusionAfterDecay (no DiQaDiQ reduction)
770 G4int tf=0; // Type combination flag
771 G4int af=0; // Annihilatio combination flag
772 if (cST==2 && pST==2) // QaQ + QaQ = DiQaDiQ (always)
773 {
774 tf=1;
775 af=1;
776 }
777 else if(cST==2 && pST==3) // QaQ + QDiQ/aQaDiQ = QDiQ/aQaDiQ (s)
778 {
779 tf=2;
780 if (pLPDG > 7 &&
781 ( (cLPDG<0 && (-cLPDG==pL1 || -cLPDG==pL2 || -cLPDG==pRPDG) ) ||
782 (cRPDG<0 && (-cRPDG==pL1 || -cRPDG==pL2 || -cRPDG==pRPDG) )
783 )
784 ) af=1;
785 else if(pRPDG > 7 &&
786 ( (cLPDG<0 && (-cLPDG==pR1 || -cLPDG==pR2 || -cLPDG==pLPDG) ) ||
787 (cRPDG<0 && (-cRPDG==pR1 || -cRPDG==pR2 || -cRPDG==pLPDG) )
788 )
789 ) af=2;
790 else if(pLPDG <-7 &&
791 ( (cLPDG>0 && ( cLPDG==pL1 || cLPDG==pL2 || cLPDG==-pRPDG) ) ||
792 (cRPDG>0 && ( cRPDG==pL1 || cRPDG==pL2 || cRPDG==-pRPDG) )
793 )
794 ) af=3;
795 else if(pRPDG <-7 &&
796 ( (cLPDG>0 && ( cLPDG==pR1 || cLPDG==pR2 || cLPDG==-pLPDG) ) ||
797 (cRPDG>0 && ( cRPDG==pR1 || cRPDG==pR2 || cRPDG==-pLPDG) )
798 )
799 ) af=4;
800#ifdef debug
801 else G4cout<<"G4QIonIonCollision::Constr:2(QaQ+QDiQ/aQaDiQ) Can't fuse"<<G4endl;
802#endif
803 }
804 else if(cST==3 && pST==2) // QDiQ/aQaDiQ + QaQ = QDiQ/aQaDiQ (s)
805 {
806 tf=3;
807 if (cLPDG > 7 &&
808 ( (pLPDG<0 && (-pLPDG==L1 || -pLPDG==L2) ) ||
809 (pRPDG<0 && (-pRPDG==L1 || -pRPDG==L2) )
810 )
811 ) af=1;
812 else if(cRPDG > 7 &&
813 ( (pLPDG<0 && (-pLPDG==R1 || -pLPDG==R2) ) ||
814 (pRPDG<0 && (-pRPDG==R1 || -pRPDG==R2) )
815 )
816 ) af=2;
817 else if(cLPDG <-7 &&
818 ( (pLPDG>0 && ( pLPDG==L1 || pLPDG==L2) ) ||
819 (pRPDG>0 && ( pRPDG==L1 || pRPDG==L2) )
820 )
821 ) af=3;
822 else if(cRPDG <-7 &&
823 ( (pLPDG>0 && ( pLPDG==R1 || pLPDG==R2) ) ||
824 (pRPDG>0 && ( pRPDG==R1 || pRPDG==R2) )
825 )
826 ) af=4;
827#ifdef debug
828 else G4cout<<"G4QIonIonCollision::Constr:3(QDiQ/aQaDiQ+QaQ) Can't fuse"<<G4endl;
829#endif
830 }
831 else if(cST==2 && pST==4) // QaQ + aDiQDiQ = QaQ (double)
832 {
833 tf=4;
834 if (pLPDG > 7) // pRPDG <-7
835 {
836 if ( (-cLPDG==pL1 || -cLPDG==pL2) && (cRPDG==pR1 || cRPDG==pR2) ) af=1;
837 else if( (-cRPDG==pL1 || -cRPDG==pL2) && (cLPDG==pR1 || cLPDG==pR2) ) af=2;
838 }
839 else if(pRPDG > 7) // pLPDG <-7
840 {
841 if ( (-cRPDG==pR1 || -cRPDG==pR2) && (cLPDG==pL1 || cLPDG==pL2) ) af=3;
842 else if( (-cLPDG==pR1 || -cLPDG==pR2) && (cRPDG==pL1 || cRPDG==pL2) ) af=4;
843 }
844#ifdef debug
845 else G4cout<<"-G4QIonIonCollision::Constr: 4 (QaQ+aQDiQDiQ) Can't fuse"<<G4endl;
846#endif
847 }
848 else if(cST==4 && pST==2) // aDiQDiQ + QaQ = QaQ (double)
849 {
850 tf=5;
851 if (cLPDG > 7) // cRPDG<-7
852 {
853 if ( (-pLPDG==L1 || -pLPDG==L2) && (pRPDG==R1 || pRPDG==R2) ) af=1;
854 else if( (-pRPDG==L1 || -pRPDG==L2) && (pLPDG==R1 || pLPDG==R2) ) af=2;
855 }
856 else if(cRPDG > 7) // cLPDG<-7
857 {
858 if ( (-pRPDG==R1 || -pRPDG==R2) && (pLPDG==L1 || pLPDG==L2) ) af=3;
859 else if( (-pLPDG==R1 || -pLPDG==R2) && (pRPDG==L1 || pRPDG==L2) ) af=4;
860 }
861#ifdef debug
862 else G4cout<<"-G4QIonIonCollision::Constr: 5 (aQDiQDiQ+QaQ) Can't fuse"<<G4endl;
863#endif
864 }
865 else if(cST==3 && pST==3) // QDiQ + aQaDiQ = QaQ (double)
866 {
867 tf=6;
868 if(pLPDG > 7)
869 {
870 if (cLPDG<-7 && (-cRPDG==pL1 || -cRPDG==pL2) && (pRPDG==L1 || pRPDG==L2))
871 af=1;
872 else if(cRPDG<-7 && (-cLPDG==pL1 || -cLPDG==pL2) && (pRPDG==R1 || pRPDG==R2))
873 af=2;
874 }
875 else if(pRPDG > 7)
876 {
877 if (cLPDG<-7 && (-cRPDG==pR1 || -cRPDG==pR2) && (pLPDG==L1 || pLPDG==L2))
878 af=3;
879 else if(cRPDG<-7 && (-cLPDG==pR1 || -cLPDG==pR2) && (pLPDG==R1 || pLPDG==R2))
880 af=4;
881 }
882 else if(cLPDG > 7)
883 {
884 if (pLPDG<-7 && (-pRPDG==L1 || -pRPDG==L2) && (cRPDG==pL1 || cRPDG==pL2))
885 af=5;
886 else if(pRPDG<-7 && (-pLPDG==L1 || -pLPDG==L2) && (cRPDG==pR1 || cRPDG==pR2))
887 af=6;
888 }
889 else if(cRPDG > 7)
890 {
891 if (pLPDG<-7 && (-pRPDG==R1 || -pRPDG==R2) && (cLPDG==pL1 || cLPDG==pL2))
892 af=7;
893 else if(pRPDG<-7 && (-pLPDG==R1 || -pLPDG==R2) && (cLPDG==pR1 || cLPDG==pR2))
894 af=8;
895 }
896#ifdef debug
897 else G4cout<<"-G4QIonIonCollision::Constr: 6 (QDiQ+aQaDiQ) Can't fuse"<<G4endl;
898#endif
899 }
900#ifdef debug
901 G4cout<<"G4QIonIonCollision::Constr: **Possibility**, tf="<<tf<<",af="<<af<<G4endl;
902#endif
903 if(tf && af)
904 {
905 // Strings can be fused, update the max excess and remember usefull switches
906 G4int order=0; // LL/RR (1) or LR/RL (2) PartonFusion
907 switch (tf)
908 {
909 case 1: // ------------------------------------> QaQ + QaQ = DiQaDiQ (always)
910 if (cLPDG > 0 && pLPDG > 0)
911 {
912 order= 1;
913 if (cLPDG > pLPDG) nLPDG=cLPDG*1000+pLPDG*100+1;
914 else if(cLPDG < pLPDG) nLPDG=pLPDG*1000+cLPDG*100+1;
915 else nLPDG=pLPDG*1000+cLPDG*100+3;
916 if (cRPDG < pRPDG) nRPDG=cRPDG*1000+pRPDG*100-1;
917 else if(cRPDG > pRPDG) nRPDG=pRPDG*1000+cRPDG*100-1;
918 else nRPDG=pRPDG*1000+cRPDG*100-3;
919 }
920 else if(cLPDG < 0 && pLPDG < 0)
921 {
922 order= 1;
923 if (cRPDG > pRPDG) nRPDG=cRPDG*1000+pRPDG*100+1;
924 else if(cRPDG < pRPDG) nRPDG=pRPDG*1000+cRPDG*100+1;
925 else nRPDG=pRPDG*1000+cRPDG*100+3;
926 if (cLPDG < pLPDG) nLPDG=cLPDG*1000+pLPDG*100-1;
927 else if(cLPDG > pLPDG) nLPDG=pLPDG*1000+cLPDG*100-1;
928 else nLPDG=pLPDG*1000+cLPDG*100-3;
929 }
930 else if(cRPDG > 0 && pLPDG > 0)
931 {
932 order=-1;
933 if (cRPDG > pLPDG) nLPDG=cRPDG*1000+pLPDG*100+1;
934 else if(cRPDG < pLPDG) nLPDG=pLPDG*1000+cRPDG*100+1;
935 else nLPDG=pLPDG*1000+cRPDG*100+3;
936 if (cLPDG < pRPDG) nRPDG=cLPDG*1000+pRPDG*100-1;
937 else if(cLPDG > pRPDG) nRPDG=pRPDG*1000+cLPDG*100-1;
938 else nRPDG=pRPDG*1000+cLPDG*100-3;
939 }
940 else if(cRPDG < 0 && pLPDG < 0)
941 {
942 order=-1;
943 if (cLPDG > pRPDG) nRPDG=cLPDG*1000+pRPDG*100+1;
944 else if(cLPDG < pRPDG) nRPDG=pRPDG*1000+cLPDG*100+1;
945 else nRPDG=pRPDG*1000+cLPDG*100+3;
946 if (cRPDG < pLPDG) nLPDG=cRPDG*1000+pLPDG*100-1;
947 else if(cRPDG > pLPDG) nLPDG=pLPDG*1000+cRPDG*100-1;
948 else nLPDG=pLPDG*1000+cRPDG*100-3;
949 }
950 break;
951 case 2: // ------------------------> QaQ + QDiQ/aQaDiQ = QDiQ/aQaDiQ (single)
952 switch (af)
953 {
954 case 1: // ....................... pLPDG > 7
955 if(cLPDG < 0)
956 {
957 order= 1;
958 if(-cLPDG==pRPDG)
959 {
960 nLPDG=pLPDG;
961 nRPDG=cRPDG;
962 }
963 else
964 {
965 if (cRPDG > pRPDG) nRPDG=cRPDG*1000+pRPDG*100+1;
966 else if(cRPDG < pRPDG) nRPDG=pRPDG*1000+cRPDG*100+1;
967 else nRPDG=pRPDG*1000+cRPDG*100+3;
968 if (-cLPDG == pL1) nLPDG=pL2;
969 else nLPDG=pL1; // -cLPDG == pL2
970 }
971 }
972 else // cRPDG < 0
973 {
974 order=-1;
975 if(-cRPDG==pRPDG)
976 {
977 nLPDG=pLPDG;
978 nRPDG=cLPDG;
979 }
980 else
981 {
982 if (cLPDG > pRPDG) nRPDG=cLPDG*1000+pRPDG*100+1;
983 else if(cLPDG < pRPDG) nRPDG=pRPDG*1000+cLPDG*100+1;
984 else nRPDG=pRPDG*1000+cLPDG*100+3;
985 if (-cRPDG == pL1) nLPDG=pL2;
986 else nLPDG=pL1; // -cRPDG == pL2
987 }
988 }
989 break;
990 case 2: // ....................... pRPDG > 7
991 if(cLPDG < 0)
992 {
993 order=-1;
994 if(-cLPDG==pLPDG)
995 {
996 nLPDG=cRPDG;
997 nRPDG=pRPDG;
998 }
999 else
1000 {
1001 if (cRPDG > pLPDG) nLPDG=cRPDG*1000+pLPDG*100+1;
1002 else if(cRPDG < pLPDG) nLPDG=pLPDG*1000+cRPDG*100+1;
1003 else nLPDG=pLPDG*1000+cRPDG*100+3;
1004 if (-cLPDG == pR1) nRPDG=pR2;
1005 else nRPDG=pR1; // -cLPDG == pR2
1006 }
1007 }
1008 else // cRPDG < 0
1009 {
1010 order= 1;
1011 if(-cRPDG==pLPDG)
1012 {
1013 nLPDG=cLPDG;
1014 nRPDG=pRPDG;
1015 }
1016 else
1017 {
1018 if (cLPDG > pLPDG) nLPDG=cLPDG*1000+pLPDG*100+1;
1019 else if(cLPDG < pLPDG) nLPDG=pLPDG*1000+cLPDG*100+1;
1020 else nLPDG=pLPDG*1000+cLPDG*100+3;
1021 if (-cRPDG == pR1) nRPDG=pR2;
1022 else nRPDG=pR1; // -cRPDG == pR2
1023 }
1024 }
1025 break;
1026 case 3: // ....................... pLPDG <-7
1027 if(cLPDG > 0)
1028 {
1029 order= 1;
1030 if(cLPDG==-pRPDG)
1031 {
1032 nLPDG=pLPDG;
1033 nRPDG=cRPDG;
1034 }
1035 else
1036 {
1037 if (cRPDG < pRPDG) nRPDG=cRPDG*1000+pRPDG*100-1;
1038 else if(cRPDG > pRPDG) nRPDG=pRPDG*1000+cRPDG*100-1;
1039 else nRPDG=pRPDG*1000+cRPDG*100-3;
1040 if ( cLPDG == pL1) nLPDG=-pL2;
1041 else nLPDG=-pL1; // cLPDG == pL2
1042 }
1043 }
1044 else // cRPDG > 0
1045 {
1046 order=-1;
1047 if(cRPDG==-pRPDG)
1048 {
1049 nLPDG=pLPDG;
1050 nRPDG=cLPDG;
1051 }
1052 else
1053 {
1054 if (cLPDG < pRPDG) nRPDG=cLPDG*1000+pRPDG*100-1;
1055 else if(cLPDG > pRPDG) nRPDG=pRPDG*1000+cLPDG*100-1;
1056 else nRPDG=pRPDG*1000+cLPDG*100-3;
1057 if ( cRPDG == pL1) nLPDG=-pL2;
1058 else nLPDG=-pL1; // cRPDG == pL2
1059 }
1060 }
1061 break;
1062 case 4: // ....................... pRPDG <-7
1063 if(cLPDG > 0)
1064 {
1065 order=-1;
1066 if(cLPDG==-pLPDG)
1067 {
1068 nLPDG=cRPDG;
1069 nRPDG=pRPDG;
1070 }
1071 else
1072 {
1073 if (cRPDG < pLPDG) nLPDG=cRPDG*1000+pLPDG*100-1;
1074 else if(cRPDG > pLPDG) nLPDG=pLPDG*1000+cRPDG*100-1;
1075 else nLPDG=pLPDG*1000+cRPDG*100-3;
1076 if ( cLPDG == pR1) nRPDG=-pR2;
1077 else nRPDG=-pR1; // cLPDG == pR2
1078 }
1079 }
1080 else // cRPDG > 0
1081 {
1082 order= 1;
1083 if(cRPDG==-pLPDG)
1084 {
1085 nLPDG=cLPDG;
1086 nRPDG=pRPDG;
1087 }
1088 else
1089 {
1090 if (cLPDG < pLPDG) nLPDG=cLPDG*1000+pLPDG*100-1;
1091 else if(cLPDG > pLPDG) nLPDG=pLPDG*1000+cLPDG*100-1;
1092 else nLPDG=pLPDG*1000+cLPDG*100-3;
1093 if ( cRPDG == pR1) nRPDG=-pR2;
1094 else nRPDG=-pR1; // cRPDG == pR2
1095 }
1096 }
1097 break;
1098 }
1099 break;
1100 case 3: // ------------------------> QDiQ/aQaDiQ + QaQ = QDiQ/aQaDiQ (single)
1101 switch (af)
1102 {
1103 case 1: // ....................... cLPDG > 7
1104 if(pLPDG < 0)
1105 {
1106 order= 1;
1107 if (pRPDG > cRPDG) nRPDG=pRPDG*1000+cRPDG*100+1;
1108 else if(pRPDG < cRPDG) nRPDG=cRPDG*1000+pRPDG*100+1;
1109 else nRPDG=cRPDG*1000+pRPDG*100+3;
1110 if (-pLPDG == L1) nLPDG=L2;
1111 else nLPDG=L1; // -pLPDG == L2
1112 }
1113 else // pRPDG < 0
1114 {
1115 order=-1;
1116 if (pLPDG > cRPDG) nLPDG=pLPDG*1000+cRPDG*100+1;
1117 else if(pLPDG < cRPDG) nLPDG=cRPDG*1000+pLPDG*100+1;
1118 else nLPDG=cRPDG*1000+pLPDG*100+3;
1119 if (-pRPDG == L1) nRPDG=L2;
1120 else nRPDG=L1; // -pRPDG == L2
1121 }
1122 break;
1123 case 2: // ....................... cRPDG > 7
1124 if(pLPDG < 0)
1125 {
1126 order=-1;
1127 if (pRPDG > cLPDG) nRPDG=pRPDG*1000+cLPDG*100+1;
1128 else if(pRPDG < cLPDG) nRPDG=cLPDG*1000+pRPDG*100+1;
1129 else nRPDG=cLPDG*1000+pRPDG*100+3;
1130 if (-pLPDG == R1) nLPDG=R2;
1131 else nLPDG=R1; // -pLPDG == R2
1132 }
1133 else // pRPDG < 0
1134 {
1135 order= 1;
1136 if (pLPDG > cLPDG) nLPDG=pLPDG*1000+cLPDG*100+1;
1137 else if(pLPDG < cLPDG) nLPDG=cLPDG*1000+pLPDG*100+1;
1138 else nLPDG=cLPDG*1000+pLPDG*100+3;
1139 if (-pRPDG == R1) nRPDG=R2;
1140 else nRPDG=R1; // -pRPDG == R2
1141 }
1142 break;
1143 case 3: // ....................... cLPDG <-7 (cRPDG <0)
1144 if(pLPDG > 0)
1145 {
1146 order= 1;
1147 if (pRPDG < cRPDG) nRPDG=pRPDG*1000+cRPDG*100-1;
1148 else if(pRPDG > cRPDG) nRPDG=cRPDG*1000+pRPDG*100-1;
1149 else nRPDG=cRPDG*1000+pRPDG*100-3;
1150 if ( pLPDG == L1) nLPDG=-L2;
1151 else nLPDG=-L1; // pLPDG == L2
1152 }
1153 else // pRPDG > 0
1154 {
1155 order=-1;
1156 if (pLPDG < cRPDG) nLPDG=pLPDG*1000+cRPDG*100-1;
1157 else if(pLPDG > cRPDG) nLPDG=cRPDG*1000+pLPDG*100-1;
1158 else nLPDG=cRPDG*1000+pLPDG*100-3;
1159 if ( pRPDG == L1) nRPDG=-L2;
1160 else nRPDG=-L1; // pRPDG == L2
1161 }
1162 break;
1163 case 4: // ....................... cRPDG <-7 (cLPDG <0)
1164 if(pLPDG > 0) // pRPDG & cLPDG are anti-quarks
1165 {
1166 order=-1;
1167 if (pRPDG < cLPDG) nRPDG=pRPDG*1000+cLPDG*100-1;
1168 else if(pRPDG > cLPDG) nRPDG=cLPDG*1000+pRPDG*100-1;
1169 else nRPDG=cLPDG*1000+pRPDG*100-3;
1170 if ( pLPDG == R1) nLPDG=-R2;
1171 else nLPDG=-R1; // pLPDG == R2
1172 }
1173 else // pRPDG > 0
1174 {
1175 order= 1;
1176 if (pLPDG < cLPDG) nLPDG=pLPDG*1000+cLPDG*100-1;
1177 else if(pLPDG > cLPDG) nLPDG=cLPDG*1000+pLPDG*100-1;
1178 else nLPDG=cLPDG*1000+pLPDG*100-3;
1179 if ( pRPDG == R1) nRPDG=-R2;
1180 else nRPDG=-R1; // pRPDG == R2
1181 }
1182 break;
1183 }
1184 break;
1185 case 4: // ------------------------------------> QaQ + aDiQDiQ = QaQ (double)
1186 switch (af)
1187 {
1188 case 1: // ....................... pLPDG > 7 && pRPDG <-7 === LL/RR
1189 order= 1;
1190 if(-cLPDG == pL1) nLPDG= pL2;
1191 else nLPDG= pL1;
1192 if( cRPDG == pR1) nRPDG=-pR2;
1193 else nRPDG=-pR1;
1194 break;
1195 case 2: // ...................... pLPDG > 7 && pRPDG <-7 === LR/RL
1196 order=-1;
1197 if(-cRPDG == pL1) nLPDG= pL2;
1198 else nLPDG= pL1;
1199 if( cLPDG == pR1) nRPDG=-pR2;
1200 else nRPDG=-pR1;
1201 break;
1202 case 3: // ...................... pRPDG > 7 && pLPDG <-7 === LL/RR
1203 order= 1;
1204 if( cLPDG == pL1) nLPDG=-pL2;
1205 else nLPDG=-pL1;
1206 if(-cRPDG == pR1) nRPDG= pR2;
1207 else nRPDG= pR1;
1208 break;
1209 case 4: // ...................... pRPDG > 7 && pLPDG <-7 === LR/RL
1210 order=-1;
1211 if( cRPDG == pL1) nLPDG=-pL2;
1212 else nLPDG=-pL1;
1213 if(-cLPDG == pR1) nRPDG= pR2;
1214 else nRPDG= pR1;
1215 break;
1216 }
1217 break;
1218 case 5: // ------------------------------------> aDiQDiQ + QaQ = QaQ (double)
1219 switch (af)
1220 {
1221 case 1: // ...................... cLPDG > 7 && cRPDG <-7 === LL/RR
1222 order= 1;
1223 if(-pLPDG == L1) nLPDG= L2;
1224 else nLPDG= L1;
1225 if( pRPDG == R1) nRPDG=-R2;
1226 else nRPDG=-R1;
1227 break;
1228 case 2: // ...................... cLPDG > 7 && cRPDG <-7 === LR/RL
1229 order=-1;
1230 if(-pRPDG == L1) nRPDG= L2;
1231 else nRPDG= L1;
1232 if( pLPDG == R1) nLPDG=-R2;
1233 else nLPDG=-R1;
1234 break;
1235 case 3: // ...................... cRPDG > 7 && cLPDG <-7 === LL/RR
1236 order= 1;
1237 if( pLPDG == L1) nLPDG=-L2;
1238 else nLPDG=-L1;
1239 if(-pRPDG == R1) nRPDG= R2;
1240 else nRPDG= R1;
1241 break;
1242 case 4: // ...................... cRPDG > 7 && cLPDG <-7 === LR/RL
1243 order=-1;
1244 if( pRPDG == L1) nRPDG=-L2;
1245 else nRPDG=-L1;
1246 if(-pLPDG == R1) nLPDG= R2;
1247 else nLPDG= R1;
1248 break;
1249 }
1250 break;
1251 case 6: // ------------------------------------> QDiQ + aQaDiQ = QaQ (double)
1252 switch (af)
1253 {
1254 case 1:
1255 order=-1;
1256 if(-cRPDG == pL1) nLPDG= pL2;
1257 else nLPDG= pL1;
1258 if( pRPDG == L1) nRPDG= -L2;
1259 else nRPDG= -L1;
1260 break;
1261 case 2:
1262 order= 1;
1263 if(-cLPDG == pL1) nLPDG= pL2;
1264 else nLPDG= pL1;
1265 if( pRPDG == R1) nRPDG= -R2;
1266 else nRPDG= -R1;
1267 break;
1268 case 3:
1269 order= 1;
1270 if(-cRPDG == pR1) nRPDG= pR2;
1271 else nRPDG= pR1;
1272 if( pLPDG == L1) nLPDG= -L2;
1273 else nLPDG= -L1;
1274 break;
1275 case 4:
1276 order=-1;
1277 if(-cLPDG == pR1) nRPDG= pR2;
1278 else nRPDG= pR1;
1279 if( pLPDG == R1) nLPDG= -R2;
1280 else nLPDG= -R1;
1281 break;
1282 case 5:
1283 order=-1;
1284 if(-pRPDG == L1) nRPDG= L2;
1285 else nRPDG= L1;
1286 if( cRPDG == pL1) nLPDG=-pL2;
1287 else nLPDG=-pL1;
1288 break;
1289 case 6:
1290 order= 1;
1291 if(-pLPDG == L1) nLPDG= L2;
1292 else nLPDG= L1;
1293 if( cRPDG == pR1) nRPDG=-pR2;
1294 else nRPDG=-pR1;
1295 break;
1296 case 7:
1297 order= 1;
1298 if(-pRPDG == R1) nRPDG= R2;
1299 else nRPDG= R1;
1300 if( cLPDG == pL1) nLPDG=-pL2;
1301 else nLPDG=-pL1;
1302 break;
1303 case 8:
1304 order=-1;
1305 if(-pLPDG == R1) nLPDG= R2;
1306 else nLPDG= R1;
1307 if( cLPDG == pR1) nRPDG=-pR2;
1308 else nRPDG=-pR1;
1309 break;
1310 }
1311 break;
1312 }
1313 if(!order) G4cerr<<"-Warning-G4QInel::Constr: t="<<tf<<", a="<<af<<", cL="<<cLPDG
1314 <<", cR="<<cRPDG<<", pL="<<pLPDG<<", pR="<<pRPDG<<G4endl;
1315 else
1316 {
1317 // With theNewHypotheticalPartons the min mass must be calculated & compared
1318 G4int LT=1;
1319 if(std::abs(nLPDG) > 7) ++LT;
1320 G4int RT=1;
1321 if(std::abs(nRPDG) > 7) ++RT;
1322 G4double minM=0.;
1323 G4bool sing=true;
1324 if(cLT==2 && cRT==2)
1325 {
1326 aLPDG=0;
1327 aRPDG=0;
1328 if(cLPDG>0)
1329 {
1330 aLPDG=nLPDG/100;
1331 aRPDG=(-nRPDG)/100;
1332 }
1333 else //cRPDG>0
1334 {
1335 aRPDG=nRPDG/100;
1336 aLPDG=(-nLPDG)/100;
1337 }
1338 G4int nL1=aLPDG/10;
1339 G4int nL2=aLPDG%10;
1340 G4int nR1=aRPDG/10;
1341 G4int nR2=aRPDG%10;
1342 if(nL1!=nR1 && nL1!=nR2 && nL2!=nR1 && nL2!=nR2) // Unreducable DiQ-aDiQ
1343 {
1344#ifdef debug
1345 G4cout<<"G4QIonIonCollis::Constr:aLPDG="<<aLPDG<<", aRPDG="<<aRPDG<<G4endl;
1346#endif
1347 sing=false;
1348 G4QPDGCode tmp;
1349 std::pair<G4int,G4int> pB=tmp.MakeTwoBaryons(nL1, nL2, nR1, nR2);
1350 minM=G4QPDGCode(pB.first).GetMass()+G4QPDGCode(pB.second).GetMass();
1351 }
1352 }
1353 if(sing)
1354 {
1355 std::pair<G4int,G4int> newPair = std::make_pair(nLPDG,nRPDG);
1356 G4QContent newStQC(newPair); // NewString QuarkContent
1357#ifdef debug
1358 G4cout<<"G4QIn::Con: LPDG="<<nLPDG<<",RPDG="<<nRPDG<<",QC="<<newStQC<<G4endl;
1359#endif
1360 G4int minPDG=newStQC.GetSPDGCode(); // PDG of the Lightest Hadron=String
1361 minM=G4QPDGCode(minPDG).GetMass() + mPi0; // Min SingleHadron=String Mass
1362 }
1363 // Compare this mass
1364 G4bool win=false;
1365 G4double pSM=0.;
1366 if(pSM2>0.) pSM=std::sqrt(pSM2);
1367 if(minC && pSM2 > maxiM2) // Up to now any positive mass is good
1368 {
1369 maxiM2=pSM2;
1370 win=true;
1371 }
1372 else if(!minC || pSM > minM)
1373 {
1374 pExcess=pSM-minM;
1375 if(minC || pExcess > excess)
1376 {
1377 minC=false;
1378 excess=pExcess;
1379 win=true;
1380 }
1381 }
1382 if(win)
1383 {
1384 sst=pst;
1385 sLPDG=nLPDG;
1386 sRPDG=nRPDG;
1387 sOrd=order;
1388 }
1389 } // End of IF(new partons are created)
1390 } // End of IF(compatible partons)
1391 //} // End of positive squared mass of the fused string
1392 } // End of the LOOP over the possible partners (with the exclusive if for itself)
1393 if(sOrd) // The best pStringCandidate was found
1394 {
1395 G4LorentzVector cL4M=cLeft->Get4Momentum();
1396 G4LorentzVector cR4M=cRight->Get4Momentum();
1397 G4QParton* pLeft=(*sst)->GetLeftParton();
1398 G4QParton* pRight=(*sst)->GetRightParton();
1399 G4LorentzVector pL4M=pLeft->Get4Momentum();
1400 G4LorentzVector pR4M=pRight->Get4Momentum();
1401#ifdef debug
1402 G4cout<<"G4QIonIonCollis::Constr:cS4M="<<cS4M<<" fused/w pS4M="<<pL4M+pR4M<<G4endl;
1403#endif
1404 if(sOrd>0)
1405 {
1406 pL4M+=cL4M;
1407 pR4M+=cR4M;
1408 }
1409 else
1410 {
1411 pL4M+=cR4M;
1412 pR4M+=cL4M;
1413 }
1414 pLeft->SetPDGCode(sLPDG);
1415 pLeft->Set4Momentum(pL4M);
1416 pRight->SetPDGCode(sRPDG);
1417 pRight->Set4Momentum(pR4M);
1418 delete (*ist);
1419 strings.erase(ist);
1420#ifdef debug
1421 G4LorentzVector ss4M=pL4M+pR4M;
1422 G4cout<<"G4QIonIonCollision::Constr:Created,4M="<<ss4M<<",m2="<<ss4M.m2()<<G4endl;
1423#endif
1424 if( ist != strings.begin() ) // To avoid going below begin() (for Windows)
1425 {
1426 ist--;
1427 con=false;
1428#ifdef debug
1429 G4cout<<"G4QIonIonCollision::Construct: *** IST Decremented ***"<<G4endl;
1430#endif
1431 }
1432 else
1433 {
1434 con=true;
1435#ifdef debug
1436 G4cout<<"G4QIonIonCollision::Construct: *** IST Begin ***"<<G4endl;
1437#endif
1438 break;
1439 }
1440 } // End of the IF(the best partnerString candidate was found)
1441 else
1442 {
1443#ifdef debug
1444 G4cout<<"-Warning-G4QInel::Const: S4M="<<cS4M<<",M2="<<cSM2<<" Leave ASIS"<<G4endl;
1445#endif
1446 ++problem;
1447 con=false;
1448 }
1449 }
1450 else con=false;
1451 } // End of loop over ist iterator
1452#ifdef debug
1453 G4cout<<"G4QIonIonCollision::Construct: *** IST While *** , con="<<con<<G4endl;
1454#endif
1455 } // End of "con" while
1456#ifdef edebug
1457 // This print has meaning only if something appear between it and the StringFragmLOOP
1458 G4LorentzVector t4M=theProjNucleus.Get4Momentum()+theTargNucleus.Get4Momentum();//NucInLS
1459 G4int rC=totChg-theProjNucleus.GetZ()-theTargNucleus.GetZ();
1460 G4int rB=totBaN-theProjNucleus.GetA()-theTargNucleus.GetA();
1461 G4int nStr=strings.size();
1462 G4cout<<"-EMCLS-G4QIn::Const: AfterSUPPRESION #ofS="<<nStr<<",tNuc4M(E=M)="<<sum<<G4endl;
1463 for(G4int i=0; i<nStr; i++)
1464 {
1465 G4LorentzVector strI4M=strings[i]->Get4Momentum();
1466 t4M+=strI4M;
1467 G4int sChg=strings[i]->GetCharge();
1468 rC-=sChg;
1469 G4int sBaN=strings[i]->GetBaryonNumber();
1470 rB-=sBaN;
1471 G4cout<<"-EMCLS-G4QIonIonCollision::Construct: St#"<<i<<", 4M="<<strI4M<<", M="
1472 <<strI4M.m()<<", C="<<sChg<<", B="<<sBaN<<G4endl;
1473 }
1474 G4cout<<"-EMCLS-G4QInel::Construct: r4M="<<t4M-totLS4M<<", rC="<<rC<<", rB="<<rB<<G4endl;
1475#endif
1476 //
1477 // --- If a problem is foreseen then the DiQaDiQ strings should be reduced if possible --
1478 //
1479#ifdef debug
1480 G4cout<<"G4QIonIonCollision::Construct: problem="<<problem<<G4endl;
1481#endif
1482 if(problem)
1483 {
1484 G4int nOfStr=strings.size();
1485#ifdef debug
1486 G4cout<<"G4QIonIonCollision::Constr:SecurityDiQaDiQReduction, #OfStr="<<nOfStr<<G4endl;
1487#endif
1488 for (G4int astring=0; astring < nOfStr; astring++)
1489 {
1490 G4QString* curString=strings[astring];
1491 G4QParton* cLeft=curString->GetLeftParton();
1492 G4QParton* cRight=curString->GetRightParton();
1493 G4int LT=cLeft->GetType();
1494 G4int RT=cRight->GetType();
1495 G4int sPDG=cLeft->GetPDGCode();
1496 G4int nPDG=cRight->GetPDGCode();
1497 if(LT==2 && RT==2)
1498 {
1499#ifdef debug
1500 G4cout<<"G4QIonIonCollis::Constr:TrySelfReduString,L="<<sPDG<<",R="<<nPDG<<G4endl;
1501#endif
1502 if( cLeft->ReduceDiQADiQ(cLeft, cRight) ) // DiQ-aDiQ pair was successfully reduced
1503 {
1504 sPDG=cLeft->GetPDGCode();
1505 nPDG=cRight->GetPDGCode();
1506#ifdef debug
1507 G4cout<<"+G4QInel::Const:#"<<astring<<" Reduced, L="<<sPDG<<", R="<<nPDG<<G4endl;
1508#endif
1509 }
1510#ifdef debug
1511 else G4cout<<"--*--G4QInel::Const: #"<<astring<<" DQ-aDQ reduction Failed"<<G4endl;
1512#endif
1513 } // End of the found DiQ/aDiQ pair
1514 else if(sPDG==3 && nPDG==-3)
1515 {
1516 sPDG= 1;
1517 nPDG=-1;
1518 cLeft->SetPDGCode(sPDG);
1519 cRight->SetPDGCode(nPDG);
1520 }
1521 else if(sPDG==-3 && nPDG==3)
1522 {
1523 sPDG=-1;
1524 nPDG= 1;
1525 cLeft->SetPDGCode(sPDG);
1526 cRight->SetPDGCode(nPDG);
1527 }
1528 }
1529 SwapPartons();
1530 } // End of IF(problem)
1531#ifdef edebug
1532 G4LorentzVector u4M=theProjNucleus.Get4Momentum()+theTargNucleus.Get4Momentum();//NucInLS
1533 G4int rCh=totChg-theProjNucleus.GetZ()-theTargNucleus.GetZ();
1534 G4int rBa=totBaN-theProjNucleus.GetA()-theTargNucleus.GetA();
1535 G4int nStri=strings.size();
1536 G4cout<<"-EMCLS-G4QIn::Const: FinalConstruct, #ofSt="<<nStri<<",tN4M(E=M)="<<t4M<<G4endl;
1537 for(G4int i=0; i<nStri; i++)
1538 {
1539 G4LorentzVector strI4M=strings[i]->Get4Momentum();
1540 u4M+=strI4M;
1541 G4int sChg=strings[i]->GetCharge();
1542 rCh-=sChg;
1543 G4int sBaN=strings[i]->GetBaryonNumber();
1544 rBa-=sBaN;
1545 G4cout<<"-EMCLS-G4QIonIonCollision::Construct: St#"<<i<<", 4M="<<strI4M<<", M="
1546 <<strI4M.m()<<", C="<<sChg<<", B="<<sBaN<<G4endl;
1547 }
1548 G4cout<<"-EMCLS-G4QInel::Construct: r4M="<<u4M-totLS4M<<",rC="<<rCh<<",rB="<<rBa<<G4endl;
1549#endif
1550}
LT
@ LT
Definition: Evaluator.cc:66
FatalException
@ FatalException
Definition: G4ExceptionSeverity.hh:60
G4QHadronVector
std::vector< G4QHadron * > G4QHadronVector
Definition: G4QHadronVector.hh:56
G4QInteractionVector
std::vector< G4QInteraction * > G4QInteractionVector
Definition: G4QInteractionVector.hh:54
G4QPartonPairVector
std::vector< G4QPartonPair * > G4QPartonPairVector
Definition: G4QPartonPairVector.hh:46
G4double
double G4double
Definition: G4Types.hh:64
G4int
int G4int
Definition: G4Types.hh:66
G4bool
bool G4bool
Definition: G4Types.hh:67
G4endl
#define G4endl
Definition: G4ios.hh:52
G4cerr
G4DLLIMPORT std::ostream G4cerr
G4cout
G4DLLIMPORT std::ostream G4cout
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:53
CLHEP::Hep3Vector::x
double x() const
CLHEP::Hep3Vector::y
double y() const
CLHEP::HepLorentzRotation::rotateY
HepLorentzRotation & rotateY(double delta)
Definition: LorentzRotation.cc:162
CLHEP::HepLorentzRotation::rotateZ
HepLorentzRotation & rotateZ(double delta)
Definition: LorentzRotation.cc:174
CLHEP::HepLorentzRotation::inverse
HepLorentzRotation inverse() const
CLHEP::HepLorentzVector::theta
double theta() const
CLHEP::HepLorentzVector::boost
HepLorentzVector & boost(double, double, double)
Definition: LorentzVector.cc:59
G4QCHIPSWorld::Get
static G4QCHIPSWorld * Get()
Definition: G4QCHIPSWorld.cc:79
G4QHadron::Get4Momentum
G4LorentzVector Get4Momentum() const
Definition: G4QHadron.hh:79
G4QHadron::GetMass
G4double GetMass() const
Definition: G4QHadron.hh:176
G4QHadron::GetBaryonNumber
G4int GetBaryonNumber() const
Definition: G4QHadron.hh:181
G4QHadron::IncrementCollisionCount
void IncrementCollisionCount(G4int aCount)
Definition: G4QHadron.hh:104
G4QHadron::GetNextParton
G4QParton * GetNextParton()
Definition: G4QHadron.cc:1614
G4QHadron::GetPDGCode
G4int GetPDGCode() const
Definition: G4QHadron.hh:170
G4QHadron::GetPosition
const G4ThreeVector & GetPosition() const
Definition: G4QHadron.hh:182
G4QHadron::GetColor
std::list< G4QParton * > GetColor()
Definition: G4QHadron.hh:93
G4QHadron::GetAntiColor
std::list< G4QParton * > GetAntiColor()
Definition: G4QHadron.hh:94
G4QHadron::GetNextAntiParton
G4QParton * GetNextAntiParton()
Definition: G4QHadron.cc:1622
G4QHadron::GetQC
G4QContent GetQC() const
Definition: G4QHadron.hh:173
G4QHadron::Set4Momentum
void Set4Momentum(const G4LorentzVector &aMom)
Definition: G4QHadron.hh:187
G4QInteraction::GetNumberOfSoftCollisions
G4int GetNumberOfSoftCollisions()
Definition: G4QInteraction.hh:59
G4QInteraction::SetNumberOfSoftCollisions
void SetNumberOfSoftCollisions(G4int nofSoft)
Definition: G4QInteraction.hh:65
G4QInteraction::GetProjectile
G4QHadron * GetProjectile() const
Definition: G4QInteraction.hh:57
G4QInteraction::GetTarget
G4QHadron * GetTarget() const
Definition: G4QInteraction.hh:58
G4QInteraction::SetTarget
void SetTarget(G4QHadron *aTarget)
Definition: G4QInteraction.hh:63
G4QInteraction::SetNumberOfDINRCollisions
void SetNumberOfDINRCollisions(G4int nofDINR)
Definition: G4QInteraction.hh:66
G4QInteraction::SetNumberOfDiffractiveCollisions
void SetNumberOfDiffractiveCollisions(G4int nofDiff)
Definition: G4QInteraction.hh:68
G4QIonIonCollision::ExciteSingDiffParticipants
G4bool ExciteSingDiffParticipants(G4QHadron *aPartner, G4QHadron *bPartner) const
Definition: G4QIonIonCollision.cc:3415
G4QIonIonCollision::ExciteDiffParticipants
G4bool ExciteDiffParticipants(G4QHadron *aPartner, G4QHadron *bPartner) const
Definition: G4QIonIonCollision.cc:3308
G4QNucleus::InitByPDG
void InitByPDG(G4int newPDG)
Definition: G4QNucleus.cc:371
G4QNucleus::GetN
G4int GetN() const
Definition: G4QNucleus.hh:71
G4QNucleus::GetZ
G4int GetZ() const
Definition: G4QNucleus.hh:70
G4QNucleus::ChooseImpactXandY
std::pair< G4double, G4double > ChooseImpactXandY(G4double maxImpact)
Definition: G4QNucleus.cc:3551
G4QNucleus::GetOuterRadius
G4double GetOuterRadius()
Definition: G4QNucleus.cc:3947
G4QNucleus::StartLoop
G4bool StartLoop()
Definition: G4QNucleus.cc:3982
G4QNucleus::SubtractNucleon
void SubtractNucleon(G4QHadron *pNucleon)
Definition: G4QNucleus.cc:612
G4QNucleus::DoLorentzBoost
void DoLorentzBoost(const G4LorentzVector &theBoost)
Definition: G4QNucleus.hh:155
G4QNucleus::GetA
G4int GetA() const
Definition: G4QNucleus.hh:73
G4QNucleus::DoLorentzRotation
void DoLorentzRotation(const G4LorentzRotation &theLoRot)
Definition: G4QNucleus.hh:160
G4QNucleus::DeleteNucleons
void DeleteNucleons()
Definition: G4QNucleus.cc:697
G4QNucleus::GetThickness
G4double GetThickness(G4double b)
Definition: G4QNucleus.cc:4048
G4QNucleus::GetNucleons4Momentum
G4LorentzVector GetNucleons4Momentum()
Definition: G4QNucleus.hh:107
G4QNucleus::Init3D
void Init3D()
Definition: G4QNucleus.cc:3910
G4QNucleus::GetNextNucleon
G4QHadron * GetNextNucleon()
Definition: G4QNucleus.hh:100
G4QNucleus::GetGSMass
G4double GetGSMass() const
Definition: G4QNucleus.hh:82
G4QPDGCode::GetMass
G4double GetMass()
Definition: G4QPDGCode.cc:693
G4QPDGCode::MakeTwoBaryons
std::pair< G4int, G4int > MakeTwoBaryons(G4int L1, G4int L2, G4int R1, G4int R2)
Definition: G4QPDGCode.cc:2413
G4QPartonPair::GetCollisionType
G4int GetCollisionType()
Definition: G4QPartonPair.hh:72
G4QParton::SetPDGCode
void SetPDGCode(G4int aPDG)
Definition: G4QParton.cc:130
G4QParton::GetQC
G4QContent GetQC() const
Definition: G4QParton.hh:82
G4QParton::GetType
const G4int & GetType() const
Definition: G4QParton.hh:88
G4QParton::GetPDGCode
G4int GetPDGCode() const
Definition: G4QParton.hh:81
G4QParton::Get4Momentum
const G4LorentzVector & Get4Momentum() const
Definition: G4QParton.hh:84
G4QParton::Set4Momentum
void Set4Momentum(const G4LorentzVector &aMomentum)
Definition: G4QParton.hh:75
G4QParton::ReduceDiQADiQ
G4bool ReduceDiQADiQ(G4QParton *d1, G4QParton *d2)
Definition: G4QParton.cc:246
G4QString::GetLeftParton
G4QParton * GetLeftParton() const
Definition: G4QString.hh:87
G4QString::GetCharge
G4int GetCharge() const
Definition: G4QString.hh:92
G4QString::Get4Momentum
G4LorentzVector Get4Momentum() const
Definition: G4QString.cc:124
G4QString::GetBaryonNumber
G4int GetBaryonNumber() const
Definition: G4QString.hh:93
G4QString::GetRightParton
G4QParton * GetRightParton() const
Definition: G4QString.hh:88
G4QString::Boost
void Boost(G4ThreeVector &Velocity)
Definition: G4QString.cc:152
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
sqr
T sqr(const T &x)
Definition: templates.hh:145
DBL_MAX
#define DBL_MAX
Definition: templates.hh:83

◆ ~G4QIonIonCollision()

G4QIonIonCollision::~G4QIonIonCollision ( )

Definition at line 1552 of file G4QIonIonCollision.cc.

1553{
1554 std::for_each(strings.begin(), strings.end(), DeleteQString() );
1555}

Member Function Documentation

◆ AnnihilationOrder()

G4int G4QIonIonCollision::AnnihilationOrder ( G4int  LS,
G4int  MS,
G4int  uP,
G4int  mP,
G4int  sP,
G4int  nP 
)
protected

Definition at line 3671 of file G4QIonIonCollision.cc.

3673{// ^L^^ Curent ^^R^
3674 G4int Ord=0;
3675 // Curent Partner
3676 if (LS==2 && MPS==2 ) // Fuse 2 Q-aQ strings to DiQ-aDiQ
3677 {
3678#ifdef debug
3679 G4cout<<"G4QIonIonCollision::AnnihilationOrder:QaQ/QaQ->DiQ-aDiQ, uPDG="<<uPDG<<G4endl;
3680#endif
3681 if ( (uPDG>0 && sPDG>0 && mPDG<0 && nPDG<0) ||
3682 (uPDG<0 && sPDG<0 && mPDG>0 && nPDG>0) ) Ord= 1; // LL/RR
3683 else if( (uPDG>0 && nPDG>0 && mPDG<0 && sPDG<0) ||
3684 (uPDG<0 && nPDG<0 && mPDG>0 && sPDG>0) ) Ord=-1; // LR/RL
3685 else G4cerr<<"-Warning-G4QIonIonCollision::AnnihilationOrder: Wrong 22 fusion, L="
3686 <<uPDG<<",R="<<mPDG<<",cL="<<sPDG<<",cR="<<nPDG<<G4endl;
3687 }
3688 else if ( LS==2 && MPS==3 )
3689 {
3690 if (uPDG > 7) // Fuse pLDiQ
3691 {
3692#ifdef debug
3693 G4cout<<"G4QIonIonCollision::AnnihOrder:pLDiQ, sPDG="<<sPDG<<", nPDG="<<nPDG<<G4endl;
3694#endif
3695 if ( sPDG<0 && (-sPDG==uPDG/1000 || -sPDG==(uPDG/100)%10) ) Ord= 1; // LL/RR
3696 else if( nPDG<0 && (-nPDG==uPDG/1000 || -nPDG==(uPDG/100)%10) ) Ord=-1; // LR/RL
3697 else G4cerr<<"-Warning-G4QIonIonCollision::AnnihilationOrder: Wrong pLDiQ, L="<<uPDG
3698 <<",R="<<mPDG<<",cL="<<sPDG<<",cR="<<nPDG<<G4endl;
3699 }
3700 else if (mPDG > 7) // Fuse pRDiQ
3701 {
3702#ifdef debug
3703 G4cout<<"G4QIonIonCollision::AnnihOrder:pRDiQ, sPDG="<<sPDG<<", nPDG="<<nPDG<<G4endl;
3704#endif
3705 if ( sPDG<0 && (-sPDG==mPDG/1000 || -sPDG==(mPDG/100)%10) ) Ord=-1; // LR/RL
3706 else if( nPDG<0 && (-nPDG==mPDG/1000 || -nPDG==(mPDG/100)%10) ) Ord= 1; // LL/RR
3707 else G4cerr<<"-Warning-G4QIonIonCollision::AnnihilationOrder: Wrong pRDiQ, L="<<uPDG
3708 <<",R="<<mPDG<<",cL="<<sPDG<<",cR="<<nPDG<<G4endl;
3709 }
3710 else if (uPDG <-7) // Fuse pLaDiQ
3711 {
3712#ifdef debug
3713 G4cout<<"G4QIonIonCollision::AnnihOrder:pLaDiQ, sPDG="<<sPDG<<",nPDG="<<nPDG<<G4endl;
3714#endif
3715 if ( sPDG>0 && (sPDG==(-uPDG)/1000 || sPDG==((-uPDG)/100)%10) ) Ord= 1; // LL/RR
3716 else if( nPDG>0 && (nPDG==(-uPDG)/1000 || nPDG==((-uPDG)/100)%10) ) Ord=-1; // LR/RL
3717 else G4cerr<<"-Warning-G4QIonIonCollision::AnnihilationOrder: Wrong pLaDiQ, L="<<uPDG
3718 <<", R="<<mPDG<<", cL="<<sPDG<<", cR="<<nPDG<<G4endl;
3719 }
3720 else if (mPDG <-7) // Fuse pRaDiQ
3721 {
3722#ifdef debug
3723 G4cout<<"G4QIonIonCollision::AnnihOrder:pRaDiQ, sPDG="<<sPDG<<",nPDG="<<nPDG<<G4endl;
3724#endif
3725 if ( sPDG>0 && (sPDG==(-mPDG)/1000 || sPDG==((-mPDG)/100)%10) ) Ord=-1; // LR/RL
3726 else if( nPDG>0 && (nPDG==(-mPDG)/1000 || nPDG==((-mPDG)/100)%10) ) Ord= 1; // LL/RR
3727 else G4cerr<<"-Warning-G4QIonIonCollision::AnnihilationOrder: Wrong pRaDiQ, L="<<uPDG
3728 <<", R="<<mPDG<<", cL="<<sPDG<<", cR="<<nPDG<<G4endl;
3729 }
3730 else if( (sPDG<0 && (-sPDG==mPDG || -sPDG==uPDG) ) ||
3731 (nPDG<0 && (-nPDG==mPDG || -nPDG==uPDG) ) ) Ord= 2; // @@ Annihilation fusion
3732#ifdef debug
3733 else G4cout<<"-Warning-G4QIonIonCollision::AnnihilatOrder:Wrong23StringFusion"<<G4endl;
3734 G4cout<<"G4QIonIonCollision::AnnihilationOrder: Ord="<<Ord<<",sPDG="<<sPDG<<",nPDG="
3735 <<nPDG<<", uPDG="<<uPDG<<",mPDG="<<mPDG<<G4endl;
3736#endif
3737 }
3738 else if ( LS==3 && MPS==2 )
3739 {
3740 if (sPDG > 7) // Fuse cLDiQ
3741 {
3742#ifdef debug
3743 G4cout<<"G4QIonIonCollision::AnnihOrder:cLDiQ, uPDG="<<uPDG<<", mPDG="<<mPDG<<G4endl;
3744#endif
3745 if ( uPDG<0 && (-uPDG==sPDG/1000 || -uPDG==(sPDG/100)%10) ) Ord= 1; // LL/RR
3746 else if( mPDG<0 && (-mPDG==sPDG/1000 || -mPDG==(sPDG/100)%10) ) Ord=-1; // LR/RL
3747 else G4cerr<<"-Warning-G4QIonIonCollision::AnnihilationOrder: Wrong cLDiQ, L="<<uPDG
3748 <<",R="<<mPDG<<",cL="<<sPDG<<",cR="<<nPDG<<G4endl;
3749 }
3750 else if (nPDG > 7) // Fuse cRDiQ
3751 {
3752#ifdef debug
3753 G4cout<<"G4QIonIonCollision::AnnihOrder:cRDiQ, uPDG="<<uPDG<<", mPDG="<<mPDG<<G4endl;
3754#endif
3755 if ( uPDG<0 && (-uPDG==nPDG/1000 || -uPDG==(nPDG/100)%10) ) Ord=-1; // LR/RL
3756 else if( mPDG<0 && (-mPDG==nPDG/1000 || -mPDG==(nPDG/100)%10) ) Ord= 1; // LL/RR
3757 else G4cerr<<"-Warning-G4QIonIonCollision::AnnihilationOrder: Wrong cRDiQ, L="<<uPDG
3758 <<",R="<<mPDG<<",cL="<<sPDG<<",cR="<<nPDG<<G4endl;
3759 }
3760 else if (sPDG <-7) // Fuse cLaDiQ
3761 {
3762#ifdef debug
3763 G4cout<<"G4QIonIonCollision::AnnihOrder:cLaDiQ, uPDG="<<uPDG<<",mPDG="<<mPDG<<G4endl;
3764#endif
3765 if ( uPDG>0 && (uPDG==(-sPDG)/1000 || uPDG==((-sPDG)/100)%10) ) Ord= 1; // LL/RR
3766 else if( mPDG>0 && (mPDG==(-sPDG)/1000 || mPDG==((-sPDG)/100)%10) ) Ord=-1; // LR/RL
3767 else G4cerr<<"-Warning-G4QIonIonCollision::AnnihilationOrder: Wrong cLaDiQ, L="<<uPDG
3768 <<", R="<<mPDG<<", cL="<<sPDG<<", cR="<<nPDG<<G4endl;
3769 }
3770 else if (nPDG <-7) // Fuse cRaDiQ
3771 {
3772#ifdef debug
3773 G4cout<<"G4QIonIonCollision::AnnihOrder:cRaDiQ, uPDG="<<uPDG<<",mPDG="<<mPDG<<G4endl;
3774#endif
3775 if ( uPDG>0 && (uPDG==(-nPDG)/1000 || uPDG==((-nPDG)/100)%10) ) Ord=-1; // LR/RL
3776 else if( mPDG>0 && (mPDG==(-nPDG)/1000 || mPDG==((-nPDG)/100)%10) ) Ord= 1; // LL/RR
3777 else G4cerr<<"-Warning-G4QIonIonCollision::AnnihilationOrder: Wrong cRaDiQ, L="<<uPDG
3778 <<", R="<<mPDG<<", cL="<<sPDG<<", cR="<<nPDG<<G4endl;
3779 }
3780 else if( (uPDG<0 && (-uPDG==sPDG || -uPDG==nPDG) ) ||
3781 (mPDG<0 && (-mPDG==sPDG || -mPDG==nPDG) ) ) Ord=2; // @@ Annihilation fusion
3782#ifdef debug
3783 else G4cout<<"-Warning-G4QIonIonCollision::AnnihilatOrder:Wrong32StringFusion"<<G4endl;
3784 G4cout<<"G4QIonIonCollision::AnnihilationOrder: Ord="<<Ord<<",sPDG="<<sPDG<<",nPDG="
3785 <<nPDG<<", uPDG="<<uPDG<<",mPDG="<<mPDG<<G4endl;
3786#endif
3787 }
3788 else if ( (LS==2 && MPS==4) || (LS==4 && MPS==2) )
3789 {
3790 if (uPDG > 7) // Annihilate pLDiQ
3791 {
3792#ifdef debug
3793 G4cout<<"G4QIonIonCollision::AnnihilOrder:pLDiQ,sPDG="<<sPDG<<",nPDG="<<nPDG<<G4endl;
3794#endif
3795 if ( sPDG<0 && (-sPDG==uPDG/1000 || -sPDG==(uPDG/100)%10) &&
3796 (nPDG==(-mPDG)/1000 || nPDG==((-mPDG)/100)%10) ) Ord= 1; // LL/RR
3797 else if( nPDG<0 && (-nPDG==uPDG/1000 || -nPDG==(uPDG/100)%10) &&
3798 (sPDG==(-mPDG)/1000 || sPDG==((-mPDG)/100)%10) ) Ord=-1; // LR/RL
3799 else G4cerr<<"-Warning-G4QIonIonCollision::AnnihilationOrder: Wrong 24 pLDiQ, L="
3800 <<uPDG<<",R="<<mPDG<<",cL="<<sPDG<<",cR="<<nPDG<<G4endl;
3801 }
3802 else if (mPDG >7) // Annihilate pRDiQ
3803 {
3804#ifdef debug
3805 G4cout<<"G4QIonIonCollision::AnnihilOrder:PRDiQ,sPDG="<<sPDG<<",nPDG="<<nPDG<<G4endl;
3806#endif
3807 if ( sPDG<0 && (-sPDG==mPDG/1000 || -sPDG==(mPDG/100)%10) &&
3808 (nPDG==(-uPDG)/1000 || nPDG==((-uPDG)/100)%10) ) Ord=-1; // LR/RL
3809 else if( nPDG<0 && (-nPDG==mPDG/1000 || -nPDG==(mPDG/100)%10) &&
3810 (sPDG==(-uPDG)/1000 || sPDG==((-uPDG)/100)%10) ) Ord= 1; // LL/RR
3811 else G4cerr<<"-Warning-G4QIonIonCollision::AnnihilationOrder: Wrong 24 pLDiQ, L="
3812 <<uPDG<<",R="<<mPDG<<",cL="<<sPDG<<",cR="<<nPDG<<G4endl;
3813 }
3814 else if (sPDG > 7) // Annihilate cLDiQ
3815 {
3816#ifdef debug
3817 G4cout<<"G4QIonIonCollision::AnnihilOrder:cLDiQ,uPDG="<<uPDG<<",mPDG="<<mPDG<<G4endl;
3818#endif
3819 if ( uPDG<0 && (-uPDG==sPDG/1000 || -uPDG==(sPDG/100)%10) &&
3820 (mPDG==(-nPDG)/1000 || mPDG==((-nPDG)/100)%10) ) Ord= 1; // LL/RR
3821 else if( mPDG<0 && (-mPDG==sPDG/1000 || -mPDG==(sPDG/100)%10) &&
3822 (uPDG==(-nPDG)/1000 || uPDG==((-nPDG)/100)%10) ) Ord=-1; // LR/RL
3823 else G4cerr<<"-Warning-G4QIonIonCollision::AnnihilationOrder: Wrong 24 cLDiQ, L="
3824 <<uPDG<<",R="<<mPDG<<",cL="<<sPDG<<",cR="<<nPDG<<G4endl;
3825 }
3826 else if (nPDG > 7) // Annihilate cRDiQ
3827 {
3828#ifdef debug
3829 G4cout<<"G4QIonIonCollision::AnnihilOrder:cRDiQ,uPDG="<<uPDG<<",mPDG="<<mPDG<<G4endl;
3830#endif
3831 if ( uPDG<0 && (-uPDG==nPDG/1000 || -uPDG==(nPDG/100)%10) &&
3832 (mPDG==(-sPDG)/1000 || mPDG==((-sPDG)/100)%10) ) Ord=-1; // LR/RL
3833 else if( mPDG<0 && (-mPDG==nPDG/1000 || -mPDG==(nPDG/100)%10) &&
3834 (uPDG==(-sPDG)/1000 || uPDG==((-sPDG)/100)%10) ) Ord= 1; // LL/RR
3835 else G4cerr<<"-Warning-G4QIonIonCollision::AnnihilationOrder: Wrong 24 cRDiQ, L="
3836 <<uPDG<<",R="<<mPDG<<",cL="<<sPDG<<",cR="<<nPDG<<G4endl;
3837 }
3838#ifdef debug
3839 else G4cout<<"-Warning-G4QIonIonCollision::AnnihilOrder: Wrong 24StringFusion"<<G4endl;
3840 G4cout<<"G4QIonIonCollision::AnnihilationOrder: Ord="<<Ord<<",sPDG="<<sPDG<<",nPDG="
3841 <<nPDG<<", uPDG="<<uPDG<<",mPDG="<<mPDG<<G4endl;
3842#endif
3843 }
3844 else if ( LS==3 && MPS==3 )
3845 {
3846 if (uPDG > 7) // Annihilate pLDiQ
3847 {
3848#ifdef debug
3849 G4cout<<"G4QIonIonCollision::AnnihilOrder:pLDiQ,sPDG="<<sPDG<<",nPDG="<<nPDG<<G4endl;
3850#endif
3851 if ( sPDG<-7 && (-nPDG==uPDG/1000 || -nPDG==(uPDG/100)%10) &&
3852 (mPDG==(-sPDG)/1000 || mPDG==((-sPDG)/100)%10) ) Ord=-1; // LR/RL
3853 else if( nPDG<-7 && (-sPDG==uPDG/1000 || -sPDG==(uPDG/100)%10) &&
3854 (mPDG==(-nPDG)/1000 || mPDG==((-nPDG)/100)%10) ) Ord= 1; // LL/RR
3855 else G4cerr<<"-Warning-G4QIonIonCollision::AnnihilationOrder: Wrong 33 pLDiQ, L="
3856 <<uPDG<<",R="<<mPDG<<",cL="<<sPDG<<",cR="<<nPDG<<G4endl;
3857 }
3858 else if(mPDG > 7) // Annihilate pRDiQ
3859 {
3860#ifdef debug
3861 G4cout<<"G4QIonIonCollision::AnnihilOrder:pRDiQ,sPDG="<<sPDG<<",nPDG="<<nPDG<<G4endl;
3862#endif
3863 if ( sPDG<-7 && (-nPDG==mPDG/1000 || -nPDG==(mPDG/100)%10) &&
3864 (uPDG==(-sPDG)/1000 || uPDG==((-sPDG)/100)%10) ) Ord= 1; // LL/RR
3865 else if( nPDG<-7 && (-sPDG==mPDG/1000 || -sPDG==(mPDG/100)%10) &&
3866 (uPDG==(-nPDG)/1000 || uPDG==((-nPDG)/100)%10) ) Ord=-1; // LR/RL
3867 else G4cerr<<"-Warning-G4QIonIonCollision::AnnihilationOrder: Wrong33pRDiQ, L="<<uPDG
3868 <<",R="<<mPDG<<",cL="<<sPDG<<",cR="<<nPDG<<G4endl;
3869 }
3870 else if(sPDG > 7) // Annihilate cLDiQ
3871 {
3872#ifdef debug
3873 G4cout<<"G4QIonIonCollision::AnnihilOrder:cLDiQ,uPDG="<<uPDG<<",mPDG="<<mPDG<<G4endl;
3874#endif
3875 if ( uPDG<-7 && (-mPDG==sPDG/1000 || -mPDG==(sPDG/100)%10) &&
3876 (nPDG==(-uPDG)/1000 || nPDG==((-uPDG)/100)%10) ) Ord=-1; // LR/RL
3877 else if( mPDG<-7 && (-uPDG==sPDG/1000 || -uPDG==(sPDG/100)%10) &&
3878 (nPDG==(-mPDG)/1000 || nPDG==((-mPDG)/100)%10) ) Ord= 1; // LL/RR
3879 else G4cerr<<"-Warning-G4QIonIonCollision::AnnihilationOrder: Wrong 33 cLDiQ, L="
3880 <<uPDG<<",R="<<mPDG<<",cL="<<sPDG<<",cR="<<nPDG<<G4endl;
3881 }
3882 else if(nPDG > 7) // Annihilate cRDiQ
3883 {
3884#ifdef debug
3885 G4cout<<"G4QIonIonCollision::AnnihilOrder:cRDiQ,uPDG="<<uPDG<<",mPDG="<<mPDG<<G4endl;
3886#endif
3887 if ( uPDG<-7 && (-mPDG==nPDG/1000 || -mPDG==(nPDG/100)%10) &&
3888 (nPDG==(-uPDG)/1000 || nPDG==((-uPDG)/100)%10) ) Ord= 1; // LL/RR
3889 else if( mPDG<-7 && (-uPDG==nPDG/1000 || -sPDG==(nPDG/100)%10) &&
3890 (sPDG==(-mPDG)/1000 || sPDG==((-mPDG)/100)%10) ) Ord=-1; // LR/RL
3891 else G4cerr<<"-Warning-G4QIonIonCollision::AnnihilationOrder: Wrong 33 cRDiQ, L="
3892 <<uPDG<<",R="<<mPDG<<",cL="<<sPDG<<",cR="<<nPDG<<G4endl;
3893 }
3894#ifdef debug
3895 else G4cout<<"-Warning-G4QIonIonCollision::AnnihilOrder: Wrong 33StringFusion"<<G4endl;
3896 G4cout<<"G4QIonIonCollision::AnnihilationOrder: Ord="<<Ord<<",sPDG="<<sPDG<<",nPDG="
3897 <<nPDG<<", uPDG="<<uPDG<<",mPDG="<<mPDG<<G4endl;
3898#endif
3899 }
3900 return Ord;
3901}

Referenced by Breeder().

◆ Breeder()

void G4QIonIonCollision::Breeder ( )
protected

! Try to fragment the new String !!

Try to fragment the new String !

Definition at line 1791 of file G4QIonIonCollision.cc.

1792{ // This is the member function, which returns the resulting vector of Hadrons & Quasmons
1793 static const G4double eps = 0.001; // Tolerance in MeV
1794 //
1795 // ------------ At this point the strings are fragmenting to hadrons in LS -------------
1796 //
1797#ifdef edebug
1798 G4LorentzVector totLS4M=theProjNucleus.Get4Momentum()+theTargNucleus.Get4Momentum(); //LS
1799 G4int totChg=theProjNucleus.GetZ()+theTargNucleus.GetZ();
1800 G4int totBaN=theProjNucleus.GetA()+theTargNucleus.GetA();
1801 G4int nStri=strings.size();
1802 G4cout<<"-EMCLS-G4QIn::Breed: CHECKRecovery #ofS="<<nStri<<",N4M(E=M)="<<totLS4M<<G4endl;
1803 for(G4int i=0; i<nStri; i++)
1804 {
1805 G4LorentzVector strI4M=strings[i]->Get4Momentum();
1806 totLS4M+=strI4M;
1807 G4int sChg=strings[i]->GetCharge();
1808 totChg+=sChg;
1809 G4int sBaN=strings[i]->GetBaryonNumber();
1810 totBaN+=sBaN;
1811 G4cout<<"-EMCLS-G4QIonIonCollision::Breeder: St#"<<i<<", 4M="<<strI4M<<", M="
1812 <<strI4M.m()<<", C="<<sChg<<", B="<<sBaN<<G4endl;
1813 }
1814#endif
1815 G4int nOfStr=strings.size();
1816#ifdef debug
1817 G4cout<<"G4QIonIonCollision::Breeder: BeforeFragmentation, #OfStr="<<nOfStr<<G4endl;
1818#endif
1819 G4LorentzVector ft4M(0.,0.,0.,0.);
1820 G4QContent ftQC(0,0,0,0,0,0);
1821 G4bool ftBad=false;
1822 for(G4int i=0; i < nOfStr; ++i)
1823 {
1824 G4LorentzVector pS4M=strings[i]->Get4Momentum(); // String 4-momentum
1825 ft4M+=pS4M;
1826 G4QContent pSQC=strings[i]->GetQC(); // String Quark Content
1827 ftQC+=pSQC;
1828 if(pS4M.m2() < 0.) ftBad=true;
1829#ifdef debug
1830 G4cout<<">G4QIonIonCollision::Breed:1stTest,S#"<<i<<",4M="<<pS4M<<",QC="<<pSQC<<G4endl;
1831#endif
1832 }
1833 if(ftBad)
1834 {
1835 G4Quasmon* stringQuasmon = new G4Quasmon(ftQC, ft4M);
1836#ifdef debug
1837 G4cout<<"->G4QIonIonCollision::Breed:*TotQ*,QC="<<ftQC<<",4M="<<ft4M<<ft4M.m()<<G4endl;
1838#endif
1839 theQuasmons.push_back(stringQuasmon);
1840 G4LorentzVector rp4M=theProjNucleus.Get4Momentum(); // Nucleus 4-momentum in LS
1841 G4int rpPDG=theProjNucleus.GetPDG();
1842 G4QHadron* resPNuc = new G4QHadron(rpPDG,rp4M);
1843 theResult->push_back(resPNuc); // Fill the residual projectile nucleus
1844 G4LorentzVector rt4M=theTargNucleus.Get4Momentum(); // Nucleus 4-momentum in LS
1845 G4int rtPDG=theTargNucleus.GetPDG();
1846 G4QHadron* resTNuc = new G4QHadron(rtPDG,rt4M);
1847 theResult->push_back(resTNuc); // Fill the residual target nucleus
1848 return;
1849 }
1850 for (G4int astring=0; astring < nOfStr; astring++)
1851 {
1852#ifdef edebug
1853 G4LorentzVector sum=theProjNucleus.Get4Momentum()+theTargNucleus.Get4Momentum(); //inLS
1854 G4int rChg=totChg-theProjNucleus.GetZ()-theTargNucleus.GetZ();
1855 G4int rBaN=totBaN-theProjNucleus.GetA()-theTargNucleus.GetA();
1856 G4int nOfHadr=theResult->size();
1857 G4cout<<"-EMCLS-G4QIonIonCollision::Breed:#ofSt="<<nOfStr<<",#ofHad="<<nOfHadr<<G4endl;
1858 for(G4int i=astring; i<nOfStr; i++)
1859 {
1860 G4LorentzVector strI4M=strings[i]->Get4Momentum();
1861 sum+=strI4M;
1862 G4int sChg=strings[i]->GetCharge();
1863 rChg-=sChg;
1864 G4int sBaN=strings[i]->GetBaryonNumber();
1865 rBaN-=sBaN;
1866 G4cout<<"-EMCLS-G4QI::Breed:S#"<<i<<",4M="<<strI4M<<",C="<<sChg<<",B="<<sBaN<<G4endl;
1867 }
1868 for(G4int i=0; i<nOfHadr; i++)
1869 {
1870 G4LorentzVector hI4M=(*theResult)[i]->Get4Momentum();
1871 sum+=hI4M;
1872 G4int hChg=(*theResult)[i]->GetCharge();
1873 rChg-=hChg;
1874 G4int hBaN=(*theResult)[i]->GetBaryonNumber();
1875 rBaN-=hBaN;
1876 G4cout<<"-EMCLS-G4QIn::Breed: H#"<<i<<",4M="<<hI4M<<",C="<<hChg<<",B="<<hBaN<<G4endl;
1877 }
1878 G4cout<<"....-EMCLS-G4QInel::Br:r4M="<<sum-totLS4M<<",rC="<<rChg<<",rB="<<rBaN<<G4endl;
1879#endif
1880 G4QString* curString=strings[astring];
1881 if(!curString->GetDirection()) continue; // Historic for the dead strings: DoesNotWork
1882#ifdef edebug
1883 G4int curStrChg = curString->GetCharge();
1884 G4int curStrBaN = curString->GetBaryonNumber();
1885#endif
1886 G4LorentzVector curString4M = curString->Get4Momentum();
1887#ifdef debug
1888 G4cout<<"=--=>G4QIonIonCollision::Breeder: String#"<<astring<<",s4M/m="<<curString4M
1889 <<curString4M.m()<<", LPart="<<curString->GetLeftParton()->GetPDGCode()
1890 <<", RPart="<<curString->GetRightParton()->GetPDGCode()<<G4endl;
1891#endif
1892 G4QHadronVector* theHadrons = 0; // Prototype of theStringFragmentationOUTPUT
1893 theHadrons=curString->FragmentString(true);// !! Fragmenting the String !!
1894 if (!theHadrons) // The string can not be fragmented
1895 {
1896 // First try to correct the diQ-antiDiQ strings, converting them to Q-antiQ
1897 G4QParton* cLeft=curString->GetLeftParton();
1898 G4QParton* cRight=curString->GetRightParton();
1899 G4int sPDG=cLeft->GetPDGCode();
1900 G4int nPDG=cRight->GetPDGCode();
1901 G4int LT=cLeft->GetType();
1902 G4int RT=cRight->GetType();
1903 G4int LS=LT+RT;
1904 if(LT==2 && RT==2)
1905 {
1906#ifdef debug
1907 G4cout<<"G4QIonIonCollision::Breed:TryReduceString, L="<<sPDG<<",R="<<nPDG<<G4endl;
1908#endif
1909 if( cLeft->ReduceDiQADiQ(cLeft, cRight) ) // DiQ-aDiQ pair was successfully reduced
1910 {
1911 LT=1;
1912 RT=1;
1913 LS=2;
1914 sPDG=cLeft->GetPDGCode();
1915 nPDG=cRight->GetPDGCode();
1916#ifdef debug
1917 G4cout<<"G4QIonIonCollision::Breed:AfterReduction,L="<<sPDG<<",R="<<nPDG<<G4endl;
1918#endif
1919 theHadrons=curString->FragmentString(true);//!! Try to fragment the new String !!
1920 cLeft=curString->GetLeftParton();
1921 cRight=curString->GetRightParton();
1922#ifdef debug
1923 G4cout<<"G4QInel::Breed:L="<<cLeft->Get4Momentum()<<",R="<<cRight->Get4Momentum()
1924 <<G4endl;
1925#endif
1926 }
1927#ifdef debug
1928 else G4cout<<"^G4QIonIonCollision::Breed: DQ-aDQ reduction to Q-aQ Failed"<<G4endl;
1929#endif
1930 } // End of the SelfReduction
1931#ifdef debug
1932 G4cout<<"G4QIonIonCollision::Breed:AfterRedAttempt, theH="<<theHadrons<<", L4M="
1933 <<cLeft->Get4Momentum()<<", R4M="<<cRight->Get4Momentum()<<G4endl;
1934#endif
1935 unsigned next=astring+1; // The next string position
1936 if (!theHadrons) // The string can not be fragmented
1937 {
1938 G4int fusionDONE=0; // StringFusion didn't happen (1=Fuse L+L/R+R, -1=Fuse L+R/R+L)
1939 if(next < strings.size()) // TheString isn't theLastString can fuse
1940 {
1941 G4int fustr=0; // The found partner index (never can be 0)
1942 G4int swap=0; // working interger for swapping parton PDG
1943 G4double Vmin=DBL_MAX; // Prototype of the found Velocity Distance
1944 G4int dPDG=nPDG;
1945 G4int qPDG=sPDG;
1946 if(dPDG<-99 || (dPDG>0&&dPDG<7) || qPDG>99 || (qPDG<0 && qPDG>-7))
1947 {
1948 swap=qPDG;
1949 qPDG=dPDG;
1950 dPDG=swap;
1951 }
1952 if(dPDG>99) dPDG/=100;
1953 if(qPDG<-99) qPDG=-(-qPDG)/100;
1954#ifdef debug
1955 G4cout<<"G4QIonIonCollision::Breed:TryFuseStringS, q="<<qPDG<<", a="<<dPDG
1956 <<", n="<<next<<G4endl;
1957#endif
1958 G4ThreeVector curV=curString4M.vect()/curString4M.e();
1959 G4int reduce=0; // a#of reduced Q-aQ pairs
1960 G4int restr=0; // To use beyon the LOOP for printing
1961 G4int MPS=0; // PLS for the selected string
1962 for (restr=next; restr < nOfStr; restr++)
1963 {
1964 reduce=0;
1965 G4QString* reString=strings[restr];
1966 G4QParton* Left=reString->GetLeftParton();
1967 G4QParton* Right=reString->GetRightParton();
1968 G4int uPDG=Left->GetPDGCode();
1969 G4int mPDG=Right->GetPDGCode();
1970 G4int PLT =Left->GetType();
1971 G4int PRT =Right->GetType();
1972 G4int aPDG=mPDG;
1973 G4int rPDG=uPDG;
1974 if(aPDG<-99 || (aPDG>0 && aPDG<7) || rPDG>99 || (rPDG<0 && rPDG>-7))
1975 {
1976 swap=rPDG;
1977 rPDG=aPDG;
1978 aPDG=swap;
1979 }
1980 if(aPDG > 99) aPDG/=100;
1981 if(rPDG <-99) rPDG=-(-rPDG)/100;
1982 // Try to reduce two DQ-aDQ strings
1983#ifdef debug
1984 G4cout<<"G4Qnel::Breed: TryReduce #"<<restr<<", q="<<rPDG<<",a="<<aPDG<<G4endl;
1985#endif
1986 if(LT==2 && RT==2 && PLT==2 && PRT==2) // Have a chance for the reduction
1987 {
1988 G4int cQ1=(-qPDG)/10;
1989 G4int cQ2=(-qPDG)%10;
1990 G4int cA1=dPDG/10;
1991 G4int cA2=dPDG%10;
1992 G4int pQ1=(-rPDG)/10;
1993 G4int pQ2=(-rPDG)%10;
1994 G4int pA1=aPDG/10;
1995 G4int pA2=aPDG%10;
1996#ifdef debug
1997 G4cout<<"G4QIonIonCollision::Breeder: cQ="<<cQ1<<","<<cQ2<<", cA="<<cA1<<","
1998 <<cA2<<", pQ="<<pQ1<<","<<pQ2<<", pA="<<pA1<<","<<pA2<<G4endl;
1999#endif
2000 G4bool iQA = (cA1==pQ1 || cA1==pQ2 || cA2==pQ1 || cA2==pQ2);
2001 G4bool iAQ = (cQ1==pA1 || cQ1==pA2 || cQ2==pA1 || cQ2==pA2);
2002 if(iQA) reduce++;
2003 if(iAQ) reduce++;
2004 if (reduce==2) // Two quark pairs can be reduced
2005 {
2006 if(sPDG>0 && uPDG<0) // LL/RR Reduction
2007 {
2008 std::pair<G4int,G4int> resLL=ReducePair(sPDG/100, (-uPDG)/100);
2009 G4int newCL=resLL.first;
2010 G4int newPL=resLL.second;
2011 if(!newCL || !newPL)
2012 {
2013 G4cerr<<"*G4QIonIonCollision::Breed:CL="<<newCL<<",PL="<<newPL<<G4endl;
2014 G4Exception("G4QIonIonCollision::Breeder:","72",FatalException,"2LL-");
2015 }
2016 std::pair<G4int,G4int> resRR=ReducePair((-nPDG)/100, mPDG/100);
2017 G4int newCR=resRR.first;
2018 G4int newPR=resRR.second;
2019 if(!newCR || !newPR)
2020 {
2021 G4cerr<<"*G4QIonIonCollision::Breed:CR="<<newCR<<",PR="<<newPR<<G4endl;
2022 G4Exception("G4QIonIonCollision::Breeder:","72",FatalException,"2RR-");
2023 }
2024 cLeft->SetPDGCode(newCL); // Reset the left quark of curString
2025 cRight->SetPDGCode(-newCR); // Reset the right quark of curString
2026 Left->SetPDGCode(-newPL); // Reset the left quark of reString
2027 Right->SetPDGCode(newPR); // Reset the right quark of reString
2028 break; // Break out of the reString internal LOOP
2029 }
2030 else if(sPDG<0 && uPDG>0) // LL/RR Reduction
2031 {
2032 std::pair<G4int,G4int> resLL=ReducePair((-sPDG)/100, uPDG/100);
2033 G4int newCL=resLL.first;
2034 G4int newPL=resLL.second;
2035 if(!newCL || !newPL)
2036 {
2037 G4cerr<<"*G4QIonIonCollision::Breed:CL="<<newCL<<",PL="<<newPL<<G4endl;
2038 G4Exception("G4QIonIonCollision::Breeder:","72",FatalException,"2LL+");
2039 }
2040 std::pair<G4int,G4int> resRR=ReducePair(nPDG/100, (-mPDG)/100);
2041 G4int newCR=resRR.first;
2042 G4int newPR=resRR.second;
2043 if(!newCR || !newPR)
2044 {
2045 G4cerr<<"*G4QIonIonCollision::Breed:CR="<<newCR<<",PR="<<newPR<<G4endl;
2046 G4Exception("G4QIonIonCollision::Breeder:","72",FatalException,"2RR+");
2047 }
2048 cLeft->SetPDGCode(-newCL); // Reset the left quark of curString
2049 cRight->SetPDGCode(newCR); // Reset the right quark of curString
2050 Left->SetPDGCode(newPL); // Reset the left quark of reString
2051 Right->SetPDGCode(-newPR); // Reset the right quark of reString
2052 break; // Break out of the reString internal LOOP
2053 }
2054 else if(sPDG>0 && mPDG<0) // LR Reduction
2055 {
2056 std::pair<G4int,G4int> resLL=ReducePair(sPDG/100, (-mPDG)/100);
2057 G4int newCL=resLL.first;
2058 G4int newPR=resLL.second;
2059 if(!newCL || !newPR)
2060 {
2061 G4cerr<<"*G4QIonIonCollision::Breed:CL="<<newCL<<",PR="<<newPR<<G4endl;
2062 G4Exception("G4QIonIonCollision::Breeder:","72",FatalException,"2-LR");
2063 }
2064 std::pair<G4int,G4int> resRR=ReducePair((-nPDG)/100, uPDG/100);
2065 G4int newCR=resRR.first;
2066 G4int newPL=resRR.second;
2067 if(!newCR || !newPR)
2068 {
2069 G4cerr<<"*G4QIonIonCollision::Breed:CR="<<newCR<<",PL="<<newPL<<G4endl;
2070 G4Exception("G4QIonIonCollision::Breeder:","72",FatalException,"2-LR");
2071 }
2072 cLeft->SetPDGCode(newCL); // Reset the left quark of curString
2073 cRight->SetPDGCode(-newCR); // Reset the right quark of curString
2074 Left->SetPDGCode(newPL); // Reset the left quark of reString
2075 Right->SetPDGCode(-newPR); // Reset the right quark of reString
2076 break; // Break out of the reString internal LOOP
2077 }
2078 else // RL Reduction
2079 {
2080 std::pair<G4int,G4int> resLL=ReducePair((-sPDG)/100, mPDG/100);
2081 G4int newCL=resLL.first;
2082 G4int newPR=resLL.second;
2083 if(!newCL || !newPR)
2084 {
2085 G4cerr<<"*G4QIonIonCollision::Breed:CL="<<newCL<<",PR="<<newPR<<G4endl;
2086 G4Exception("G4QIonIonCollision::Breeder:","72",FatalException,"2-RL");
2087 }
2088 std::pair<G4int,G4int> resRR=ReducePair(nPDG/100, (-uPDG)/100);
2089 G4int newCR=resRR.first;
2090 G4int newPL=resRR.second;
2091 if(!newCR || !newPR)
2092 {
2093 G4cerr<<"*G4QIonIonCollision::Breed:CR="<<newCR<<",PL="<<newPL<<G4endl;
2094 G4Exception("G4QIonIonCollision::Breeder:","72",FatalException,"2-RL");
2095 }
2096 cLeft->SetPDGCode(-newCL); // Reset the left quark of curString
2097 cRight->SetPDGCode(newCR); // Reset the right quark of curString
2098 Left->SetPDGCode(-newPL); // Reset the left quark of reString
2099 Right->SetPDGCode(newPR); // Reset the right quark of reString
2100 break; // Break out of the reString internal LOOP
2101 }
2102 } // End of IF(possible reduction)
2103 }
2104 // Check StringsCanBeFused: all QaQ+QaQ(22), single QaQ+QDiQ/aQaDtQ(23/32),
2105 // double QaQ+DiQaDiQ(24/42), QDiQ+aDiQaQ(34/43)
2106#ifdef debug
2107 G4cout<<"G4QInel::Breed: TryFuse/w #"<<restr<<",q="<<rPDG<<",a="<<aPDG<<G4endl;
2108#endif
2109 G4int PLS=PLT+PRT;
2110 if( (LS==2 && PLS==2) || // QaQ+QaQ always to DiQaDiQ
2111 ( ( (LS==2 && PLS==3) || (LS==3 && PLS==2) ) &&// QaQ w QDiQ/aQaDiQ(single)
2112 ( (aPDG> 7 && (-dPDG==aPDG/10 || -dPDG==aPDG%10) ) || // cAQ w DQ
2113 (dPDG> 7 && (-aPDG==dPDG/10 || -aPDG==dPDG%10) ) || // AQ w cDQ
2114 (rPDG<-7 && (qPDG==(-rPDG)/10 || qPDG==(-rPDG)%10) ) || // cQ w ADQ
2115 (qPDG<-7 && (rPDG==(-qPDG)/10 || rPDG==(-qPDG)%10) ) // Q w cADQ
2116 //|| (aPDG< 0 && -aPDG==qPDG) || (dPDG< 0 && -dPDG==rPDG) // aQ w Q
2117 )
2118 ) || // -----------------------
2119 ( ( LS==2 && PLS==4 &&// QaQ w DiQaDiQ (double)
2120 (aPDG> 7 && (-dPDG == aPDG/10 || -dPDG == aPDG%10) ) &&
2121 (rPDG<-7 && (qPDG==(-rPDG)/10 || qPDG==(-rPDG)%10) )
2122 ) ||
2123 ( LS==4 && PLS==2 &&// DiQaDiQ w QaQ (double)
2124 (dPDG> 7 && (-aPDG == dPDG/10 || -aPDG == dPDG%10) ) &&
2125 (qPDG<-7 && (rPDG==(-qPDG)/10 || rPDG==(-qPDG)%10) )
2126 )
2127 ) || // -----------------------
2128 ( LS==3 && PLS==3 &&// QDiQ w aDiQaQ (double)
2129 ( (aPDG> 7 && (-dPDG == aPDG/10 || -dPDG == aPDG%10) &&
2130 qPDG<-7 && (rPDG==(-qPDG)/10 || rPDG==(-qPDG)%10)
2131 ) ||
2132 (dPDG> 7 && (-aPDG == dPDG/10 || -aPDG == dPDG%10) &&
2133 rPDG<-7 && (dPDG==(-rPDG)/10 || dPDG==(-rPDG)%10)
2134 )
2135 )
2136 )
2137 )
2138 {
2139 G4LorentzVector reString4M = reString->Get4Momentum();
2140 G4ThreeVector reV = reString4M.vect()/reString4M.e();
2141 G4double dV=(curV-reV).mag2(); // Squared difference between relVelocities
2142#ifdef debug
2143 G4cout<<"G4QIonIonCollision::Breeder: StringCand#"<<restr<<", q="<<rPDG
2144 <<", a="<<aPDG<<", L="<<uPDG<<", R="<<mPDG<<",dV="<<dV<<G4endl;
2145#endif
2146 if(dV < Vmin)
2147 {
2148 Vmin=dV;
2149 fustr=restr;
2150 MPS=PLS;
2151 }
2152 }
2153 }
2154 if(reduce==2) // String mutual reduction happened
2155 {
2156#ifdef debug
2157 G4cout<<"-G4QIonIonCollision::Breed:Reduced #"<<astring<<" & #"<<restr<<G4endl;
2158#endif
2159 astring--; // String was QCreduced using another String
2160 continue; // Repeat fragmentation of the same string
2161 }
2162 if(fustr) // The partner was found -> fuse strings
2163 {
2164#ifdef debug
2165 G4cout<<"G4QInel::Breeder: StPartner#"<<fustr<<", LT="<<LT<<",RT="<<RT<<G4endl;
2166#endif
2167 G4QString* fuString=strings[fustr];
2168 G4QParton* Left=fuString->GetLeftParton();
2169 G4QParton* Right=fuString->GetRightParton();
2170 G4int uPDG=Left->GetPDGCode();
2171 G4int mPDG=Right->GetPDGCode();
2172 G4int rPDG=uPDG;
2173 G4int aPDG=mPDG;
2174 if(aPDG<-99 || (aPDG>0 && aPDG<7) || rPDG>99 || (rPDG<0 && rPDG>-7))
2175 {
2176 swap=rPDG;
2177 rPDG=aPDG;
2178 aPDG=swap;
2179 }
2180 if(aPDG > 99) aPDG/=100;
2181 if(rPDG <-99) rPDG=-(-rPDG)/100;
2182 // Check that the strings can fuse (no anti-diquarks assumed)
2183 G4LorentzVector resL4M(0.,0.,0.,0.);
2184 G4LorentzVector resR4M(0.,0.,0.,0.);
2185 G4LorentzVector L4M=cLeft->Get4Momentum();
2186 G4LorentzVector R4M=cRight->Get4Momentum();
2187#ifdef debug
2188 G4cout<<"G4QIonIonCollision::Breeder:BeforeFuDir,sL="<<sPDG<<",nR="<<nPDG
2189 <<",uL="<<uPDG<<",mR="<<mPDG<<",L4M="<<L4M<<",R4M="<<R4M<<G4endl;
2190#endif
2191 G4LorentzVector PL4M=Left->Get4Momentum();
2192 G4LorentzVector PR4M=Right->Get4Momentum();
2193 fusionDONE=AnnihilationOrder(LS, MPS, uPDG, mPDG, sPDG, nPDG);
2194 if (fusionDONE>0)
2195 {
2196 if(fusionDONE>1) // Anihilation algorithm
2197 {
2198 if ( (uPDG<0 || nPDG<0) && -uPDG==nPDG ) Left->SetPDGCode(sPDG);
2199 else if( (mPDG<0 || sPDG<0) && -mPDG==sPDG ) Right->SetPDGCode(nPDG);
2200 else if( (uPDG<0 || sPDG<0) && -uPDG==sPDG ) Left->SetPDGCode(nPDG);
2201 else if( (mPDG<0 || nPDG<0) && -mPDG==nPDG ) Right->SetPDGCode(sPDG);
2202 }
2203 {
2204 Left->SetPDGCode(SumPartonPDG(uPDG, sPDG));
2205 Right->SetPDGCode(SumPartonPDG(mPDG, nPDG));
2206 }
2207 Left->Set4Momentum(L4M+PL4M);
2208 Right->Set4Momentum(R4M+PR4M);
2209#ifdef debug
2210 G4cout<<"G4QIonIonCollision::Breeder:LL/RR s4M="<<fuString->Get4Momentum()
2211 <<",S="<<L4M+PL4M+R4M+PR4M<<", L="<<Left->Get4Momentum()<<", R="
2212 <<Right->Get4Momentum()<<G4endl;
2213#endif
2214 }
2215 else if(fusionDONE<0)
2216 {
2217 Left->SetPDGCode(SumPartonPDG(uPDG, nPDG));
2218 Left->Set4Momentum(L4M+PR4M);
2219 Right->SetPDGCode(SumPartonPDG(mPDG, sPDG));
2220 Right->Set4Momentum(R4M+PL4M);
2221#ifdef debug
2222 G4cout<<"G4QIonIonCollision::Breeder:LR/RL s4M="<<fuString->Get4Momentum()
2223 <<",S="<<L4M+PL4M+R4M+PR4M<<", L="<<Left->Get4Momentum()<<", R="
2224 <<Right->Get4Momentum()<<G4endl;
2225#endif
2226 }
2227#ifdef debug
2228 else G4cout<<"-Warning-G4QIonIonCollision::Breeder: WrongStringFusion"<<G4endl;
2229#endif
2230#ifdef edebug
2231 G4cout<<"#EMC#G4QIonIonCollision::Breed:StringFused,F="<<fusionDONE<<",L="<<L4M
2232 <<",R="<<R4M<<",pL="<<PL4M<<",pR="<<PR4M<<",nL="<<Left->Get4Momentum()
2233 <<",nR="<<Right->Get4Momentum()<<",S="<<fuString->Get4Momentum()<<G4endl;
2234#endif
2235 if(fusionDONE)
2236 {
2237#ifdef debug
2238 G4cout<<"###G4QIonIonCollision::Breed: Str#"<<astring<<" fused/w Str#"<<fustr
2239 <<"->"<<fuString->Get4Momentum()<<fuString->Get4Momentum().m()<<G4endl;
2240#endif
2241 continue; // @@ killing of the curString is excluded
2242 }
2243 }
2244#ifdef debug
2245 else
2246 {
2247
2248 G4cerr<<"**G4QIonIonCollision::Breed:*NoPart*M="<<curString->Get4Momentum().m()
2249 <<", F="<<fusionDONE<<", LPDG="<<curString->GetLeftParton()->GetPDGCode()
2250 <<", RPDG="<<curString->GetRightParton()->GetPDGCode()<<G4endl;
2251 }
2252#endif
2253 }
2254 if(!fusionDONE) // Fusion of theString failed, try hadrons
2255 {
2256 G4int nHadr=theResult->size(); // #of collected resulting hadrons upToNow
2257#ifdef debug
2258 G4cout<<"+++4QFragmentation::Breeder:*TryHadr* M="<<curString->Get4Momentum().m()
2259 <<", nH="<<nHadr<<", LPDG="<<curString->GetLeftParton()->GetPDGCode()
2260 <<", RPDG="<<curString->GetRightParton()->GetPDGCode()<<G4endl;
2261#endif
2262 // The only solution now is to try fusion with the secondary hadron
2263 while( (nHadr=theResult->size()) && !theHadrons)
2264 {
2265#ifdef edebug
2266 for(G4int i=0; i<nHadr; i++)
2267 {
2268 G4LorentzVector h4M=(*theResult)[i]->Get4Momentum();
2269 G4int hPDG=(*theResult)[i]->GetPDGCode();
2270 G4QContent hQC=(*theResult)[i]->GetQC();
2271 G4cout<<"-EMC-G4QInel::Breed:H#"<<i<<",4M="<<h4M<<hQC<<",PDG="<<hPDG<<G4endl;
2272 }
2273#endif
2274 G4int fusDONE=0; // String+Hadron fusion didn't happen
2275 G4int fuhad=-1; // The found hadron index
2276 G4int newPDG=0; // PDG ofTheParton afterMerging with Hadr
2277 G4int secPDG=0; // Second PDG oParton afterMerging w/Hadr
2278 G4double maM2=-DBL_MAX; // Prototype of the max ResultingStringM2
2279 G4LorentzVector selH4M(0.,0.,0.,0.); // 4-mom of the selected hadron
2280 G4QHadron* selHP=0; // Pointer to the used hadron for erasing
2281 G4QString* cString=strings[astring]; // Must be the last string by definition
2282 G4LorentzVector cString4M = cString->Get4Momentum();
2283 cLeft=cString->GetLeftParton();
2284 cRight=cString->GetRightParton();
2285 G4int sumT=cLeft->GetType()+cRight->GetType();
2286 sPDG=cLeft->GetPDGCode();
2287 nPDG=cRight->GetPDGCode();
2288 G4int cMax=0; // Both or only one cuark can merge
2289 for (G4int reh=0; reh < nHadr; reh++)
2290 {
2291 G4QHadron* curHadr=(*theResult)[reh];
2292 G4QContent curQC=curHadr->GetQC(); // Quark content of the hadron
2293 G4int partPDG1=curQC.AddParton(sPDG);
2294 G4int partPDG2=curQC.AddParton(nPDG);
2295#ifdef debug
2296 G4cout<<"G4QIonIonCollision::Breeder: Hadron # "<<reh<<", QC="<<curQC
2297 <<", P1="<<partPDG1<<", P2="<<partPDG2<<G4endl;
2298#endif
2299 if(partPDG1 || partPDG2) // Hadron can merge at least w/one parton
2300 {
2301 G4int cCur=1;
2302 if(sumT>3 && partPDG1 && partPDG2) cCur=2;
2303 G4LorentzVector curHadr4M = curHadr->Get4Momentum();
2304 G4double M2=(cString4M+curHadr4M).m2();// SquaredMass of theResultingString
2305#ifdef debug
2306 G4cout<<"G4QIonIonCollision::Breeder:*IN*Hadron#"<<reh<<",M2="<<M2<<G4endl;
2307#endif
2308 if( (sumT<4 || cCur>=cMax) && M2 > maM2) // DeepAnnihilation for DiQ-aDiQ
2309 {
2310 maM2=M2;
2311 fuhad=reh;
2312 if(partPDG1)
2313 {
2314 fusDONE= 1;
2315 newPDG=partPDG1;
2316 if(partPDG2)
2317 {
2318 secPDG=partPDG2;
2319 cMax=cCur;
2320 }
2321 }
2322 else
2323 {
2324 fusDONE=-1;
2325 newPDG=partPDG2;
2326 }
2327#ifdef debug
2328 G4cout<<"G4QInel::Br:*Selected*,P1="<<partPDG1<<",P2="<<partPDG2<<G4endl;
2329#endif
2330 selH4M=curHadr4M;
2331 selHP=curHadr;
2332 } // End of IF(update selection)
2333 } // End of IF(HadronCanMergeWithTheString)
2334 } // End of the LOOP over Hadrons
2335#ifdef debug
2336 G4cout<<"G4QIonIonCollision::Breeder: fuh="<<fuhad<<",fus="<<fusDONE<<G4endl;
2337#endif
2338 if(fuhad>-1) // The partner was found -> fuse H&S
2339 {
2340 if (fusDONE>0) // Fuse hadron with the left parton
2341 {
2342 cLeft->SetPDGCode(newPDG);
2343 G4LorentzVector newLeft=cLeft->Get4Momentum()+selH4M;
2344 cLeft->Set4Momentum(newLeft);
2345 if(secPDG && cMax>1)
2346 {
2347#ifdef debug
2348 G4cout<<"G4QInel::Br:TryReduce, nPDG="<<newPDG<<",sPDG="<<secPDG<<G4endl;
2349#endif
2350 cLeft->ReduceDiQADiQ(cLeft, cRight);
2351 }
2352#ifdef debug
2353 G4cout<<"G4QIonIonCollision::Breed: Left, s4M="<<curString->Get4Momentum()
2354 <<", L4M="<<newLeft<<", R4M="<<cRight->Get4Momentum()<<", h4M="
2355 <<selH4M<<", newL="<<newPDG<<", oldR="<<cRight->GetPDGCode()<<G4endl;
2356#endif
2357 }
2358 else if(fusDONE<0) // Fuse hadron with the right parton
2359 {
2360 cRight->SetPDGCode(newPDG);
2361 G4LorentzVector newRight=cRight->Get4Momentum()+selH4M;
2362 cRight->Set4Momentum(newRight);
2363#ifdef debug
2364 G4cout<<"G4QIonIonCollision::Breed: Right, s4M="<<curString->Get4Momentum()
2365 <<", L4M="<<cLeft->Get4Momentum()<<", R4M="<<newRight<<", h4M="
2366 <<selH4M<<", newR="<<newPDG<<", oldL="<<cLeft->GetPDGCode()<<G4endl;
2367#endif
2368 }
2369#ifdef debug
2370 else G4cout<<"-G4QIonIonCollision::Breed: Wrong String+HadronFusion"<<G4endl;
2371#endif
2372#ifdef debug
2373 if(fusDONE) G4cout<<"####G4QIonIonCollision::Breeder: String #"<<astring
2374 <<" is fused with Hadron #"<<fuhad
2375 <<", new4Mom="<<curString->Get4Momentum()
2376 <<", M2="<<curString->Get4Momentum().m2()
2377 <<", QC="<<curString->GetQC()<<G4endl;
2378#endif
2379 }
2380 else
2381 {
2382#ifdef debug
2383 G4cerr<<"**G4QIonIonCollision::Breed:*NoH*,M="<<curString->Get4Momentum().m()
2384 <<", LPDG="<<curString->GetLeftParton()->GetPDGCode()
2385 <<", RPDG="<<curString->GetRightParton()->GetPDGCode()<<G4endl;
2386 // @@ Temporary exception for the future solution
2387 //G4Exception("G4QIonIonCollision::Breed:","72",FatalException,"SHNotFused");
2388#endif
2389 break; // Breake the While LOOP
2390 } // End of the namespace where both Fusion and reduction have failed
2391 // The fused hadron must be excluded from theResult
2392#ifdef debug
2393 G4cout<<"G4QIonIonCollision::Breed: before HR, nH="<<theResult->size()<<G4endl;
2394 G4int icon=0; // Loop counter
2395#endif
2396 G4QHadronVector::iterator ih;
2397#ifdef debug
2398 G4bool found=false;
2399#endif
2400 for(ih = theResult->begin(); ih != theResult->end(); ih++)
2401 {
2402#ifdef debug
2403 G4cout<<"G4QInelast::Breeder:#"<<icon<<", i="<<(*ih)<<", sH="<<selHP<<G4endl;
2404 icon++;
2405#endif
2406 if((*ih)==selHP)
2407 {
2408#ifdef debug
2409 G4cout<<"G4QInel::Breed: *HadronFound*, PDG="<<selHP->GetPDGCode()<<G4endl;
2410#endif
2411 G4LorentzVector p4M=selHP->Get4Momentum();
2412 curString4M+=p4M;
2413#ifdef edebug
2414 G4int Chg=selHP->GetCharge();
2415 G4int BaN=selHP->GetBaryonNumber();
2416 curStrChg+=Chg;
2417 curStrBaN+=BaN;
2418 G4cout<<"-EMC->>>G4QIonIonCollision::Breed: S+=H, 4M="<<curString4M<<", M="
2419 <<curString4M.m()<<", Charge="<<curStrChg<<", B="<<curStrBaN<<G4endl;
2420#endif
2421 delete selHP; // delete the Hadron
2422 theResult->erase(ih); // erase the Hadron from theResult
2423#ifdef debug
2424 found=true;
2425#endif
2426 break; // beak the LOOP over hadrons
2427 }
2428 } // End of the LOOP over hadrons
2429#ifdef debug
2430 if(!found) G4cout<<"*G4QIonIonCollision::Breed:nH="<<theResult->size()<<G4endl;
2431#endif
2432 // New attempt of the string decay
2433 theHadrons=curString->FragmentString(true);//! Try to fragment the new String !
2434#ifdef debug
2435 G4cout<<"G4QInel::Breeder: tH="<<theHadrons<<",nH="<<theResult->size()<<G4endl;
2436#endif
2437 } // End of the while LOOP over the fusion with hadrons
2438#ifdef debug
2439 G4cout<<"*G4QIonIonCollision::Breed: *CanTryToDecay?* nH="<<theHadrons<<", next="
2440 <<next<<" =? nS="<<strings.size()<<", nR="<<theResult->size()<<G4endl;
2441#endif
2442 if(!theHadrons && next == strings.size() && !(theResult->size()))// TryToDecay
2443 {
2444 G4QContent miQC=curString->GetQC(); // QContent of the Lightest Hadron
2445 G4int miPDG=miQC.GetSPDGCode(); // PDG of the Lightest Hadron
2446 if(miPDG == 10) // ==> Decay Hadron-Chipolino
2447 {
2448 G4QChipolino QCh(miQC); // define theTotalResidual as aChipolino
2449 G4QPDGCode h1QPDG=QCh.GetQPDG1(); // QPDG of the first hadron
2450 G4double h1M =h1QPDG.GetMass();// Mass of the first hadron
2451 G4QPDGCode h2QPDG=QCh.GetQPDG2(); // QPDG of the second hadron
2452 G4double h2M =h2QPDG.GetMass();// Mass of the second hadron
2453 G4double ttM =curString4M.m(); // Real Mass of the Chipolino
2454 if(h1M+h2M<ttM+eps) // Two particles decay of Chipolino
2455 {
2456 G4LorentzVector h14M(0.,0.,0.,h1M);
2457 G4LorentzVector h24M(0.,0.,0.,h2M);
2458 if(std::fabs(ttM-h1M-h2M)<=eps)
2459 {
2460 G4double part1=h1M/(h1M+h2M);
2461 h14M=part1*curString4M;
2462 h24M=curString4M-h14M;
2463 }
2464 else
2465 {
2466 if(!G4QHadron(curString4M).DecayIn2(h14M,h24M))
2467 {
2468 G4cerr<<"***G4QIonIonCollision::Breeder: tM="<<ttM<<"->h1="<<h1QPDG
2469 <<"("<<h1M<<")+h2="<<h1QPDG<<"("<<h2M<<")="<<h1M+h2M<<G4endl;
2470 G4Exception("G4QIonIonCollision::Breeder:","72",FatalException,"QDec");
2471 }
2472 }
2473 G4QHadron* h1H = new G4QHadron(h1QPDG.GetPDGCode(),h14M);
2474 theResult->push_back(h1H); // (delete equivalent)
2475#ifdef debug
2476 G4LorentzVector f4M=h1H->Get4Momentum();
2477 G4int fPD=h1H->GetPDGCode();
2478 G4int fCg=h1H->GetCharge();
2479 G4int fBN=h1H->GetBaryonNumber();
2480 G4cout<<"-EMC->>G4QIonIonCollision::Breed:String=HadrChiPro1's filled,f4M="
2481 <<f4M<<", fPDG="<<fPD<<", fCg="<<fCg<<", fBN="<<fBN<<G4endl;
2482#endif
2483 G4QHadron* h2H = new G4QHadron(h2QPDG.GetPDGCode(),h24M);
2484 theResult->push_back(h2H); // (delete equivalent)
2485#ifdef debug
2486 G4LorentzVector s4M=h2H->Get4Momentum();
2487 G4int sPD=h2H->GetPDGCode();
2488 G4int sCg=h2H->GetCharge();
2489 G4int sBN=h2H->GetBaryonNumber();
2490 G4cout<<"-EMC->>G4QIonIonCollision::Breed:String=HadrChiPro2's filled,s4M="
2491 <<s4M<<", sPDG="<<sPD<<", sCg="<<sCg<<", sBN="<<sBN<<G4endl;
2492#endif
2493#ifdef edebug
2494 G4cout<<"-EMC-..Chi..G4QIonIonCollision::Breeder: DecayCHECK, Ch4M="
2495 <<curString4M<<", d4M="<<curString4M-h14M-h24M<<G4endl;
2496#endif
2497 break; // Go out of the main StringDecay LOOP
2498 }
2499 else
2500 {
2501 G4Quasmon* stringQuasmon = new G4Quasmon(miQC, curString4M);// String->Quas
2502 theQuasmons.push_back(stringQuasmon);
2503 break; // Go out of the main StringDecay LOOP
2504 }
2505 }
2506 else if(miPDG) // Decay Hadron as a Resonans
2507 {
2508 if (miPDG>0 && miPDG%10 < 3) miPDG+=2; // Convert to Resonans
2509 else if(miPDG<0 && (-miPDG)%10< 3) miPDG-=2; // Convert to antiResonans
2510 G4Quasmon Quasm;
2511 G4QHadron* sHad = new G4QHadron(miPDG,curString4M);
2512 G4QHadronVector* tmpQHadVec=Quasm.DecayQHadron(sHad); // It deletes sHad
2513 G4int tmpN=tmpQHadVec->size();
2514#ifdef debug
2515 G4cout<<"G4QIonIonCollision::Breeder: Decay the Last, Res#H="<<tmpN<<G4endl;
2516#endif
2517 if(tmpN>1)
2518 {
2519 for(G4int aH=0; aH < tmpN; aH++)
2520 {
2521 theResult->push_back((*tmpQHadVec)[aH]);//TheDecayProdOfHadronDecIsFilled
2522#ifdef debug
2523 G4QHadron* prodH =(*tmpQHadVec)[aH];
2524 G4LorentzVector p4M=prodH->Get4Momentum();
2525 G4int PDG=prodH->GetPDGCode();
2526 G4int Chg=prodH->GetCharge();
2527 G4int BaN=prodH->GetBaryonNumber();
2528 G4cout<<"-EMC->>G4QIonIonCollis::Breed:String=Hadr,H#"<<aH<<" filled,4M="
2529 <<p4M<<", PDG="<<PDG<<", Chg="<<Chg<<", BaN="<<BaN<<G4endl;
2530#endif
2531 }
2532 }
2533 else
2534 {
2535 G4Quasmon* stringQuasmon = new G4Quasmon(miQC, curString4M);// String->Quas
2536#ifdef debug
2537 G4cout<<"G4QIonIonCollision::Breeder:==> to Quasm="<<miQC<<curString4M
2538 <<", pNuc="<<theProjNucleus<<theProjNucleus.Get4Momentum()<<", tNuc="
2539 <<theTargNucleus<<theTargNucleus.Get4Momentum()<<", NString="
2540 <<strings.size()<<", nR="<<theResult->size()<<", nQ="
2541 <<theQuasmons.size()<<G4endl;
2542#endif
2543 theQuasmons.push_back(stringQuasmon);
2544 delete sHad;
2545 tmpQHadVec->clear();
2546 delete tmpQHadVec; // WhoCallsDecayQHadron is responsible for clear&delete
2547 break; // Go out of the main StringDecay LOOP
2548 }
2549 tmpQHadVec->clear();
2550 delete tmpQHadVec; // Who calls DecayQHadron is responsible for clear&delete
2551 break; // Go out of the main String Decay LOOP
2552 }
2553 } // End of the DecayOfTheLast
2554 } // End of IF(String-Hadron fusion)
2555 } // End of IF(NO_Hadrons) for String-String and String-Hadron fusion
2556 // The last hope is to CORREC the string, using other strings (ForwardInLOOP)
2557#ifdef debug
2558 G4cout<<"G4QIonIonCollision::Breeder: theH="<<theHadrons<<"?=0, next="<<next<<G4endl;
2559#endif
2560 if(!theHadrons && next < strings.size()) // ForwardInLOOP strings exist
2561 {
2562 // @@ string can be not convertable to one hadron (2,0.0,0,2,0) - To be improved
2563 G4QContent miQC=curString->GetQC(); // QContent of the Lightest Hadron
2564 G4int miPDG=miQC.GetSPDGCode();// PDG of the Lightest Hadron
2565#ifdef debug
2566 G4cout<<">>>G4QIonIonCollision::Breeder: SQC="<<miQC<<", miSPDG="<<miPDG<<G4endl;
2567#endif
2568 G4double miM=0.; // Prototype of the Mass of the Cur LightestHadron
2569 if(miPDG!=10) miM=G4QPDGCode(miPDG).GetMass(); // Mass of theCurLightestOneHadron
2570 else
2571 {
2572 G4QChipolino QCh(miQC); // define the TotalString as a Chipolino
2573 miM=QCh.GetQPDG1().GetMass()+QCh.GetQPDG2().GetMass();//MinMass of theStringChipo
2574 }
2575 G4double cM2=curString4M.m2(); // Actual squared mass of the Cur String
2576#ifdef debug
2577 G4cout<<">>>G4QIonIonCollision::Breeder: minMass="<<miM<<", realM2="<<cM2<<G4endl;
2578#endif
2579 G4double cM=0.;
2580 if(cM2>0.)
2581 {
2582 cM=std::sqrt(cM2);
2583 if(std::fabs(cM-miM) < eps) // Convert to hadron(2 hadrons) w/o calculations
2584 {
2585 if(miPDG!=10)
2586 {
2587 G4QHadron* sHad = new G4QHadron(miPDG,curString4M);
2588 theResult->push_back(sHad);// Fill the curString as a hadron
2589#ifdef debug
2590 G4cout<<">>>G4QIonIonCollision::Breeder:S->H="<<miPDG<<curString4M<<G4endl;
2591#endif
2592 }
2593 else
2594 {
2595 G4QChipolino QCh(miQC); // define TotalResidual as a Chipolino
2596 G4QPDGCode h1QPDG=QCh.GetQPDG1(); // QPDG of the first hadron
2597 G4double h1M =h1QPDG.GetMass(); // Mass of the first hadron
2598 G4QPDGCode h2QPDG=QCh.GetQPDG2(); // QPDG of the second hadron
2599 G4double h2M =h2QPDG.GetMass(); // Mass of the second hadron
2600 G4double pt1 =h1M/(h1M+h2M); // 4-mom part of the first hadron
2601 G4LorentzVector h14M=pt1*curString4M; // 4-mom of the first hadron
2602 G4LorentzVector h24M=curString4M-h14M;// 4-mom of the second hadron
2603 G4QHadron* h1H = new G4QHadron(h1QPDG.GetPDGCode(),h14M);
2604 theResult->push_back(h1H); // (delete equivalent)
2605#ifdef debug
2606 G4LorentzVector f4M=h1H->Get4Momentum();
2607 G4int fPD=h1H->GetPDGCode();
2608 G4int fCg=h1H->GetCharge();
2609 G4int fBN=h1H->GetBaryonNumber();
2610 G4cout<<"-EMC->>G4QIonIonCollision::Breed:Str=2HadrAR Prod-F is filled, f4M="
2611 <<f4M<<", fPDG="<<fPD<<", fCg="<<fCg<<", fBN="<<fBN<<G4endl;
2612#endif
2613 G4QHadron* h2H = new G4QHadron(h2QPDG.GetPDGCode(),h24M);
2614 theResult->push_back(h2H); // (delete equivalent)
2615#ifdef debug
2616 G4LorentzVector s4M=h2H->Get4Momentum();
2617 G4int sPD=h2H->GetPDGCode();
2618 G4int sCg=h2H->GetCharge();
2619 G4int sBN=h2H->GetBaryonNumber();
2620 G4cout<<"-EMC->>G4QIonIonCollision::Breed:Str=2HadrAR Prod-S is filled, s4M="
2621 <<s4M<<", sPDG="<<sPD<<", sCg="<<sCg<<", sBN="<<sBN<<G4endl;
2622#endif
2623 }
2624 continue; // Continue the LOOP over the curStrings
2625 }
2626 else // Try to recover (+/-) to min Mass
2627 {
2628 G4ThreeVector cP=curString4M.vect(); // Momentum of the curString
2629 G4double cE=curString4M.e(); // Energy of the curString
2630 G4ThreeVector curV=cP/cE; // curRelativeVelocity
2631 G4double miM2=miM*miM;
2632 G4int restr=0; // To use beyon the LOOP for printing
2633 G4int fustr=0; // Selected String-partner (0 = NotFound)
2634 G4double selX=0.; // Selected value of x
2635 G4double maD=-DBL_MAX; // Maximum Free Mass
2636 G4double Vmin=DBL_MAX; // min Velocity Distance
2637 G4LorentzVector s4M(0.,0.,0.,0.); // Selected 4-mom of the hadron
2638#ifdef debug
2639 G4cout<<"G4QIonIonCollision::Breeder: TryRecover, cV="<<curV<<G4endl;
2640#endif
2641 nOfStr=strings.size();
2642 for(restr=next; restr < nOfStr; ++restr) if(restr != astring)
2643 {
2644 G4QString* pString=strings[restr];
2645 G4LorentzVector p4M=pString->Get4Momentum();
2646 G4ThreeVector pP=p4M.vect(); // Momentum of the partnerString
2647 G4double pE=p4M.e(); // Energy of the partnerString
2648 G4double D2=cE*pE-cP.dot(pP);
2649 G4double pM2=p4M.m2();
2650 G4double dM4=pM2*(miM2-cM2);
2651 G4double D=D2*D2+dM4;
2652 G4double x=-1.; // Bad preexpectation
2653 if(D >= 0. && pM2>.01) x=(std::sqrt(D)-D2)/pM2; // what we should get from p
2654#ifdef debug
2655 else G4cout<<"G4QIonIonCollis::Breed:D="<<D<<",D2="<<D2<<",dM="<<dM4<<G4endl;
2656 G4cout<<"G4QIonIonCollision::Breed: pM2="<<pM2<<",D2="<<D2<<",x="<<x<<G4endl;
2657#endif
2658 if(x > 0. && x < 1.) // We are getting x part of p4M
2659 {
2660 G4QContent pQC=pString->GetQC(); // Quark Content of The Partner
2661 G4int pPDG=pQC.GetSPDGCode();// PDG of The Lightest Hadron for the Partner
2662 G4double pM=0.; // Mass of the LightestHadron
2663 if(pPDG==10)
2664 {
2665 G4QChipolino QCh(pQC); // define the TotalString as a Chipolino
2666 pM=QCh.GetQPDG1().GetMass()+QCh.GetQPDG2().GetMass();// Mass of Chipolino
2667 }
2668 else pM=G4QPDGCode(pPDG).GetMass();// Mass of theLightestHadron for Partner
2669 G4double rM=std::sqrt(pM2); // Real mass of the string-partner
2670 G4double delta=(1.-x)*rM-pM;// @@ Minimum CM disterbance measure
2671 if(delta > 0. && delta > maD)
2672 {
2673 maD=delta;
2674#ifdef debug
2675 G4cout<<"G4QIonIonCollision::Breed: Subtr,S#"<<restr<<",d="<<maD<<G4endl;
2676#endif
2677 fustr=restr;
2678 selX=x;
2679 s4M=p4M;
2680 }
2681 }
2682 else if(x <= 0.) // We are adding to p4M, so use RelVelocity
2683 {
2684 G4ThreeVector pV=pP/pE; // curRelativeVelocity
2685 G4double dV=(curV-pV).mag2();// SquaredDifferenceBetweenRelVelocities
2686 if(dV < Vmin)
2687 {
2688#ifdef debug
2689 G4cout<<"G4QIonIonCollision::Breed: FreeAdd,S#"<<restr<<",x="<<x<<G4endl;
2690#endif
2691 Vmin=dV;
2692 fustr=restr;
2693 selX=x;
2694 s4M=p4M;
2695 }
2696 }
2697#ifdef debug
2698 G4cout<<"G4QIonIonCollision::Breed:EndOfLOOP r="<<restr<<"<"<<nOfStr<<G4endl;
2699#endif
2700 } // End of the LOOP over string-partners for Correction
2701#ifdef debug
2702 G4cout<<"G4QIonIonCollision::Breeder: AfterLOOP fustr="<<fustr<<G4endl;
2703#endif
2704 if(fustr)
2705 {
2706#ifdef edebug
2707 G4LorentzVector sum4M=s4M+curString4M;
2708 G4cout<<"G4QIonIonCollision::Breeder: Found Sum4M="<<sum4M<<G4endl;
2709#endif
2710 G4QString* pString=strings[fustr];
2711 curString4M+=selX*s4M;
2712 if(std::abs(miPDG)%10 > 2) // Decay String-Hadron-Resonance
2713 {
2714 G4Quasmon Quasm;
2715 G4QHadron* sHad = new G4QHadron(miPDG,curString4M);
2716 G4QHadronVector* tmpQHadVec=Quasm.DecayQHadron(sHad); // It deletes sHad
2717#ifdef debug
2718 G4cout<<"G4QIonIonCollision::Breed:DecStH,nH="<<tmpQHadVec->size()<<G4endl;
2719#endif
2720 for(unsigned aH=0; aH < tmpQHadVec->size(); aH++)
2721 {
2722 theResult->push_back((*tmpQHadVec)[aH]);//TheDecayProdOfHadron is filled
2723#ifdef debug
2724 G4QHadron* prodH =(*tmpQHadVec)[aH];
2725 G4LorentzVector p4M=prodH->Get4Momentum();
2726 G4int PDG=prodH->GetPDGCode();
2727 G4int Chg=prodH->GetCharge();
2728 G4int BaN=prodH->GetBaryonNumber();
2729 G4cout<<"-EMC->>G4QIonIonCollision::Breed:St=Had,pH#"<<aH<<" filled, 4M="
2730 <<p4M<<", PDG="<<PDG<<", Chg="<<Chg<<", BaN="<<BaN<<G4endl;
2731#endif
2732 }
2733 tmpQHadVec->clear();
2734 delete tmpQHadVec; // Who calls DecayQHadron is responsibleRorClear&Delete
2735 }
2736 else if(miPDG == 10) // ==> Decay Hadron-Chipolino
2737 {
2738 G4QChipolino QCh(miQC); // define theTotalResid as aChipolino
2739 G4QPDGCode h1QPDG=QCh.GetQPDG1(); // QPDG of the first hadron
2740 G4double h1M =h1QPDG.GetMass();// Mass of the first hadron
2741 G4QPDGCode h2QPDG=QCh.GetQPDG2(); // QPDG of the second hadron
2742 G4double h2M =h2QPDG.GetMass();// Mass of the second hadron
2743 G4double ttM =curString4M.m(); // Real Mass of the Chipolino
2744 if(h1M+h2M<ttM+eps) // Two particles decay of Chipolino
2745 {
2746 G4LorentzVector h14M(0.,0.,0.,h1M);
2747 G4LorentzVector h24M(0.,0.,0.,h2M);
2748 if(std::fabs(ttM-h1M-h2M)<=eps)
2749 {
2750 G4double part1=h1M/(h1M+h2M);
2751 h14M=part1*curString4M;
2752 h24M=curString4M-h14M;
2753 }
2754 else
2755 {
2756 if(!G4QHadron(curString4M).DecayIn2(h14M,h24M))
2757 {
2758 G4cerr<<"***G4QIonIonCollision::Breeder: tM="<<ttM<<"->h1="<<h1QPDG
2759 <<"("<<h1M<<")+h2="<<h1QPDG<<"("<<h2M<<")="<<h1M+h2M<<G4endl;
2760 G4Exception("G4QIonIonCollision::Breeder:","72",FatalException,"CD");
2761 }
2762 }
2763 G4QHadron* h1H = new G4QHadron(h1QPDG.GetPDGCode(),h14M);
2764 theResult->push_back(h1H); // (delete equivalent)
2765#ifdef debug
2766 G4LorentzVector f4M=h1H->Get4Momentum();
2767 G4int fPD=h1H->GetPDGCode();
2768 G4int fCg=h1H->GetCharge();
2769 G4int fBN=h1H->GetBaryonNumber();
2770 G4cout<<"-EMC->>>G4QIonIonCollision::Breed:Str=Hadr Prod-F's filled,f4M="
2771 <<f4M<<", fPDG="<<fPD<<", fCg="<<fCg<<", fBN="<<fBN<<G4endl;
2772#endif
2773 G4QHadron* h2H = new G4QHadron(h2QPDG.GetPDGCode(),h24M);
2774 theResult->push_back(h2H); // (delete equivalent)
2775#ifdef debug
2776 G4LorentzVector s4M=h2H->Get4Momentum();
2777 G4int sPD=h2H->GetPDGCode();
2778 G4int sCg=h2H->GetCharge();
2779 G4int sBN=h2H->GetBaryonNumber();
2780 G4cout<<"-EMC->>>G4QIonIonCollision::Breed:Str=Hadr Prod-S's filled,s4M="
2781 <<s4M<<", sPDG="<<sPD<<", sCg="<<sCg<<", sBN="<<sBN<<G4endl;
2782#endif
2783#ifdef edebug
2784 G4cout<<"-EMC-Chipo.G4QIonIonCollision::Breed:DecCHECK,c4M="<<curString4M
2785 <<", ChQC="<<miQC<<", d4M="<<curString4M-h14M-h24M<<G4endl;
2786#endif
2787 }
2788 else
2789 {
2790 G4cerr<<"***G4QInel::Breeder: tM="<<ttM<<miQC<<"->h1="<<h1QPDG<<"(" <<h1M
2791 <<")+h2="<<h1QPDG<<"("<<h2M<<") = "<<h1M+h2M<<G4endl;
2792 G4Exception("G4QIonIonCollision::Breeder:","72",FatalException,"ChiDec");
2793 }
2794 }
2795 else
2796 {
2797 G4QHadron* sHad = new G4QHadron(miPDG,curString4M);
2798 theResult->push_back(sHad); // The original string-hadron is filled
2799#ifdef debug
2800 G4cout<<"-EMC->>>G4QIonIonCollision::Breeder:Str=Hadr Filled, 4M="
2801 <<curString4M<<", PDG="<<miPDG<<G4endl;
2802#endif
2803 }
2804 G4double corF=1-selX;
2805 G4QParton* Left=pString->GetLeftParton();
2806 G4QParton* Right=pString->GetRightParton();
2807 Left->Set4Momentum(corF*Left->Get4Momentum());
2808 Right->Set4Momentum(corF*Right->Get4Momentum());
2809#ifdef edebug
2810 G4cout<<"-EMC-...Cor...G4QIonIonCollision::Breeder:CorCHECK Sum="<<sum4M
2811 <<" =? "<<curString4M+pString->Get4Momentum()<<", M="<<miM<<" ?= "
2812 <<curString4M.m()<<G4endl;
2813#endif
2814#ifdef debug
2815 G4cout<<">>>G4QIonIonCollision::Breeder:*Corrected* String->Hadr="<<miPDG
2816 <<curString4M<<" by String #"<<fustr<<G4endl;
2817#endif
2818 continue; // Continue the LOOP over the curStrings
2819 } // End of Found combination for the String on string correction
2820 } // End of the Try-to-recover String+String correction algorithm
2821 } // End of IF(CM2>0.)
2822 } // End of IF(Can try to correct by String-String)
2823#ifdef debug
2824 else G4cerr<<"***G4QIonIonCollision::Breed:**No SSCorrection**, next="<<next<<G4endl;
2825#endif
2826 // ------------ At this point we can reduce the 3/-3 meson to 1/-1 meson ------------
2827 G4QParton* lpcS=curString->GetLeftParton();
2828 G4QParton* rpcS=curString->GetRightParton();
2829 G4int lPDGcS=lpcS->GetPDGCode();
2830 G4int rPDGcS=rpcS->GetPDGCode();
2831 if (lPDGcS==3 && rPDGcS==-3)
2832 {
2833 lpcS->SetPDGCode( 1);
2834 rpcS->SetPDGCode(-1);
2835 }
2836 else if(lPDGcS==-3 && rPDGcS==3)
2837 {
2838 lpcS->SetPDGCode(-1);
2839 rpcS->SetPDGCode( 1);
2840 }
2841 // -------- Now the only hope is Correction, using the already prodused Hadrons -----
2842 G4int nofRH=theResult->size(); // #of resulting Hadrons
2843#ifdef debug
2844 G4cout<<"G4QIonIonCollision::Breeder: theH="<<theHadrons<<", #OfH="<<nofRH<<G4endl;
2845#endif
2846 if(!theHadrons && nofRH) // Hadrons are existing for SH Correction
2847 {
2848#ifdef debug
2849 G4cerr<<"!G4QIonIonCollision::Breeder:CanTrySHCor, nH="<<theResult->size()<<G4endl;
2850#endif
2851 // @@ string can be not convertable to one hadron (2,0.0,0,2,0) - To be improved
2852 G4QContent miQC=curString->GetQC(); // QContent of the Lightest Hadron
2853 G4int miPDG=miQC.GetSPDGCode(); // PDG of the Lightest Hadron
2854 G4double miM=0.; // Prototype ofMass of theCurLightestHadron
2855 if(miPDG==10) // Mass of the Cur Lightest ChipolinoHadron
2856 {
2857 G4QChipolino QCh(miQC); // define the TotalString as a Chipolino
2858 miM=QCh.GetQPDG1().GetMass()+QCh.GetQPDG2().GetMass();//MinMass of theStringChipo
2859 }
2860 else miM=G4QPDGCode(miPDG).GetMass(); // Mass of the Cur Lightest Hadron
2861 G4double spM=0.; // Mass of the selected Hadron-Partner
2862 G4ThreeVector cP=curString4M.vect(); // Momentum of the curString
2863 G4double cE=curString4M.e(); // Energy of the curString
2864 G4ThreeVector curV=cP/cE; // curRelativeVelocity
2865 G4int reha=0; // Hadron # to use beyon the LOOP
2866 G4int fuha=0; // Selected Hadron-partner (0 = NotFound)
2867 G4double dMmin=DBL_MAX; // min Excess of the mass
2868 G4LorentzVector s4M(0.,0.,0.,0.); // Selected 4-mom of the Hadron+String
2869 G4double sM=0.; // Selected Mass of the Hadron+String
2870 for (reha=next; reha < nofRH; reha++) // LOOP over already collected hadrons
2871 {
2872 G4QHadron* pHadron=(*theResult)[reha];// Pointer to the current Partner-Hadron
2873 G4LorentzVector p4M=pHadron->Get4Momentum();
2874 G4double pM=p4M.m(); // Mass of the Partner-Hadron
2875 G4LorentzVector t4M=p4M+curString4M; // Total momentum of the compound
2876 G4double tM2=t4M.m2(); // Squared total mass of the compound
2877 if(tM2 >= sqr(pM+miM+eps)) // Condition of possible correction
2878 {
2879 G4double tM=std::sqrt(tM2); // Mass of the Hadron+String compound
2880 G4double dM=tM-pM-miM; // Excess of the compound mass
2881 if(dM < dMmin)
2882 {
2883 dMmin=dM;
2884 fuha=reha;
2885 spM=pM;
2886 s4M=t4M;
2887 sM=tM;
2888 }
2889 }
2890 } // End of the LOOP over string-partners for Correction
2891 if(fuha) // The hadron-partner was found
2892 {
2893 G4QHadron* pHadron=(*theResult)[fuha];// Necessary for update
2894 G4LorentzVector mi4M(0.,0.,0.,miM); // Prototype of the new String=Hadron
2895 if(miM+spM<sM+eps) // Decay into CorrectedString+theSameHadron
2896 {
2897 G4LorentzVector sp4M(0.,0.,0.,spM);
2898 if(std::fabs(sM-miM-spM)<=eps)
2899 {
2900 G4double part1=miM/(miM+spM);
2901 mi4M=part1*s4M;
2902 sp4M=s4M-mi4M;
2903 }
2904 else
2905 {
2906 if(!G4QHadron(s4M).DecayIn2(mi4M,sp4M))
2907 {
2908 G4cerr<<"***G4QIonIonCollision::Breeder: *SH*, tM="<<sM<<"->h1=("<<miPDG
2909 <<")"<<miM<<" + h2="<<spM<<" = "<<miM+spM<<G4endl;
2910 G4Exception("G4QIonIonCollision::Breeder:","72",FatalException,"SHChiDec");
2911 }
2912 }
2913 pHadron->Set4Momentum(sp4M);
2914#ifdef debug
2915 G4cout<<"-EMC->...>G4QIonIonCollision::Breed: H# "<<fuha<<" is updated, new4M="
2916 <<sp4M<<G4endl;
2917#endif
2918 }
2919 else
2920 {
2921 G4cerr<<"***G4QInel::Breeder: HS Failed, tM="<<sM<<"->h1M=("<<miPDG<<")"<<miM
2922 <<"+h2M="<<spM<<" = "<<miM+spM<<G4endl;
2923 G4Exception("G4QIonIonCollision::Breeder:","72",FatalException,"HSChipoliDec");
2924 }
2925 if(std::abs(miPDG)%10 > 2) // Decay Hadron-Resonans
2926 {
2927 G4Quasmon Quasm;
2928 G4QHadron* sHad = new G4QHadron(miPDG,mi4M);
2929 G4QHadronVector* tmpQHadVec=Quasm.DecayQHadron(sHad); // It deletes sHad
2930#ifdef debug
2931 G4cout<<"G4QInelast::Breeder: *HS* DecStrHad, nH="<<tmpQHadVec->size()<<G4endl;
2932#endif
2933 for(unsigned aH=0; aH < tmpQHadVec->size(); aH++)
2934 {
2935 theResult->push_back((*tmpQHadVec)[aH]);// TheDecayProductOfTheHadronIsFilled
2936#ifdef debug
2937 G4QHadron* prodH =(*tmpQHadVec)[aH];
2938 G4LorentzVector p4M=prodH->Get4Momentum();
2939 G4int PDG=prodH->GetPDGCode();
2940 G4int Chg=prodH->GetCharge();
2941 G4int BaN=prodH->GetBaryonNumber();
2942 G4cout<<"-EMC->>>G4QIonIonCollision::Breed:Str+Hadr PrH#"<<aH<<" filled, 4M="
2943 <<p4M<<", PDG="<<PDG<<", Chg="<<Chg<<", BaN="<<BaN<<G4endl;
2944#endif
2945 }
2946 tmpQHadVec->clear();
2947 delete tmpQHadVec; // Who calls DecayQHadron is responsible for clear & delete
2948 }
2949 else if(miPDG == 10) // ==> Decay Hadron-Chipolino
2950 {
2951 G4QChipolino QCh(miQC); // define the TotalResidual as a Chipolino
2952 G4QPDGCode h1QPDG=QCh.GetQPDG1(); // QPDG of the first hadron
2953 G4double h1M =h1QPDG.GetMass();// Mass of the first hadron
2954 G4QPDGCode h2QPDG=QCh.GetQPDG2(); // QPDG of the second hadron
2955 G4double h2M =h2QPDG.GetMass();// Mass of the second hadron
2956 G4double ttM =curString4M.m(); // Real Mass of the Chipolino
2957 if(h1M+h2M<miM+eps) // Two particles decay of Chipolino
2958 {
2959 G4LorentzVector h14M(0.,0.,0.,h1M);
2960 G4LorentzVector h24M(0.,0.,0.,h2M);
2961 if(std::fabs(ttM-h1M-h2M)<=eps)
2962 {
2963 G4double part1=h1M/(h1M+h2M);
2964 h14M=part1*mi4M;
2965 h24M=mi4M-h14M;
2966 }
2967 else
2968 {
2969 if(!G4QHadron(mi4M).DecayIn2(h14M,h24M))
2970 {
2971 G4cerr<<"***G4QIonIonCollision::Breeder: HS tM="<<ttM<<"->h1="<<h1QPDG
2972 <<"("<<h1M<<")+h2="<<h1QPDG<<"("<<h2M<<")="<<h1M+h2M<<G4endl;
2973 G4Exception("G4QIonIonCollision::Breeder:","72",FatalException,"ChiDec");
2974 }
2975 }
2976 G4QHadron* h1H = new G4QHadron(h1QPDG.GetPDGCode(),h14M);
2977 theResult->push_back(h1H); // (delete equivalent)
2978#ifdef debug
2979 G4LorentzVector f4M=h1H->Get4Momentum();
2980 G4int fPD=h1H->GetPDGCode();
2981 G4int fCg=h1H->GetCharge();
2982 G4int fBN=h1H->GetBaryonNumber();
2983 G4cout<<"-EMC->>G4QIonIonCollision::Breed: CorStrHadr Prod-1 is filled, f4M="
2984 <<f4M<<", fPDG="<<fPD<<", fCg="<<fCg<<", fBN="<<fBN<<G4endl;
2985#endif
2986 G4QHadron* h2H = new G4QHadron(h2QPDG.GetPDGCode(),h24M);
2987 theResult->push_back(h2H); // (delete equivalent)
2988#ifdef debug
2989 G4LorentzVector n4M=h2H->Get4Momentum();
2990 G4int nPD=h2H->GetPDGCode();
2991 G4int nCg=h2H->GetCharge();
2992 G4int nBN=h2H->GetBaryonNumber();
2993 G4cout<<"-EMC->>>G4QIonIonCollision::Breed:CorStrHadr Prod-2 is filled, n4M="
2994 <<n4M<<", nPDG="<<nPD<<", nCg="<<nCg<<", nBN="<<nBN<<G4endl;
2995#endif
2996#ifdef edebug
2997 G4cout<<"-EMC-...HS-Chipo...G4QIonIonCollision::Breeder:DecCHECK, Ch4M="
2998 <<mi4M<<", ChQC="<<miQC<<", d4M="<<mi4M-h14M-h24M<<G4endl;
2999#endif
3000 }
3001 }
3002 else
3003 {
3004 G4QHadron* sHad = new G4QHadron(miPDG, mi4M);
3005 theResult->push_back(sHad); // The original string=hadron is filled
3006#ifdef debug
3007 G4cout<<">>..>>G4QIonIonCollision::Breeder: CorStr=Hadr is Filled, 4M="
3008 <<curString4M<<", StPDG="<<miPDG<<G4endl;
3009#endif
3010 }
3011#ifdef edebug
3012 G4cout<<"-EMC-...Cor...G4QIonIonCollision::Breeder:StHadCor CHECK Sum="<<s4M
3013 <<" =? "<<mi4M+pHadron->Get4Momentum()<<G4endl;
3014#endif
3015#ifdef debug
3016 G4cout<<">>>G4QIonIonCollision::Breeder:*Corrected* String+Hadr="<<miPDG
3017 <<mi4M<<" by Hadron #"<<reha<<G4endl;
3018#endif
3019 continue; // Continue the LOOP over the curStrings
3020 }
3021 else
3022 {
3023#ifdef debug
3024 G4cout<<"-EMC->>>G4QIonIonCollision::Breeder: Str+Hadr Failed, 4M="<<curString4M
3025 <<", PDG="<<miPDG<<G4endl;
3026#endif
3027 }
3028 // @@@ convert string to Quasmon with curString4M
3029 G4QContent curStringQC=curString->GetQC();
3030 G4Quasmon* stringQuasmon = new G4Quasmon(curStringQC, curString4M);
3031 theQuasmons.push_back(stringQuasmon);
3032 continue; // Continue the LOOP over the curStrings
3033 } // End of IF(Can try the String-Hadron correction
3034 } // End of IF(NO_Hadrons) = Problem solution namespace
3035 G4Quasmon tmpQ; // @@ an issue of Q to decay resonances
3036 G4int nHfin=0;
3037 if(theHadrons) nHfin=theHadrons->size();
3038 else // !! Sum Up all strings and convert them in a Quasmon (Exception for development)
3039 {
3040 G4LorentzVector ss4M(0.,0.,0.,0.);
3041 G4QContent ssQC(0,0,0,0,0,0);
3042 G4int tnSt=strings.size();
3043 for(G4int i=astring; i < tnSt; ++i)
3044 {
3045 G4LorentzVector pS4M=strings[i]->Get4Momentum(); // String 4-momentum
3046 ss4M+=pS4M;
3047 G4QContent pSQC=strings[i]->GetQC(); // String Quark Content
3048 ssQC+=pSQC;
3049#ifdef debug
3050 G4cout<<"=--=>G4QIonIonCollision::Breed:S#"<<i<<",4M="<<pS4M<<",QC="<<pSQC<<G4endl;
3051#endif
3052 }
3053#ifdef debug
3054 G4cout<<"==>G4QIonIonCollision::Breed:AllStrings are summed up in a Quasmon"<<G4endl;
3055#endif
3056 G4Quasmon* stringQuasmon = new G4Quasmon(ssQC, ss4M);
3057 theQuasmons.push_back(stringQuasmon);
3058 break; // break the LOOP over Strings
3059 }
3060#ifdef debug
3061 G4cout<<"G4QIonIonCollision::Breeder: Trying to decay hadrons #ofHRes="<<nHfin<<G4endl;
3062#endif
3063 for(G4int aTrack=0; aTrack<nHfin; aTrack++)
3064 {
3065 G4QHadron* curHadron=(*theHadrons)[aTrack];
3066 G4int hPDG=curHadron->GetPDGCode();
3067#ifdef edebug
3068 G4LorentzVector curH4M=curHadron->Get4Momentum();
3069 G4int curHCh=curHadron->GetCharge();
3070 G4int curHBN=curHadron->GetBaryonNumber();
3071#endif
3072#ifdef debug
3073 G4cout<<">>...>>G4QIonIonCollision::Breeder:S#"<<astring<<",H#"<<aTrack<<",PDG="
3074 <<hPDG<<",4M="<<curHadron->Get4Momentum()<<G4endl;
3075#endif
3076 if(std::abs(hPDG)%10 > 2)
3077 {
3078 G4QHadronVector* tmpQHadVec=tmpQ.DecayQHadron(curHadron); // It deletes curHadron
3079#ifdef debug
3080 G4cout<<"G4QIonIonCollision::Breed:-DECAY'S DONE-,nH="<<tmpQHadVec->size()<<G4endl;
3081#endif
3082 //G4int tmpS=tmpQHadVec->size(); // "The elegant method" (tested) is commented
3083 //theResult->resize(tmpS+theResult->size()); // Resize theQHadrons length
3084 //copy(tmpQHadVec->begin(), tmpQHadVec->end(), theResult->end()-tmpS);
3085 for(unsigned aH=0; aH < tmpQHadVec->size(); aH++)
3086 {
3087 theResult->push_back((*tmpQHadVec)[aH]);// TheDecayProduct of TheHadron is filled
3088#ifdef edebug
3089 G4QHadron* prodH =(*tmpQHadVec)[aH];
3090 G4LorentzVector p4M=prodH->Get4Momentum();
3091 G4int PDG=prodH->GetPDGCode();
3092 G4int Chg=prodH->GetCharge();
3093 G4int BaN=prodH->GetBaryonNumber();
3094 curH4M-=p4M;
3095 curString4M-=p4M;
3096 curStrChg-=Chg;
3097 curStrBaN-=BaN;
3098 curHCh-=Chg;
3099 curHBN-=BaN;
3100 G4cout<<"-EMC->.>>G4QIonIonCollision::Breed:Str*Filled, 4M="<<p4M<<", PDG="<<PDG
3101 <<", Chg="<<Chg<<", BaN="<<BaN<<G4endl;
3102#endif
3103 }
3104#ifdef edebug
3105 G4cout<<"-EMC-.G4QIn::Br:Dec,r4M="<<curH4M<<",rC="<<curHCh<<",rB="<<curHBN<<G4endl;
3106#endif
3107 tmpQHadVec->clear();
3108 delete tmpQHadVec; // Who calls DecayQHadron is responsible for clear & delete
3109 }
3110 else // Chipolino is not checked here
3111 {
3112 theResult->push_back(curHadron); // The original hadron is filled
3113#ifdef edebug
3114 curString4M-=curH4M;
3115 G4int curCh=curHadron->GetCharge();
3116 G4int curBN=curHadron->GetBaryonNumber();
3117 curStrChg-=curCh;
3118 curStrBaN-=curBN;
3119 G4cout<<"-EMC->>..>>G4QIonIonCollision::Breeder: curH filled 4M="<<curH4M<<",PDG="
3120 <<curHadron->GetPDGCode()<<", Chg="<<curCh<<", BaN="<<curBN<<G4endl;
3121#endif
3122 }
3123 }
3124 // clean up (the issues are filled to theResult)
3125 if(theHadrons) delete theHadrons;
3126#ifdef edebug
3127 G4cout<<"-EMC-.......G4QIonIonCollision::Breeder: StringDecay CHECK, r4M="<<curString4M
3128 <<", rChg="<<curStrChg<<", rBaN="<<curStrBaN<<G4endl;
3129#endif
3130 } // End of the main LOOP over decaying strings
3131 G4LorentzVector rp4M=theProjNucleus.Get4Momentum(); // projNucleus 4-momentum in LS
3132 G4int rpPDG=theProjNucleus.GetPDG();
3133 G4QHadron* resPNuc = new G4QHadron(rpPDG,rp4M);
3134 theResult->push_back(resPNuc); // Fill the residual nucleus
3135 G4LorentzVector rt4M=theTargNucleus.Get4Momentum(); // Nucleus 4-momentum in LS
3136 G4int rtPDG=theTargNucleus.GetPDG();
3137 G4QHadron* resTNuc = new G4QHadron(rtPDG,rt4M);
3138 theResult->push_back(resTNuc); // Fill the residual nucleus
3139#ifdef edebug
3140 G4LorentzVector s4M(0.,0.,0.,0.); // Sum of the Result in LS
3141 G4int rCh=totChg;
3142 G4int rBN=totBaN;
3143 G4int nHadr=theResult->size();
3144 G4int nQuasm=theQuasmons.size();
3145 G4cout<<"-EMCLS-G4QInel::Breeder: #ofHadr="<<nHadr<<", #OfQ="<<nQuasm<<", rPA="<<rp4M.m()
3146 <<"="<<G4QNucleus(rpPDG).GetGSMass()<<", rTA="<<rt4M.m()<<"="
3147 <<G4QNucleus(rtPDG).GetGSMass()<<G4endl;
3148 for(G4int i=0; i<nHadr; i++)
3149 {
3150 G4LorentzVector hI4M=(*theResult)[i]->Get4Momentum();
3151 s4M+=hI4M;
3152 G4int hChg=(*theResult)[i]->GetCharge();
3153 rCh-=hChg;
3154 G4int hBaN=(*theResult)[i]->GetBaryonNumber();
3155 rBN-=hBaN;
3156 G4cout<<"-EMCLS-G4QIonIonCollision::Breeder: Hadron#"<<i<<", 4M="<<hI4M<<", PDG="
3157 <<(*theResult)[i]->GetPDGCode()<<", C="<<hChg<<", B="<<hBaN<<G4endl;
3158 }
3159 for(G4int i=0; i<nQuasm; i++)
3160 {
3161 G4LorentzVector hI4M=theQuasmons[i]->Get4Momentum();
3162 s4M+=hI4M;
3163 G4int hChg=theQuasmons[i]->GetCharge();
3164 rCh-=hChg;
3165 G4int hBaN=theQuasmons[i]->GetBaryonNumber();
3166 rBN-=hBaN;
3167 G4cout<<"-EMCLS-G4QIonIonCollision::Breeder: Quasmon#"<<i<<", 4M="<<hI4M<<", C="<<hChg
3168 <<", B="<<hBaN<<G4endl;
3169 }
3170 G4cout<<"-EMCLS-G4QInel::Breed: LS r4M="<<s4M-totLS4M<<",rC="<<rCh<<",rB="<<rBN<<G4endl;
3171#endif
3172 // Now we need to coolect particles for creation of a Quasmon @@ improve !!
3173 G4int nRes=theResult->size();
3174 if(nRes>2)
3175 {
3176 G4QHadronVector::iterator ih;
3177 G4QHadronVector::iterator nih;
3178 G4QHadronVector::iterator mih;
3179 G4QHadronVector::iterator lst=theResult->end();
3180 lst--;
3181 G4double minMesEn=DBL_MAX;
3182 G4double minBarEn=DBL_MAX;
3183 G4bool nfound=false;
3184 G4bool mfound=false;
3185 for(ih = theResult->begin(); ih < theResult->end(); ++ih) if(ih != lst)
3186 {
3187 G4LorentzVector h4M=(*ih)->Get4Momentum();
3188 G4int hPDG=(*ih)->GetPDGCode();
3189#ifdef debug
3190 G4cout<<"%->G4QIonIonCollision::Breeder: TRY hPDG="<<hPDG<<", h4M="<<h4M<<G4endl;
3191#endif
3192 if(hPDG>1111 && hPDG<3333)
3193 {
3194 G4double bE=h4M.e()-(*ih)->GetMass();
3195 if(bE < minBarEn)
3196 {
3197 minBarEn=bE;
3198 nih=ih;
3199 nfound=true;
3200 }
3201 }
3202 else if(hPDG>-1111)
3203 {
3204 G4double mE=h4M.e();
3205 if(mE < minMesEn)
3206 {
3207 minMesEn=mE;
3208 mih=ih;
3209 mfound=true;
3210 }
3211 }
3212 }
3213 if(nfound && mfound && minMesEn+minBarEn < maxEn)
3214 {
3215 G4QHadron* resNuc = (*theResult)[nRes-1]; // ResidualNucleus is theLastH
3216 theResult->pop_back(); // Fill the QHad-nucleus later
3217 G4LorentzVector q4M=(*nih)->Get4Momentum()+(*mih)->Get4Momentum();
3218 G4QContent qQC=(*nih)->GetQC()+(*mih)->GetQC();
3219 maxEn -= minBarEn+minMesEn; // Reduce the energy limit
3220#ifdef debug
3221 G4cout<<"%->G4QIonIonCollision::Breeder:Exclude,4M="<<(*nih)->Get4Momentum()
3222 <<", & 4m="<<(*mih)->Get4Momentum()<<G4endl;
3223#endif
3224 delete (*nih);
3225 delete (*mih);
3226 if(nih > mih) // Exclude nucleon first
3227 {
3228 theResult->erase(nih); // erase Baryon from theResult
3229 theResult->erase(mih); // erase Meson from theResult
3230 }
3231 else // Exclude meson first
3232 {
3233 theResult->erase(mih); // erase Baryon from theResult
3234 theResult->erase(nih); // erase Meson from theResult
3235 }
3236#ifdef debug
3237 G4cout<<"%->G4QI::Breed: BeforeLOOP, dE="<<maxEn<<", nR="<<theResult->size()<<G4endl;
3238#endif
3239 if(maxEn > 0.) // More hadrons to absorbe
3240 {
3241 for(ih = theResult->begin(); ih < theResult->end(); ih++)
3242 {
3243 G4LorentzVector h4M=(*ih)->Get4Momentum();
3244 G4int hPDG=(*ih)->GetPDGCode();
3245 G4double dE=0.;
3246 if (hPDG> 1111 && hPDG<3333) dE=h4M.e()-(*ih)->GetMass(); // Baryons
3247 else if(hPDG>-1111) dE=h4M.e(); // Mesons
3248 else continue; // Don't consider anti-baryons
3249 if(dE < maxEn)
3250 {
3251 maxEn-=dE;
3252 q4M+=h4M;
3253 qQC+=(*ih)->GetQC();
3254#ifdef debug
3255 G4cout<<"%->G4QIonIonCollision::Breed:Exclude,4M="<<h4M<<",dE="<<maxEn<<G4endl;
3256#endif
3257 delete (*ih);
3258 theResult->erase(ih); // erase Hadron from theResult
3259 --ih;
3260 }
3261 }
3262 }
3263 G4Quasmon* softQuasmon = new G4Quasmon(qQC, q4M); // SoftPart -> Quasmon
3264#ifdef debug
3265 G4cout<<"%->G4QIonIonCollision::Breed:QuasmonIsFilled,4M="<<q4M<<",QC="<<qQC<<G4endl;
3266#endif
3267 theQuasmons.push_back(softQuasmon);
3268 theResult->push_back(resNuc);
3269 }
3270#ifdef edebug
3271 G4LorentzVector f4M(0.,0.,0.,0.); // Sum of the Result in LS
3272 G4int fCh=totChg;
3273 G4int fBN=totBaN;
3274 G4int nHd=theResult->size();
3275 G4int nQm=theQuasmons.size();
3276 G4cout<<"-EMCLS-G4QIonIonCollision::Breeder:#ofHadr="<<nHd<<", #OfQ="<<nQm
3277 <<",rPA="<<rt4M.m()<<"="<<G4QNucleus(rtPDG).GetGSMass()
3278 <<",rTA="<<rt4M.m()<<"="<<G4QNucleus(rtPDG).GetGSMass()<<G4endl;
3279 for(G4int i=0; i<nHd; i++)
3280 {
3281 G4LorentzVector hI4M=(*theResult)[i]->Get4Momentum();
3282 f4M+=hI4M;
3283 G4int hChg=(*theResult)[i]->GetCharge();
3284 fCh-=hChg;
3285 G4int hBaN=(*theResult)[i]->GetBaryonNumber();
3286 fBN-=hBaN;
3287 G4cout<<"-EMCLS-G4QIonIonCollision::Breeder: Hadron#"<<i<<", 4M="<<hI4M<<", PDG="
3288 <<(*theResult)[i]->GetPDGCode()<<", C="<<hChg<<", B="<<hBaN<<G4endl;
3289 }
3290 for(G4int i=0; i<nQm; i++)
3291 {
3292 G4LorentzVector hI4M=theQuasmons[i]->Get4Momentum();
3293 f4M+=hI4M;
3294 G4int hChg=theQuasmons[i]->GetCharge();
3295 fCh-=hChg;
3296 G4int hBaN=theQuasmons[i]->GetBaryonNumber();
3297 fBN-=hBaN;
3298 G4cout<<"-EMCLS-G4QIonIonCollision::Breeder: Quasmon#"<<i<<", 4M="<<hI4M<<", C="
3299 <<hChg<<", B="<<hBaN<<G4endl;
3300 }
3301 G4cout<<"-EMCLS-G4QInel::Breed:, r4M="<<f4M-totLS4M<<", rC="<<fCh<<",rB="<<fBN<<G4endl;
3302#endif
3303 } // End of the soft Quasmon Creation
3304 return;
3305} // End of Breeder
CLHEP::Hep3Vector::dot
double dot(const Hep3Vector &) const
CLHEP::HepLorentzVector::vect
Hep3Vector vect() const
G4QContent::GetSPDGCode
G4int GetSPDGCode() const
Definition: G4QContent.cc:1204
G4QContent::AddParton
G4int AddParton(G4int pPDG) const
Definition: G4QContent.cc:1655
G4QHadron::GetCharge
G4int GetCharge() const
Definition: G4QHadron.hh:179
G4QHadron::DecayIn2
G4bool DecayIn2(G4LorentzVector &f4Mom, G4LorentzVector &s4Mom)
Definition: G4QHadron.cc:544
G4QIonIonCollision::AnnihilationOrder
G4int AnnihilationOrder(G4int LS, G4int MS, G4int uP, G4int mP, G4int sP, G4int nP)
Definition: G4QIonIonCollision.cc:3671
G4QIonIonCollision::ReducePair
std::pair< G4int, G4int > ReducePair(G4int P1, G4int P2) const
Definition: G4QIonIonCollision.cc:3651
G4QIonIonCollision::SumPartonPDG
G4int SumPartonPDG(G4int PDG1, G4int PFG2) const
Definition: G4QIonIonCollision.cc:3586
G4QNucleus::GetPDG
G4int GetPDG() const
Definition: G4QNucleus.hh:69
G4QPDGCode::GetPDGCode
G4int GetPDGCode() const
Definition: G4QPDGCode.hh:326
G4QString::FragmentString
G4QHadronVector * FragmentString(G4bool QL)
Definition: G4QString.cc:348
G4QString::GetDirection
G4int GetDirection() const
Definition: G4QString.hh:89
G4QString::GetQC
G4QContent GetQC() const
Definition: G4QString.hh:91
G4Quasmon::DecayQHadron
G4QHadronVector * DecayQHadron(G4QHadron *hadron)
Definition: G4Quasmon.cc:5779
CLHEP::swap
void swap(shared_ptr< P > &, shared_ptr< P > &)
Definition: memory.h:1247

Referenced by Fragment().

◆ ChooseX()

G4double G4QIonIonCollision::ChooseX ( G4double  Xmin,
G4double  Xmax 
) const
protected

Definition at line 3551 of file G4QIonIonCollision.cc.

3552{
3553 // choose an x between Xmin and Xmax with P(x) ~ 1/x @@ M.K. -> 1/sqrt(x)
3554 //G4double range=Xmax-Xmin;
3555 if(Xmax == Xmin) return Xmin;
3556 if( Xmin < 0. || Xmax < Xmin)
3557 {
3558 G4cerr<<"***G4QIonIonCollision::ChooseX: Xmin="<<Xmin<<", Xmax="<<Xmax<< G4endl;
3559 G4Exception("G4QIonIonCollision::ChooseX:","72",FatalException,"Bad X or X-Range");
3560 }
3561 //G4double x;
3562 //do {x=Xmin+G4UniformRand()*range;} while ( Xmin/x < G4UniformRand() );
3563 G4double sxi=std::sqrt(Xmin);
3564 G4double x=sqr(sxi+G4UniformRand()*(std::sqrt(Xmax)-sxi));
3565#ifdef debug
3566 G4cout<<"G4QIonIonCollision::ChooseX: DiffractiveX="<<x<<G4endl;
3567#endif
3568 return x;
3569} // End of ChooseX

Referenced by ExciteDiffParticipants(), and ExciteSingDiffParticipants().

◆ ExciteDiffParticipants()

G4bool G4QIonIonCollision::ExciteDiffParticipants ( G4QHadron *  aPartner,
G4QHadron *  bPartner 
) const
protected

Definition at line 3308 of file G4QIonIonCollision.cc.

3310{
3311 G4LorentzVector Pprojectile=projectile->Get4Momentum();
3312 G4double Mprojectile=projectile->GetMass();
3313 G4double Mprojectile2=Mprojectile*Mprojectile;
3314 G4LorentzVector Ptarget=target->Get4Momentum();
3315 G4double Mtarget=target->GetMass();
3316 G4double Mtarget2=Mtarget*Mtarget;
3317#ifdef debug
3318 G4cout<<"G4QInel::ExciteDiffPartici: Ep="<<Pprojectile.e()<<", Et="<<Ptarget.e()<<G4endl;
3319#endif
3320 // Transform momenta to cms and then rotate parallel to z axis;
3321 G4LorentzVector Psum=Pprojectile+Ptarget;
3322 G4LorentzRotation toCms(-Psum.boostVector()); // Boost Rotation to CMS
3323 G4LorentzVector Ptmp=toCms*Pprojectile;
3324 if(Ptmp.pz()<=0.) // "String" moving backwards in CMS, abort collision !! ? M.K.
3325 {
3326#ifdef debug
3327 G4cout<<"G4QIonIonCollision::ExciteDiffParticipants:*1* abort Collision!! *1*"<<G4endl;
3328#endif
3329 return false;
3330 }
3331 toCms.rotateZ(-Ptmp.phi());
3332 toCms.rotateY(-Ptmp.theta());
3333#ifdef debug
3334 G4cout<<"G4QIonIonCollision::ExciteDiffParticipantts:BeforBoost Pproj="<<Pprojectile
3335 <<",Ptarg="<<Ptarget<<G4endl;
3336#endif
3337 G4LorentzRotation toLab(toCms.inverse()); // Boost Rotation to LabSys (LS)
3338 Pprojectile.transform(toCms);
3339 Ptarget.transform(toCms);
3340#ifdef debug
3341 G4cout<< "G4QInelast::ExciteDiffParticipantts: AfterBoost Pproj="<<Pprojectile<<",Ptarg="
3342 <<Ptarget<<", cms4M="<<Pprojectile+Ptarget<<G4endl;
3343 G4cout<<"G4QIonIonCollis::ExciteDiffParticipants:ProjX+="<<Pprojectile.plus()<<",ProjX-="
3344 <<Pprojectile.minus()<<", tX+="<< Ptarget.plus()<<",tX-="<<Ptarget.minus()<<G4endl;
3345#endif
3346 G4LorentzVector Qmomentum(0.,0.,0.,0.);
3347 G4int whilecount=0;
3348#ifdef debug
3349 G4cout<<"G4QIonIonCollision::ExciteDiffParticipants: Before DO"<<G4endl;
3350#endif
3351 do
3352 {
3353 // Generate pt
3354 G4double maxPtSquare=sqr(Ptarget.pz());
3355#ifdef debug
3356 G4cout<<"G4QIonIonCollision::ExciteDiffParticipants: maxPtSq="<<maxPtSquare<<G4endl;
3357 if(whilecount++>=500 && whilecount%100==0) // @@ M.K. Hardwired limits
3358 G4cout<<"G4QIonIonCollision::ExciteDiffParticipantts: can loop, loopCount="<<whilecount
3359 <<", maxPtSquare="<<maxPtSquare<<G4endl;
3360#endif
3361 if(whilecount>1000) // @@ M.K. Hardwired limits
3362 {
3363#ifdef debug
3364 G4cout<<"G4QIonIonCollision::ExciteDiffParticipants: *2* abort Loop!! *2*"<<G4endl;
3365#endif
3366 return false; // Ignore this interaction
3367 }
3368 Qmomentum=G4LorentzVector(GaussianPt(widthOfPtSquare,maxPtSquare),0);
3369#ifdef debug
3370 G4cout<<"G4QIonIonCollision::ExciteDiffParticipants: generated Pt="<<Qmomentum
3371 <<", ProjPt="<<Pprojectile+Qmomentum<<", TargPt="<<Ptarget-Qmomentum<<G4endl;
3372#endif
3373 // Momentum transfer
3374 G4double Xmin=0.;
3375 G4double Xmax=1.;
3376 G4double Xplus =ChooseX(Xmin,Xmax);
3377 G4double Xminus=ChooseX(Xmin,Xmax);
3378#ifdef debug
3379 G4cout<<"G4QIonIonCollision::ExciteDiffParticipant:X+="<<Xplus<<",X-="<<Xminus<<G4endl;
3380#endif
3381 G4double pt2=Qmomentum.vect().mag2();
3382 G4double Qplus =-pt2/Xminus/Ptarget.minus();
3383 G4double Qminus= pt2/Xplus /Pprojectile.plus();
3384 Qmomentum.setPz((Qplus-Qminus)/2);
3385 Qmomentum.setE( (Qplus+Qminus)/2);
3386#ifdef debug
3387 G4cout<<"G4QInelast::ExciteDiffParticip: Qplus="<<Qplus<<", Qminus="<<Qminus<<", pt2="
3388 <<pt2<<", Qmomentum="<<Qmomentum<<", ProjM="<<(Pprojectile+Qmomentum).mag()
3389 <<", TargM="<<(Ptarget-Qmomentum).mag()<<G4endl;
3390#endif
3391 } while((Pprojectile+Qmomentum).mag2()<=Mprojectile2 ||
3392 (Ptarget-Qmomentum).mag2()<=Mtarget2);
3393 Pprojectile += Qmomentum;
3394 Ptarget -= Qmomentum;
3395#ifdef debug
3396 G4cout<<"G4QInelast::ExciteDiffParticipan: Proj(Q)="<<Pprojectile<<", Targ(Q)="<<Ptarget
3397 <<", Proj(back)="<<toLab*Pprojectile<<", Targ(bac)="<< toLab*Ptarget << G4endl;
3398#endif
3399 // Transform back and update SplitableHadron Participant.
3400 Pprojectile.transform(toLab);
3401 Ptarget.transform(toLab);
3402#ifdef debug
3403 G4cout<< "G4QIonIonCollision::ExciteDiffParticipants:TargetMass="<<Ptarget.mag()<<G4endl;
3404#endif
3405 target->Set4Momentum(Ptarget);
3406#ifdef debug
3407 G4cout<<"G4QIonIonCollision::ExciteDiffParticipant:ProjMass="<<Pprojectile.mag()<<G4endl;
3408#endif
3409 projectile->Set4Momentum(Pprojectile);
3410 return true;
3411} // End of ExciteDiffParticipants
G4LorentzVector
CLHEP::HepLorentzVector G4LorentzVector
Definition: G4LorentzVector.hh:35
CLHEP::HepLorentzVector::boostVector
Hep3Vector boostVector() const
Definition: LorentzVector.cc:177
G4QIonIonCollision::GaussianPt
G4ThreeVector GaussianPt(G4double widthSquare, G4double maxPtSquare) const
Definition: G4QIonIonCollision.cc:3572
G4QIonIonCollision::ChooseX
G4double ChooseX(G4double Xmin, G4double Xmax) const
Definition: G4QIonIonCollision.cc:3551

Referenced by G4QIonIonCollision().

◆ ExciteSingDiffParticipants()

G4bool G4QIonIonCollision::ExciteSingDiffParticipants ( G4QHadron *  aPartner,
G4QHadron *  bPartner 
) const
protected

Definition at line 3415 of file G4QIonIonCollision.cc.

3417{
3418 G4LorentzVector Pprojectile=projectile->Get4Momentum();
3419 G4double Mprojectile=projectile->GetMass();
3420 G4double Mprojectile2=Mprojectile*Mprojectile;
3421 G4LorentzVector Ptarget=target->Get4Momentum();
3422 G4double Mtarget=target->GetMass();
3423 G4double Mtarget2=Mtarget*Mtarget;
3424#ifdef debug
3425 G4cout<<"G4QInel::ExSingDiffPartici: Ep="<<Pprojectile.e()<<", Et="<<Ptarget.e()<<G4endl;
3426#endif
3427 G4bool KeepProjectile= G4UniformRand() > 0.5;
3428 // Reset minMass of the non diffractive particle to its value, (minus for rounding...)
3429 if(KeepProjectile )
3430 {
3431#ifdef debug
3432 G4cout<<"-1/2-G4QIonIonCollision::ExSingDiffParticipants: Projectile is fixed"<<G4endl;
3433#endif
3434 Mprojectile2 = projectile->GetMass2()*(1.-perCent); // Isn't it too big reduction? M.K.
3435 }
3436 else
3437 {
3438#ifdef debug
3439 G4cout<<"---1/2---G4QIonIonCollision::ExSingDiffParticipants: Target is fixed"<<G4endl;
3440#endif
3441 Mtarget2 = target->GetMass2()*(1.-perCent); // @@ Isn't it too big reduction? M.K.
3442 }
3443 // @@ From this point it repeats the Diffractional excitation (? Use flag ?)
3444 // Transform momenta to cms and then rotate parallel to z axis;
3445 G4LorentzVector Psum=Pprojectile+Ptarget;
3446 G4LorentzRotation toCms(-Psum.boostVector()); // Boost Rotation to CMS
3447 G4LorentzVector Ptmp=toCms*Pprojectile;
3448 if(Ptmp.pz()<=0.) // "String" moving backwards in CMS, abort collision !! ? M.K.
3449 {
3450#ifdef debug
3451 G4cout<<"G4QIonIonCollision::ExciteSingDiffParticipants: *1*abortCollision*1*"<<G4endl;
3452#endif
3453 return false;
3454 }
3455 toCms.rotateZ(-Ptmp.phi());
3456 toCms.rotateY(-Ptmp.theta());
3457#ifdef debug
3458 G4cout<<"G4QInel::ExciteSingDiffParticipantts: Be4Boost Pproj="<<Pprojectile<<", Ptarg="
3459 <<Ptarget<<G4endl;
3460#endif
3461 G4LorentzRotation toLab(toCms.inverse()); // Boost Rotation to LabSys (LS)
3462 Pprojectile.transform(toCms);
3463 Ptarget.transform(toCms);
3464#ifdef debug
3465 G4cout<< "G4QInelast::ExciteDiffParticipantts: AfterBoost Pproj="<<Pprojectile<<"Ptarg="
3466 <<Ptarget<<", cms4M="<<Pprojectile+Ptarget<<G4endl;
3467
3468 G4cout<<"G4QInelast::ExciteDiffParticipantts: ProjX+="<<Pprojectile.plus()<<", ProjX-="
3469 <<Pprojectile.minus()<<", TargX+="<< Ptarget.plus()<<", TargX-="<<Ptarget.minus()
3470 <<G4endl;
3471#endif
3472 G4LorentzVector Qmomentum(0.,0.,0.,0.);
3473 G4int whilecount=0;
3474 do
3475 {
3476 // Generate pt
3477 G4double maxPtSquare=sqr(Ptarget.pz());
3478 if(whilecount++>=500 && whilecount%100==0) // @@ M.K. Hardwired limits
3479#ifdef debug
3480 G4cout<<"G4QIonIonCollision::ExciteSingDiffParticipantts: can loop, loopCount="
3481 <<whilecount<<", maxPtSquare="<<maxPtSquare<<G4endl;
3482#endif
3483 if(whilecount>1000) // @@ M.K. Hardwired limits
3484 {
3485#ifdef debug
3486 G4cout<<"G4QIonIonCollision::ExciteSingDiffParticipants:*2* abortLoop!! *2*"<<G4endl;
3487#endif
3488 return false; // Ignore this interaction
3489 }
3490 Qmomentum=G4LorentzVector(GaussianPt(widthOfPtSquare,maxPtSquare),0);
3491#ifdef debug
3492 G4cout<<"G4QInelast::ExciteSingDiffParticipants: generated Pt="<<Qmomentum
3493 <<", ProjPt="<<Pprojectile+Qmomentum<<", TargPt="<<Ptarget-Qmomentum<<G4endl;
3494#endif
3495 // Momentum transfer
3496 G4double Xmin=0.;
3497 G4double Xmax=1.;
3498 G4double Xplus =ChooseX(Xmin,Xmax);
3499 G4double Xminus=ChooseX(Xmin,Xmax);
3500#ifdef debug
3501 G4cout<<"G4QInel::ExciteSingDiffPartici: X-plus="<<Xplus<<", X-minus="<<Xminus<<G4endl;
3502#endif
3503 G4double pt2=G4ThreeVector(Qmomentum.vect()).mag2();
3504 G4double Qplus =-pt2/Xminus/Ptarget.minus();
3505 G4double Qminus= pt2/Xplus /Pprojectile.plus();
3506 if (KeepProjectile)
3507 Qminus=(projectile->GetMass2()+pt2)/(Pprojectile.plus()+Qplus) - Pprojectile.minus();
3508 else Qplus=Ptarget.plus() - (target->GetMass2()+pt2)/(Ptarget.minus()-Qminus);
3509 Qmomentum.setPz((Qplus-Qminus)/2);
3510 Qmomentum.setE( (Qplus+Qminus)/2);
3511#ifdef debug
3512 G4cout<<"G4QInel::ExciteDiffParticip: Qplus="<<Qplus<<", Qminus="<<Qminus<<", pt2="
3513 <<pt2<<", Qmomentum="<<Qmomentum<<", ProjM="<<(Pprojectile+Qmomentum).mag()
3514 <<", TargM="<<(Ptarget-Qmomentum).mag()<<G4endl;
3515#endif
3516 // while is different from the Double Diffractive Excitation (@@ !)
3517 //} while((Pprojectile+Qmomentum).mag2()<= Mprojectile2 ||
3518 // (Ptarget-Qmomentum).mag2()<=Mtarget2);
3519 } while((Ptarget-Qmomentum).mag2()<=Mtarget2 ||
3520 (Pprojectile+Qmomentum).mag2()<=Mprojectile2 ||
3521 (Ptarget-Qmomentum).e() < 0. || (Pprojectile+Qmomentum).e() < 0.);
3522 Pprojectile += Qmomentum;
3523 Ptarget -= Qmomentum;
3524#ifdef debug
3525 G4cout<<"G4QIonIonCollision::ExciteSingDiffParticipan: Proj(Q)="<<Pprojectile<<"(E="
3526 <<Pprojectile.e()<<"), Targ(Q)="<<Ptarget<<"(E="<<Ptarget.e()
3527 <<"), Proj(back)="<<toLab*Pprojectile<<", Targ(bac)="<< toLab*Ptarget << G4endl;
3528#endif
3529 // Transform back and update SplitableHadron Participant.
3530 Pprojectile.transform(toLab);
3531 Ptarget.transform(toLab);
3532#ifdef debug
3533 G4cout<< "G4QIonIonCollision::ExciteSingleDiffParticipants: TgM="<<Ptarget.mag()<<G4endl;
3534#endif
3535 target->Set4Momentum(Ptarget);
3536#ifdef debug
3537 G4cout<<"G4QIonIonCollision::ExciteSingleParticipants:ProjM="<<Pprojectile.mag()<<G4endl;
3538#endif
3539 projectile->Set4Momentum(Pprojectile);
3540 return true;
3541} // End of ExciteSingleDiffParticipants
G4ThreeVector
CLHEP::Hep3Vector G4ThreeVector
Definition: G4ThreeVector.hh:42
CLHEP::Hep3Vector::mag2
double mag2() const

Referenced by G4QIonIonCollision().

◆ Fragment()

G4QHadronVector * G4QIonIonCollision::Fragment ( )

Definition at line 1557 of file G4QIonIonCollision.cc.

1558{ // This is the member function fragmenting Strings & Quasmons (in nuclear matter)
1559#ifdef debug
1560 G4cout<<"*******>G4QIonIonCollision::Fragment: ***Called***, Res="<<theResult<<G4endl;
1561#endif
1562 G4int striNum=strings.size(); // Find out if there're strings
1563 G4int hadrNum=theResult->size(); // Find out if there're hadrons
1564#ifdef edebug
1565 G4int nQm=theQuasmons.size();
1566 G4LorentzVector totLS4M=theProjNucleus.Get4Momentum()+theTargNucleus.Get4Momentum(); //LS
1567 G4int totChg=theProjNucleus.GetZ()+theTargNucleus.GetZ();
1568 G4int totBaN=theTargNucleus.GetA()+theTargNucleus.GetA();
1569 G4cout<<"-EMCLS-G4QIn::Fragment: CHECKRecovery, #ofS="<<striNum<<", Nuc4M(E=M)="<<totLS4M
1570 <<",#Q="<<nQm<<",#H="<<hadrNum<<G4endl;
1571 for(G4int i=0; i < striNum; i++)
1572 {
1573 G4LorentzVector strI4M=strings[i]->Get4Momentum();
1574 totLS4M+=strI4M;
1575 G4int sChg=strings[i]->GetCharge();
1576 totChg+=sChg;
1577 G4int sBaN=strings[i]->GetBaryonNumber();
1578 totBaN+=sBaN;
1579 G4cout<<"-EMCLS-G4QIonIonCollision::Fragment:String#"<<i<<", 4M="<<strI4M<<", M="
1580 <<strI4M.m()<<", C="<<sChg<<", B="<<sBaN<<G4endl;
1581 }
1582 for(G4int i=0; i < nQm; i++)
1583 {
1584 G4LorentzVector hI4M=theQuasmons[i]->Get4Momentum();
1585 totLS4M+=hI4M;
1586 G4int hChg=theQuasmons[i]->GetCharge();
1587 totChg+=hChg;
1588 G4int hBaN=theQuasmons[i]->GetBaryonNumber();
1589 totBaN+=hBaN;
1590 G4cout<<"-EMCLS-G4QIonIonCollision::Fragment: Quasmon#"<<i<<", 4M="<<hI4M<<", C="<<hChg
1591 <<", B="<<hBaN<<G4endl;
1592 }
1593 for(G4int i=0; i < hadrNum; i++)
1594 {
1595 G4LorentzVector hI4M=(*theResult)[i]->Get4Momentum();
1596 totLS4M+=hI4M;
1597 G4int hChg=(*theResult)[i]->GetCharge();
1598 totChg+=hChg;
1599 G4int hBaN=(*theResult)[i]->GetBaryonNumber();
1600 totBaN+=hBaN;
1601 G4cout<<"-EMCLS-G4QIn::Fragment:H#"<<i<<",4M="<<hI4M<<",C="<<hChg<<",B="<<hBaN<<G4endl;
1602 }
1603#endif
1604#ifdef debug
1605 G4cout<<"***>G4QIonIonCollision::Fragm: #OfStr="<<striNum<<", #OfRes="<<hadrNum<<G4endl;
1606#endif
1607 if(!striNum && hadrNum) // Quasi-elastic or decoupled p
1608 {
1609#ifdef debug
1610 G4cout<<"***>G4QIonIonCollision::Fragm:**Quasi-Elastic**, #OfResult="<<hadrNum<<G4endl;
1611#endif
1612 return theResult;
1613 }
1614 else if(striNum) Breeder(); // Strings fragmentation
1615 else // No strings, make HadrNucleus
1616 {
1617 if(hadrNum)
1618 {
1619 for(G4int ih=0; ih<hadrNum; ih++) delete (*theResult)[ih];
1620 theResult->clear();
1621 }
1622 G4LorentzVector rp4M=theProjNucleus.Get4Momentum(); // Nucleus 4-momentum in LS
1623 G4int rpPDG=theProjNucleus.GetPDG(); // Nuclear PDG
1624 G4QHadron* resPNuc = new G4QHadron(rpPDG,rp4M); // Nucleus -> Hadron
1625 theResult->push_back(resPNuc); // Fill the residual nucleus
1626 G4LorentzVector rt4M=theTargNucleus.Get4Momentum(); // Nucleus 4-momentum in LS
1627 G4int rtPDG=theTargNucleus.GetPDG(); // Nuclear PDG
1628 G4QHadron* resTNuc = new G4QHadron(rtPDG,rt4M); // Nucleus -> Hadron
1629 theResult->push_back(resTNuc); // Fill the residual nucleus
1630 }
1631 G4int nQuas=theQuasmons.size(); // Size of the Quasmon OUTPUT
1632 G4int theRS=theResult->size(); // Size of Hadron Output by now
1633#ifdef debug
1634 G4cout<<"***>G4QIonIonCollision::Fragm:beforeEnv, #OfQ="<<nQuas<<",#OfR="<<theRS<<G4endl;
1635#endif
1636 if(nQuas && theRS)
1637 {
1638
1639 G4QHadron* resNuc = (*theResult)[theRS-1]; // Pointer to Residual Nucleus
1640 G4LorentzVector resNuc4M = resNuc->Get4Momentum(); // 4-Momentum of the Nucleuz
1641 G4int resNucPDG= resNuc->GetPDGCode(); // PDG Code of the Nucleus
1642 G4QNucleus theEnv(resNucPDG); // NucleusHadron->NucleusAtRest
1643 delete resNuc; // Delete resNucleus as aHadron
1644 theResult->pop_back(); // Exclude the nucleus from HV
1645 --theRS; // Reduce the OUTPUT by theNucl
1646#ifdef pdebug
1647 G4cout<<"G4QIonIonCollision::Fragm:#OfRemainingHadron="<<theRS<<", A="<<theEnv<<G4endl;
1648#endif
1649 // Now we need to be sure that the compound nucleus is heavier than the Ground State
1650 for(G4int j=theRS-1; j>-2; --j) // try to reach M_compound>M_GS
1651 {
1652 G4LorentzVector qsum4M=resNuc4M; // Proto compound 4-momentum
1653 G4QContent qsumQC=theEnv.GetQCZNS(); // Proto compound Quark Content
1654#ifdef pdebug
1655 G4cout<<"G4QIonIonCollis::Fragm: rN4M"<<qsum4M<<qsum4M.m()<<",rNQC="<<qsumQC<<G4endl;
1656#endif
1657 G4Quasmon* firstQ=0; // Prototype of theFirstQuasmon
1658 G4LorentzVector first4M; // Proto of the FirstQuasmon 4M
1659 G4QContent firstQC; // Proto of the FirstQuasmon QC
1660 for(G4int i=0; i<nQuas; ++i) // LOOP over Quasmons
1661 {
1662 G4Quasmon* curQuasm=theQuasmons[i]; // current Quasmon
1663 G4LorentzVector cur4M=curQuasm->Get4Momentum(); // 4-Mom of the Quasmon
1664 G4QContent curQC=curQuasm->GetQC(); // Quark Content of the Quasmon
1665 qsum4M+=cur4M; // Add quasmon's 4-momentum
1666 qsumQC+=curQC; // Add quasmon's Quark Content
1667#ifdef pdebug
1668 G4cout<<"G4QIonIonCollis::Fragm: Q#"<<i<<", QQC="<<curQC<<", sQC="<<qsumQC<<G4endl;
1669#endif
1670 if(!i) // Remember 1-st for correction
1671 {
1672 firstQ =curQuasm;
1673 first4M=cur4M;
1674 firstQC=curQC;
1675 }
1676 }
1677 G4int miPDG=qsumQC.GetSPDGCode(); // PDG of minM of hadron/fragm.
1678 G4double gsM=0.; // Proto minM of had/frag forQC
1679 if(miPDG == 10)
1680 {
1681 G4QChipolino QCh(qsumQC); // define TotNuc as a Chipolino
1682 gsM=QCh.GetQPDG1().GetMass()+QCh.GetQPDG2().GetMass(); // Sum of Hadron Masses
1683 //gsM=theWorld->GetQParticle(QCh.GetQPDG1())->MinMassOfFragm() +
1684 // theWorld->GetQParticle(QCh.GetQPDG2())->MinMassOfFragm();
1685 }
1686 // @@ it is not clear, why it does not work ?!
1687 //else if(miPDG>80000000) // Compound Nucleus
1688 //{
1689 // G4QNucleus rtN(qsumQC); // CreatePseudoNucl for totComp
1690 // gsM=rtN.GetGSMass(); // MinMass of residQ+(Env-ParC)
1691 //}
1692 else if(miPDG < 80000000 && std::abs(miPDG)%10 > 2)
1693 gsM=theWorld->GetQParticle(G4QPDGCode(miPDG))->MinMassOfFragm();
1694 else gsM=G4QPDGCode(miPDG).GetMass(); // minM of hadron/fragm. for QC
1695 G4double reM=qsum4M.m(); // real mass of the compound
1696#ifdef pdebug
1697 G4cout<<"G4QIonIonCollision::Fragm: PDG="<<miPDG<<", rM="<<reM<<",GSM="<<gsM<<G4endl;
1698#endif
1699 if(reM > gsM) break; // CHIPS can be called
1700 if(j > -1) // Can try to add hadrons to Q0
1701 {
1702 G4QHadron* cH = (*theResult)[j]; // Pointer to the last Hadron
1703 G4LorentzVector h4M = cH->Get4Momentum(); // 4-Momentum of the Hadron
1704 G4QContent hQC = cH->GetQC(); // QC of the Hadron
1705 firstQ->Set4Momentum(first4M+h4M); // Update the Quasmon's 4-Mom
1706 firstQ->SetQC(firstQC+hQC); // Update the Quasmon's QCont
1707 delete cH; // Delete the Hadron
1708 theResult->pop_back(); // Exclude the hadron from HV
1709#ifdef pdebug
1710 G4cout<<"G4QInel::Fragm: H#"<<j<<", hQC="<<hQC<<",hPDG="<<cH->GetPDGCode()<<G4endl;
1711#endif
1712 }
1713 else
1714 {
1715 G4cerr<<"***G4QIonIonCollision::Fra:PDG="<<miPDG<<",M="<<reM<<",GSM="<<gsM<<G4endl;
1716 G4Exception("G4QIonIonCollision::Fragment:","27",FatalException,"Can'tRecoverGSM");
1717 }
1718 }
1719 G4double nucE=resNuc4M.e(); // Total energy of the nuclEnv
1720 if(nucE<1.E-12) nucE=0.; // Computer accuracy safety
1721 G4ThreeVector nucVel(0.,0.,0.); // Proto of the NucleusVelocity
1722 G4QHadronVector* output=0; // NucleusFragmentation Hadrons
1723 G4QEnvironment* pan= new G4QEnvironment(theEnv); // ---> DELETED --->----------+
1724#ifdef pdebug
1725 G4cout<<"G4QIonIonCollision::Fragm: nucE="<<nucE<<", nQ="<<nQuas<<G4endl; // |
1726#endif
1727 if(nucE) nucVel=resNuc4M.vect()/nucE; // The NucleusVelocity |
1728 for(G4int i=0; i<nQuas; ++i) // LOOP over Quasmons |
1729 { // |
1730 G4Quasmon* curQuasm=theQuasmons[i]; // current Quasmon |
1731 if(nucE) curQuasm->Boost(-nucVel); // Boost it to CMS of Nucleus |
1732 pan->AddQuasmon(curQuasm); // Fill the predefined Quasmon|
1733#ifdef pdebug
1734 G4LorentzVector cQ4M=curQuasm->Get4Momentum(); // Just for printing |
1735 G4cout<<"G4QIonIonCollision::Fragment: Quasmon# "<<i<<" added, 4M="<<cQ4M<<G4endl;//|
1736#endif
1737 } // |
1738 try // |
1739 { // |
1740 delete output; // |
1741 output = pan->Fragment();// DESTROYED after theHadrons are transferred to theResult |
1742 } // | |
1743 catch (G4QException& error) // | |
1744 { // | |
1745 G4cerr<<"***G4QIonIonCollision::Fragment: G4QE Exception is catched"<<G4endl; // | |
1746 // G4Exception("G4QIonIonCollision::Fragment:","27",FatalException,"CHIPSCrash");// | |
1747 G4Exception("G4QIonIonCollision::Fragment()", "HAD_CHPS_0027",
1748 FatalException, "CHIPSCrash");
1749 } // | |
1750 delete pan; // Delete the Nuclear Environment <-----<--+-+
1751 if(output) // Output exists |
1752 { // |
1753 G4int nOut=output->size(); // #ofHadrons in the Nuclear Fragmentation |
1754 for(G4int j=0; j<nOut; j++) // LOOP over Hadrons transferring to LS |
1755 { // |
1756 G4QHadron* curHadron=(*output)[j]; // Hadron from the nucleus fragmentation |
1757 if(nucE) curHadron->Boost(nucVel); // Boost it back to Laboratory System |
1758 theResult->push_back(curHadron); // Transfer it to the result |
1759 } // |
1760 delete output; // Delete the OUTPUT <-----<-----<-----<---+
1761 }
1762 }
1763 else if(!striNum) G4cout<<"-Warning-G4QIonIonCollision::Fragment:NothingWasDone"<<G4endl;
1764#ifdef debug
1765 G4cout<<"=--=>G4QIonIonCollision::Fragment: Final #OfResult="<<theResult->size()<<G4endl;
1766#endif
1767 G4int nQ =theQuasmons.size();
1768 if(nQ) theQuasmons.clear(); // @@ Not necesary ?
1769#ifdef edebug
1770 G4LorentzVector f4M(0.,0.,0.,0.); // Sum of the Result in LS
1771 G4int fCh=totChg;
1772 G4int fBN=totBaN;
1773 G4int nHd=theResult->size();
1774 G4cout<<"-EMCLS-G4QIonIonCollision::Fragment: #ofHadr="<<nHd<<",#OfQuasm="<<nQ<<G4endl;
1775 for(G4int i=0; i<nHd; i++)
1776 {
1777 G4LorentzVector hI4M=(*theResult)[i]->Get4Momentum();
1778 f4M+=hI4M;
1779 G4int hChg=(*theResult)[i]->GetCharge();
1780 fCh-=hChg;
1781 G4int hBaN=(*theResult)[i]->GetBaryonNumber();
1782 fBN-=hBaN;
1783 G4cout<<"-EMCLS-G4QIonIonCollision::Fragment: Hadron#"<<i<<", 4M="<<hI4M<<", PDG="
1784 <<(*theResult)[i]->GetPDGCode()<<", C="<<hChg<<", B="<<hBaN<<G4endl;
1785 }
1786 G4cout<<"-EMCLS-G4QInel::Fragm:, r4M="<<f4M-totLS4M<<", rC="<<fCh<<",rB="<<fBN<<G4endl;
1787#endif
1788 return theResult;
1789} // End of fragmentation
G4QCHIPSWorld::GetQParticle
G4QParticle * GetQParticle(G4int PDG) const
Definition: G4QCHIPSWorld.hh:87
G4QEnvironment::AddQuasmon
void AddQuasmon(G4Quasmon *Q)
Definition: G4QEnvironment.cc:10522
G4QEnvironment::Fragment
G4QHadronVector * Fragment()
Definition: G4QEnvironment.cc:4512
G4QHadron::Boost
void Boost(const G4LorentzVector &theBoost)
Definition: G4QHadron.cc:1293
G4QParticle::MinMassOfFragm
G4double MinMassOfFragm()
Definition: G4QParticle.hh:113
G4Quasmon::GetQC
G4QContent GetQC() const
Definition: G4Quasmon.hh:162
G4Quasmon::Set4Momentum
void Set4Momentum(G4LorentzVector Q4M)
Definition: G4Quasmon.hh:96
G4Quasmon::Boost
void Boost(const G4LorentzVector &theBoost)
Definition: G4Quasmon.cc:6202
G4Quasmon::Get4Momentum
G4LorentzVector Get4Momentum() const
Definition: G4Quasmon.hh:161
G4Quasmon::SetQC
void SetQC(G4QContent QQC)
Definition: G4Quasmon.hh:97

Referenced by G4QInelastic::PostStepDoIt().

◆ GaussianPt()

G4ThreeVector G4QIonIonCollision::GaussianPt ( G4double  widthSquare,
G4double  maxPtSquare 
) const
protected

Definition at line 3572 of file G4QIonIonCollision.cc.

3573{
3574#ifdef debug
3575 G4cout<<"G4QIonIonCollision::GaussianPt: w^2="<<widthSq<<",maxPt2="<<maxPtSquare<<G4endl;
3576#endif
3577 G4double pt2=0.;
3578 G4double rm=maxPtSquare/widthSq; // Negative
3579 if(rm>-.01) pt2=widthSq*(std::sqrt(1.-G4UniformRand()*(1.-sqr(1.+rm)))-1.);
3580 else pt2=widthSq*std::log(1.-G4UniformRand()*(1.-std::exp(rm)));
3581 pt2=std::sqrt(pt2); // It is not pt2, but pt now
3582 G4double phi=G4UniformRand()*twopi;
3583 return G4ThreeVector(pt2*std::cos(phi),pt2*std::sin(phi),0.);
3584} // End of GaussianPt

Referenced by ExciteDiffParticipants(), and ExciteSingDiffParticipants().

◆ IsSingleDiffractive()

G4bool G4QIonIonCollision::IsSingleDiffractive ( )
inlineprotected

Definition at line 76 of file G4QIonIonCollision.hh.

76{G4bool r=false; if(G4UniformRand()<1.) r=true; return r;}

Referenced by G4QIonIonCollision().

◆ ReducePair()

std::pair< G4int, G4int > G4QIonIonCollision::ReducePair ( G4int  P1,
G4int  P2 
) const
protected

Definition at line 3651 of file G4QIonIonCollision.cc.

3652{
3653#ifdef debug
3654 G4cout<<"G4QIonIonCollision::ReducePair: **Called** P1="<<P1<<", P2="<<P2<<G4endl;
3655#endif
3656 G4int P11=P1/10;
3657 G4int P12=P1%10;
3658 G4int P21=P2/10;
3659 G4int P22=P2%10;
3660 if (P11==P21) return std::make_pair(P12,P22);
3661 else if(P11==P22) return std::make_pair(P12,P21);
3662 else if(P12==P21) return std::make_pair(P11,P22);
3663 else if(P12==P22) return std::make_pair(P11,P21);
3664 //#ifdef debug
3665 G4cout<<"-Warning-G4QIonIonCollision::ReducePair:**Failed**,P1="<<P1<<",P2="<<P2<<G4endl;
3666 //#endif
3667 return std::make_pair(0,0); // Reduction failed
3668}

Referenced by Breeder().

◆ SetParameters()

void G4QIonIonCollision::SetParameters ( G4int  nC,
G4double  strTens,
G4double  tubeDens,
G4double  SigPt 
)
static

Definition at line 3543 of file G4QIonIonCollision.cc.

3544{
3545 nCutMax = nC; // max number of pomeron cuts
3546 stringTension = stT; // string tension for absorbed energy
3547 tubeDensity = tbD; // Flux Tube Density of nuclear nucleons
3548 widthOfPtSquare = -2*SigPt*SigPt; // width^2 of pt for string excitation
3549}

◆ SumPartonPDG()

G4int G4QIonIonCollision::SumPartonPDG ( G4int  PDG1,
G4int  PFG2 
) const
protected

Definition at line 3586 of file G4QIonIonCollision.cc.

3587{
3588 if (PDG1 < 7 && PDG1 > 0 && PDG2 < 7 && PDG2 > 0) // Sum up two Q in DiQ (S=0)
3589 {
3590 if(PDG1 > PDG2) return PDG1*1000+PDG2*100+1;
3591 else return PDG2*1000+PDG1*100+1;
3592 }
3593 else if (PDG1 >-7 && PDG1 < 0 && PDG2 >-7 && PDG2 < 0) // Sum up two AQ in ADiQ (S=0)
3594 {
3595 if(-PDG1 > -PDG2) return PDG1*1000+PDG2*100-1;
3596 else return PDG2*1000+PDG1*100-1;
3597 }
3598 else if (PDG1 <-99 && PDG2 < 7 && PDG2 > 0) // Sum up Q and ADiQ in AQ
3599 {
3600 G4int PDG=-PDG1/100;
3601 if(PDG2==PDG/10) return -PDG%10;
3602 if(PDG2==PDG%10) return -PDG/10;
3603 else
3604 {
3605 G4cerr<<"***4QFragmentation::SumPartonPDG: PDG1="<<PDG1<<", PDG2="<<PDG2<<G4endl;
3606 G4Exception("G4QIonIonCollision::SumPartonPDG:","72",FatalException,"Q&ADiQnoMatch");
3607 }
3608 }
3609 else if (PDG2 <-99 && PDG1 < 7 && PDG1 > 0) // Sum up ADiQ and Q in AQ
3610 {
3611 G4int PDG=-PDG2/100;
3612 if(PDG1==PDG/10) return -PDG%10;
3613 if(PDG1==PDG%10) return -PDG/10;
3614 else
3615 {
3616 G4cerr<<"***4QFragmentation::SumPartonPDG: PDG1="<<PDG1<<", PDG2="<<PDG2<<G4endl;
3617 G4Exception("G4QIonIonCollision::SumPartonPDG:","72",FatalException,"ADiQ&QnoMatch");
3618 }
3619 }
3620 else if (PDG1 > 99 && PDG2 >-7 && PDG2 < 0) // Sum up DiQ and AQ in Q
3621 {
3622 G4int PDG=PDG1/100;
3623 if(PDG2==-PDG/10) return PDG%10;
3624 if(PDG2==-PDG%10) return PDG/10;
3625 else
3626 {
3627 G4cerr<<"***4QFragmentation::SumPartonPDG: PDG1="<<PDG1<<", PDG2="<<PDG2<<G4endl;
3628 G4Exception("G4QIonIonCollision::SumPartonPDG:","72",FatalException,"DiQ&AQnoMatch");
3629 }
3630 }
3631 else if (PDG2 > 99 && PDG1 >-7 && PDG1 < 0) // Sum up AQ and DiQ in Q
3632 {
3633 G4int PDG=PDG2/100;
3634 if(PDG1==-PDG/10) return PDG%10;
3635 if(PDG1==-PDG%10) return PDG/10;
3636 else
3637 {
3638 G4cerr<<"***G4QIonIonCollision::SumPartonPDG: PDG1="<<PDG1<<", PDG2="<<PDG2<<G4endl;
3639 G4Exception("G4QIonIonCollision::SumPartonPDG:","72",FatalException,"AQ&DiQnoMatch");
3640 }
3641 }
3642 else
3643 {
3644 G4cerr<<"***G4QIonIonCollision::SumPartonPDG: PDG1="<<PDG1<<", PDG2="<<PDG2<<G4endl;
3645 G4Exception("G4QIonIonCollision::SumPartonPDG:","72",FatalException,"Can'tSumUpParts");
3646 }
3647 return 0;
3648} // End of SumPartonPDG

Referenced by Breeder().

◆ SwapPartons()

void G4QIonIonCollision::SwapPartons ( )
protected

Definition at line 3903 of file G4QIonIonCollision.cc.

3904{
3905 static const G4double baryM=800.; // Mass excess for baryons
3906 G4QStringVector::iterator ist;
3907 for(ist = strings.begin(); ist < strings.end(); ist++)
3908 {
3909 G4LorentzVector cS4M=(*ist)->Get4Momentum();
3910 G4double cSM2=cS4M.m2(); // Squared mass of the String
3911 if(cSM2<0.1) // Parton Swapping is needed
3912 {
3913 G4QParton* cLeft=(*ist)->GetLeftParton(); // Current String Left Parton
3914 G4QParton* cRight=(*ist)->GetRightParton(); // Current String Right Parton
3915 G4int cLPDG=cLeft->GetPDGCode();
3916 G4int cRPDG=cRight->GetPDGCode();
3917 G4LorentzVector cL4M=cLeft->Get4Momentum();
3918 G4LorentzVector cR4M=cRight->Get4Momentum();
3919 G4int cLT=cLeft->GetType();
3920 G4int cRT=cRight->GetType();
3921 G4QStringVector::iterator sst; // Selected partner string
3922 G4QStringVector::iterator pst; // LOOP iterator
3923 G4double maxM=-DBL_MAX; // Swapping providing the max mass
3924 G4int sDir=0; // Selected direction of swapping
3925 for(pst = strings.begin(); pst < strings.end(); pst++) if(pst != ist)
3926 {
3927 G4QParton* pLeft=(*pst)->GetLeftParton(); // Partner String Left Parton
3928 G4QParton* pRight=(*pst)->GetRightParton(); // Partner String Right Parton
3929 G4int pLPDG=pLeft->GetPDGCode();
3930 G4int pRPDG=pRight->GetPDGCode();
3931 G4LorentzVector pL4M=pLeft->Get4Momentum();
3932 G4LorentzVector pR4M=pRight->Get4Momentum();
3933 G4int pLT=pLeft->GetType();
3934 G4int pRT=pRight->GetType();
3935 G4double LM=0.;
3936 G4double RM=0.;
3937 if(((cLPDG<-7 || (cLPDG>0 && cLPDG< 7) ) && (pLPDG<-7 || (pLPDG>0 && pLPDG< 7) ))||
3938 ((cLPDG> 7 || (cLPDG<0 && cLPDG>-7) ) && (pLPDG> 7 || (pLPDG<0 && cLPDG>-7) )))
3939 {
3940 G4double pLM2=(cL4M+pR4M).m2(); // new partner M2
3941 G4double cLM2=(cR4M+pL4M).m2(); // new partner M2
3942 if(pLM2>0. && cLM2>0.)
3943 {
3944 G4double pLM=std::sqrt(pLM2);
3945 if(cLT+pRT==3) pLM-=baryM;
3946 G4double cLM=std::sqrt(cLM2);
3947 if(cRT+pLT==3) cLM-=baryM;
3948 LM=std::min(pLM2,cLM2);
3949 }
3950 }
3951 if(((cRPDG<-7 || (cRPDG>0 && cRPDG< 7) ) && (pRPDG<-7 || (pRPDG>0 && pRPDG< 7) ))||
3952 ((cRPDG> 7 || (cRPDG<0 && cRPDG>-7) ) && (pRPDG> 7 || (pRPDG<0 && cRPDG>-7) )) )
3953 {
3954 G4double pRM2=(cR4M+pL4M).m2(); // new partner M2
3955 G4double cRM2=(cL4M+pR4M).m2(); // new partner M2
3956 if(pRM2>0. && cRM2>0.)
3957 {
3958 G4double pRM=std::sqrt(pRM2);
3959 if(cRT+pLT==3) pRM-=baryM;
3960 G4double cRM=std::sqrt(cRM2);
3961 if(cLT+pRT==3) cRM-=baryM;
3962 RM=std::min(pRM,cRM);
3963 }
3964 }
3965 G4int dir=0;
3966 G4double sM=0.;
3967 if( LM && LM > RM )
3968 {
3969 dir= 1;
3970 sM=LM;
3971 }
3972 else if(RM)
3973 {
3974 dir=-1;
3975 sM=RM;
3976 }
3977 if(sM > maxM)
3978 {
3979 sst=pst;
3980 maxM=sM;
3981 sDir=dir;
3982 }
3983 }
3984 if(sDir)
3985 {
3986 G4QParton* pLeft=(*sst)->GetLeftParton(); // Partner String Left Parton
3987 G4QParton* pRight=(*sst)->GetRightParton(); // Partner String Right Parton
3988 G4QParton* swap=pLeft; // Prototype remember the partner Left
3989 if(sDir>0) // swap left partons
3990 {
3991 (*sst)->SetLeftParton(cLeft);
3992 (*ist)->SetLeftParton(swap);
3993 }
3994 else
3995 {
3996 swap=pRight;
3997 (*sst)->SetRightParton(cRight);
3998 (*ist)->SetRightParton(swap);
3999 }
4000 }
4001#ifdef debug
4002 else G4cout<<"***G4QIonIonCollision::SwapPartons:**Failed**,cLPDG="<<cLPDG<<",cRPDG="
4003 <<cRPDG<<",-->cM2="<<cSM2<<G4endl;
4004#endif
4005
4006 }
4007 }
4008}

Referenced by G4QIonIonCollision().

The documentation for this class was generated from the following files: