Context Navigation

source: trunk/source/processes/hadronic/models/chiral_inv_phase_space/fragmentation/src/G4QIonIonCollision.cc @ 1315

Last change on this file since 1315 was 1315, checked in by garnier, 14 years ago
update geant4-09-04-beta-cand-01 interfaces-V09-03-09 vis-V09-03-08
File size: 179.6 KB

Line
1	//
2	// ********************************************************************
3	// * License and Disclaimer *
4	// * *
5	// * The Geant4 software is copyright of the Copyright Holders of *
6	// * the Geant4 Collaboration. It is provided under the terms and *
7	// * conditions of the Geant4 Software License, included in the file *
8	// * LICENSE and available at http://cern.ch/geant4/license . These *
9	// * include a list of copyright holders. *
10	// * *
11	// * Neither the authors of this software system, nor their employing *
12	// * institutes,nor the agencies providing financial support for this *
13	// * work make any representation or warranty, express or implied, *
14	// * regarding this software system or assume any liability for its *
15	// * use. Please see the license in the file LICENSE and URL above *
16	// * for the full disclaimer and the limitation of liability. *
17	// * *
18	// * This code implementation is the result of the scientific and *
19	// * technical work of the GEANT4 collaboration. *
20	// * By using, copying, modifying or distributing the software (or *
21	// * any work based on the software) you agree to acknowledge its *
22	// * use in resulting scientific publications, and indicate your *
23	// * acceptance of all terms of the Geant4 Software license. *
24	// ********************************************************************
25	//
26	//
27	// $Id: G4QIonIonCollision.cc,v 1.8 2010/04/01 15:03:35 mkossov Exp $
28	// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
29	//
30	// -----------------------------------------------------------------------------
31	// GEANT 4 class header file
32	//
33	// History:
34	// Created by Mikhail Kossov, October 2006
35	// CHIPS QGS fragmentation class
36	// For comparison similar member functions can be found in the G4 classes:
37	// G4QGSParticipants
38	// G4QGSModels
39	// G4ExcitedStringDecay
40	// -----------------------------------------------------------------------------
41	// Short description: CHIPS QG string fragmentation class
42	// Rhe key member function is Scatter, making the interaction (see G4QCollision)
43	// -----------------------------------------------------------------------------
44	//#define debug
45	//#define edebug
46	//#define pdebug
47	//#define sdebug
48	//#define ppdebug
49
50	#include "G4QIonIonCollision.hh"
51
52	// Promoting model parameters from local variables class properties @@(? M.K.)
53
54	// Definition of static parameters
55	G4int G4QIonIonCollision::nCutMax=7;
56	G4double G4QIonIonCollision::stringTension=1.*GeV/fermi;
57	G4double G4QIonIonCollision::tubeDensity =1./fermi;
58	// Parameters of diffractional fragmentation
59	G4double G4QIonIonCollision::widthOfPtSquare=-0.75GeVGeV; // ptWidth2 forStringExcitation
60
61	G4QIonIonCollision::G4QIonIonCollision(G4QNucleus &pNucleus, const G4QNucleus &tNucleus)
62	{
63	static const G4double mProt= G4QPDGCode(2212).GetMass(); // Mass of proton
64	static const G4double mPi0= G4QPDGCode(111).GetMass(); // Mass of Pi0
65	theWorld= G4QCHIPSWorld::Get(); // Get a pointer to the CHIPS World
66	theResult = new G4QHadronVector; // Define theResultingHadronVector
67	G4bool stringsInitted=false; // Strings are initiated
68	G4LorentzRotation toZ; // Lorentz Transformation to the projectileSystem
69	G4LorentzVector proj4M=pNucleus.Get4Momentum(); // Projectile nucleus 4-momentum in LS
70	//G4double projM=pNucleus.GetGSMass(); // Ground state mass of the projectile nucleus
71	G4double targM=tNucleus.GetGSMass(); // Ground state mass of the target nucleus
72	#ifdef edebug
73	G4cout<<"G4QIn::Constr:Called,pA="<<pNucleus<<proj4M<<",tA="<<tNucleus<<targM<<G4endl;
74	#endif
75	G4int pZ=pNucleus.GetZ();
76	G4int pN=pNucleus.GetN();
77	G4int pA=pZ+pN;
78	G4int pPDG=90000000+pZ*1000+pN;
79	G4int tZ=tNucleus.GetZ();
80	G4int tN=tNucleus.GetN();
81	G4int tA=tZ+tN;
82	G4int tPDG=90000000+tZ*1000+tN;
83	toZ.rotateZ(-proj4M.phi());
84	toZ.rotateY(-proj4M.theta());
85	G4LorentzVector zProj4M=toZ*proj4M; // Proj 4-momentum in LS rotated to Z axis
86	pNucleus.Set4Momentum(zProj4M); // Now the projectile nucleus moves along Z axis
87	#ifdef edebug
88	G4int totChg=pZ+tZ; // Charge of the projectile+target for the CHECK
89	G4int totBaN=pA+tA; // Baryon Number of Proj+Targ for CHECK
90	G4LorentzVector tgLS4M(0.,0.,0.,targM); // Target 4-momentum in LS
91	G4LorentzVector totLS4M=proj4M+tgLS4M; // Total 4-momentum in LS
92	G4LorentzVector totZLS4M=zProj4M+tgLS4M;// Total 4-momentum in LS with momentum along Z
93	G4cout<<"-EMC-G4QIonIonCollision::Constr: tLS4M="<<totLS4M<<",tZLS4M="<<totZLS4M<<G4endl;
94	// === From nere all consideration is made in the rotated LS frame (proj is along Z) ===
95	#endif
96	G4LorentzRotation toLab(toZ.inverse()); // Lorentz Transfornation "ZLS"->LS (at the end)
97	G4QInteractionVector theInteractions; // A vector of interactions (taken from the Body)
98	G4QPartonPairVector thePartonPairs; // The parton pairs (taken from the Body)
99	G4int ModelMode=SOFT; // The current model type (taken from the Body)
100	theTargNucleus=G4QNucleus(tZ,tN); // Create theTargNucleus to Move From LS to CM
101	theTargNucleus.InitByPDG(tPDG); // Reinit the Nucleus for the new Attempt?
102	theTargNucleus.Init3D(); // 3D-initialisation(nucleons)ofTheTGNucleusClone
103	#ifdef debug
104	G4cout<<"G4QIonIonCollision::Constr:TargNuclMom="<<theTargNucleus.Get4Momentum()<<G4endl;
105	#endif
106	theProjNucleus=G4QNucleus(pZ,pN); // Create theProjNucleus to Move From LS to CM
107	theProjNucleus.InitByPDG(pPDG); // Reinit the Nucleus for the new Attempt?
108	theProjNucleus.Init3D(); // 3D-initialisation(nucleons)ofThePrNucleusClone
109	#ifdef debug
110	G4cout<<"G4QIonIonCollision::Constr:ProjNuclMom="<<theProjNucleus.Get4Momentum()<<G4endl;
111	#endif
112	#ifdef edebug
113	G4LorentzVector sumP1=theProjNucleus.GetNucleons4Momentum();// Sum ofPrNucleons4M inRotLS
114	G4LorentzVector sumT1=theTargNucleus.GetNucleons4Momentum();// Sum ofTgNucleons4M inRotLS
115	G4cout<<"-EMC-G4QIonIonCollision::Construct: pNuc4M="<<sumP1<<", tNuc4M="<<sumT1<<G4endl;
116	#endif
117	G4ThreeVector theCMVelocity(0.,0.,0.); // Prototype of the "CM" velocity
118	G4ThreeVector theALVelocity(0.,0.,0.); // Prototype of "CMAntiLab" velocity
119	// Very important! Here the projectile 4M is recalculated in the ZLS (previously in LS)
120	proj4M = pNucleus.Get4Momentum(); // 4-mom of theProjectileHadron inZLS
121	G4double pz_per_projectile = proj4M.pz()/pA; // Momentum per nucleon
122	G4double e_per_projectile = proj4M.e()/pA+targM/tA;// Add MassOfTargetNucleon
123	G4double vz = pz_per_projectile/e_per_projectile; // Speed of the "oneNuclenCM"
124	#ifdef debug
125	G4cout<<"G4QIonIonCollision::Construct: (ÐÑÐµÐ¯)Projectile4M="<<proj4M<<", Vz="<<vz
126	<<", pA="<<pA<<", pE="<<e_per_projectile<<G4endl;
127	#endif
128	theCMVelocity.setZ(vz); // CM (w/r to one nucleon) velosity
129	theProjNucleus.DoLorentzBoost(-theCMVelocity); // BoostProjNucleus to "CM"
130	theTargNucleus.DoLorentzBoost(-theCMVelocity); // BoostProjNucleus to "CM"
131	G4LorentzVector prN4M=theProjNucleus.Get4Momentum();
132	G4double avz=prN4M.pz()/prN4M.e(); // Speed of AntiLabSys in CMS
133	theALVelocity.setZ(avz); // CM (w/r to one nucleon) velosity
134	#ifdef edebug
135	G4LorentzVector sumP2=theProjNucleus.GetNucleons4Momentum();// Sum of Nucleons4M in RotCM
136	G4LorentzVector sumT2=theTargNucleus.GetNucleons4Momentum();// Sum of Nucleons4M in RotCM
137	G4cout<<"-EMC-G4QIonIonCollision::Construct: Boosted, vCM="<<vz<<", vAL="<<avz<<", tN4M="
138	<<sumT2<<", pN4M="<<sumP2<<G4endl;
139	#endif
140	G4int maxCuts = 7; // @@ Can be reduced for low energies
141	//
142	// >>>>>>>>>> Find collisions meeting collision conditions and the NN interaction XS
143	//
144	G4double outerRadius = theProjNucleus.GetOuterRadius()+theTargNucleus.GetOuterRadius();
145	G4QProbability theProbability(2212); // * proj/targ nucleons *
146	// Clean up all previous interactions and reset the counters
147	#ifdef debug
148	G4cout<<"G4QIonIonCollision::Construct: theIntSize="<<theInteractions.size()<<G4endl;
149	#endif
150	G4int attCnt=-1;
151	G4int maxAtt=27;
152	// From here a loop over nucleons of the projectile Nucleus (@@ Z-sorting isn't implem!!)
153	//
154	G4QHadron* pNucleon; // Proto of the Projectile Nucleon pointer
155	G4QHadron* tNucleon; // Proto of the Target Nucleon pointer
156	G4QNucleus curProjNucleus;
157	G4QNucleus curTargNucleus;
158	while(!theInteractions.size() && ++attCnt < maxAtt)// TillAtLeastOneInteractionIsCreated
159	{
160	// std::for_each(theInteractions.begin(),theInteractions.end(), DeleteQInteraction());
161	// Here we need to clean up the ProjNucleon and the TargNucleon in the Interactions !
162	G4int nint=theInteractions.size(); // For the 1st attempt should be zero
163	for(G4int ni=0; ni < nint; ++ni)
164	{
165	G4QInteraction* curInt = theInteractions[ni];
166	delete curInt->GetProjectile();
167	delete curInt->GetTarget();
168	delete curInt;
169	}
170	theInteractions.clear();
171	// Now we need to make a copy of the projectile and the target nuclei with 3D info
172	if(attCnt) // This is not theFirstAttempt->Clean
173	{
174	curProjNucleus.DeleteNucleons();
175	curTargNucleus.DeleteNucleons();
176	}
177	curProjNucleus = theProjNucleus;
178	curTargNucleus = theTargNucleus;
179	// choose random impact parameter
180	std::pair<G4double, G4double> theImpactParameter;
181	theImpactParameter = curProjNucleus.ChooseImpactXandY(outerRadius);
182	G4double impactX = theImpactParameter.first;
183	G4double impactY = theImpactParameter.second;
184	#ifdef debug
185	G4cout<<"G4QIonIonCollision::Con:At#="<<attCnt<<",X="<<impactX<<",Y="<<impactY<<G4endl;
186	#endif
187	curProjNucleus.StartLoop(); // Prepare Loop ovder nucleons
188	G4int pnCnt=0;
189	G4int pnA=curProjNucleus.GetA();
190	#ifdef debu
191	G4cout<<"G4QIonIonCollision::Constr: Before the WHILE over pNucleons,pA="<<pnA<<G4endl;
192	#endif
193	while((pNucleon=curProjNucleus.GetNextNucleon()) && pnCnt < pnA) // @@ can be for LOOP?
194	{
195	++pnCnt;
196	G4QInteractionVector curInteractions; // A temporary vector of interactions
197	G4LorentzVector pNuc4M=pNucleon->Get4Momentum();
198	G4ThreeVector pNucR=pNucleon->GetPosition(); // Position of the pNucleon WRTo pNucl
199	G4double pImpX = impactX+pNucR.x();
200	G4double pImpY = impactY+pNucR.y();
201	G4double prEn=proj4M.e(); // For mesons
202	G4int proB=pNucleus.GetBaryonNumber();
203	if (proB>0) prEn-=pNucleus.GetMass(); // For baryons
204	else if(proB<0) prEn+=mProt; // For anti-baryons
205	#ifdef debug
206	G4double impactUsed = 0.;
207	G4cout<<"G4QIonIonCollision::Construct: estimated energy, prEn="<<prEn<<G4endl;
208	#endif
209	G4int totalCuts = 0;
210	// @@ This is a fake (random) s calculation @@ can be done inside the TARG-while
211	G4int tnA=curTargNucleus.GetA();
212	G4double pImp=std::sqrt(pImpXpImpX+pImpYpImpY);
213	G4double eflen=curTargNucleus.GetThickness(pImp); // EffectiveLength
214	maxEn=eflen*stringTension; // maxAbsorbedEnergy in IndUnits=MeV
215	maxNuc=static_cast<int>(eflen*tubeDensity+0.5);// max #0f involved nuclear nucleons
216	#ifdef debug
217	G4cout<<"G4QIonIonCollision::Con:pE="<<prEn<<" < mE="<<maxEn<<",mN="<<maxNuc<<G4endl;
218	#endif
219	if(prEn < maxEn) // Create DIN interaction and go out
220	{
221	G4QHadron* aProjectile = new G4QHadron(*pNucleon);// Copy selected PrNuc for String
222	curTargNucleus.StartLoop(); // Initialize newSelection forNucleons
223	tNucleon=curTargNucleus.GetNextNucleon(); // Select a nucleon
224	G4QHadron* aTarget = new G4QHadron(*tNucleon);// Copy selected nucleon for String
225	G4QInteraction* anInteraction = new G4QInteraction(aProjectile);
226	anInteraction->SetTarget(aTarget);
227	anInteraction->SetNumberOfDINRCollisions(1); // Consider this interaction as DINR
228	curInteractions.push_back(anInteraction); //--> now the Interaction is not empty
229	curProjNucleus.DoLorentzBoost(-theALVelocity);// Boost theResPrNucleus toRotAntiLab
230	curProjNucleus.SubtractNucleon(pNucleon); // Pointer to the used ProjNucleon
231	curProjNucleus.DoLorentzBoost(theALVelocity);// Boost theResProjNucleus back to CM
232	curTargNucleus.DoLorentzBoost(theCMVelocity);// Boost theResTgNucleus to RotatedLS
233	curTargNucleus.SubtractNucleon(tNucleon); // Pointer to the used TargNucleon
234	curTargNucleus.DoLorentzBoost(-theCMVelocity);// Boost theResTargNucleus back to CM
235	#ifdef debug
236	G4cout<<"G4QIonIonCollision::Construct: DINR interaction is created"<<G4endl;
237	#endif
238	break; // Break the WHILE of the pNucleons
239	}
240	// LOOP over nuclei of the target nucleus to select collisions
241	curTargNucleus.StartLoop(); // To get the same nucleon
242	G4int tnCnt = 0;
243	#ifdef debu
244	G4cout<<"G4QIonIonCollision::Constr: Before WHILE over tNucleons, tA="<<tnA<<G4endl;
245	#endif
246	while((tNucleon=curTargNucleus.GetNextNucleon()) && tnCnt<tnA && totalCuts<maxCuts)
247	{
248	++tnCnt;
249	G4LorentzVector tNuc4M=tNucleon->Get4Momentum();
250	#ifdef debug
251	G4cout<<"G4QIonIonCollision::Constr: OuterR="<<outerRadius<<", mC="<<maxCuts
252	<<", pA="<<curProjNucleus<<", tA="<<curTargNucleus<<G4endl;
253	#endif
254	// Check the reaction threshold
255	G4double s = (tNuc4M + pNuc4M).mag2(); // Squared CM Energy of compound
256	G4double ThresholdMass = pNucleon->GetMass() + tNucleon->GetMass();
257	#ifdef debug
258	G4cout<<"G4QInel::Constr: s="<<s<<", ThreshM="<<sqr(ThresholdMass)<<G4endl;
259	#endif
260	ModelMode = SOFT; // NOT-Diffractive hadronization
261	if (s < 0.) // At ThP=0 is impossible(virtNucl)
262	{
263	G4cerr<<"*G4QInelast::Constr:s="<<s<<",pN4M="<<pNuc4M<<",tN4M="<<tNuc4M<<G4endl;
264	G4Exception("G4QIonIonCollision::Construct:","72",FatalException,"LowEn(NegS)");
265	}
266	if (s < sqr(ThresholdMass)) // --> Only diffractive interaction
267	{
268	#ifdef debug
269	G4cout<<"G4QIonIonCollision::Constr: *OnlyDiffraction*, ThM="<<ThresholdMass
270	<<">sqrt(s)="<<std::sqrt(s)<<" -> only Diffraction is possible"<<G4endl;
271	#endif
272	ModelMode = DIFFRACTIVE;
273	}
274	G4ThreeVector tNucR=tNucleon->GetPosition(); // Position of the tNucleon WRTo tNucl
275	G4double dImpX = pImpX-tNucR.x();
276	G4double dImpY = pImpY-tNucR.y();
277	G4double Distance2=dImpXdImpX+dImpYdImpY;
278	#ifdef sdebug
279	G4cout<<"G4QIonIonCollision::Construct: s="<<s<<", D2="<<Distance2<<G4endl;
280	#endif
281	// Needs to be moved to Probability class @@
282	if(s<=10000.)
283	{
284	G4cout<<"-Warning-G4QIonIonCollision::Construct: s < .01 GeV^2, p4M="
285	<<pNucleon->Get4Momentum()<<",t4M="<<tNucleon->Get4Momentum()<<G4endl;
286	continue; // skip the rest of the targetNucleons
287	}
288	G4double Probability = theProbability.GetPomInelProbability(s, Distance2);// P_INEL
289	// test for inelastic collision
290	#ifdef sdebug
291	G4cout<<"G4QIonIonCollision::Construct: Probubility="<<Probability<<G4endl;
292	#endif
293	G4double rndNumber = G4UniformRand(); // For the printing purpose
294	// ModelMode = DIFFRACTIVE;
295	#ifdef sdebug
296	G4cout<<"G4QIonIonCollision::Construct: NLOOP prob="<<Probability<<", rndm="
297	<<rndNumber<<", d="<<std::sqrt(Distance2)<<G4endl;
298	#endif
299	if (Probability > rndNumber) // Inelastic (diffractive or soft) interaction (JOB)
300	{
301	G4QHadron* aProjectile = new G4QHadron(*pNucleon);// Copy selected pNuc to String
302	G4QHadron* aTarget = new G4QHadron(*tNucleon);// Copy selected tNucleon to String
303	#ifdef edebug
304	G4cout<<"--->EMC-->G4QIonIonCollision::Construct: TargNucleon is filled, 4M/PDG="
305	<<aTarget->Get4Momentum()<<aTarget->GetPDGCode()<<G4endl;
306	#endif
307	// Now the energy of the nucleons must be updated in CMS
308	curProjNucleus.DoLorentzBoost(-theALVelocity);// Boost theResNucleus toRotatedLS
309	curProjNucleus.SubtractNucleon(pNucleon); // Pointer to the used nucleon
310	curProjNucleus.DoLorentzBoost(theALVelocity); // Boost theResNucleus back to CM
311	curTargNucleus.DoLorentzBoost(theCMVelocity); // Boost theResNucleus toRotatedLS
312	curTargNucleus.SubtractNucleon(tNucleon); // Pointer to the used nucleon
313	curTargNucleus.DoLorentzBoost(-theCMVelocity);// Boost theResNucleus back to CM
314	if((theProbability.GetPomDiffProbability(s,Distance2)/Probability >
315	G4UniformRand() && ModelMode==SOFT ) \|\| ModelMode==DIFFRACTIVE)
316	{
317	// ------------->>>> diffractive interaction @@ IsSingleDiffractive called once
318	if(IsSingleDiffractive()) ExciteSingDiffParticipants(aProjectile, aTarget);
319	else ExciteDiffParticipants(aProjectile, aTarget);
320	G4QInteraction* anInteraction = new G4QInteraction(aProjectile);
321	anInteraction->SetTarget(aTarget);
322	anInteraction->SetNumberOfDiffractiveCollisions(1); // Why not increment? M.K.?
323	curInteractions.push_back(anInteraction);//--> now theInteractions not empty
324	// @@ Why not breake the NLOOP, if only one diffractive can happend?
325	totalCuts++; // UpdateOfNucleons in't necessary
326	#ifdef debug
327	G4cout<<"G4QIonIonCollision::Constr:NLOOP DiffInteract,tC="<<totalCuts<<G4endl;
328	#endif
329	}
330	else
331	{
332	// -------------->>>>> nondiffractive = soft interaction
333	// sample nCut+1 (cut Pomerons) pairs of strings can be produced
334	G4int nCut; // Result in a chosen number of cuts
335	G4double* running = new G4double[nCutMax];// @@ This limits the max cuts
336	for(nCut = 0; nCut < nCutMax; nCut++) // Calculates multiCut probabilities
337	{
338	running[nCut]= theProbability.GetCutPomProbability(s, Distance2, nCut+1);
339	if(nCut) running[nCut] += running[nCut-1];// Sum up with the previous one
340	}
341	G4double random = running[nCutMax-1]*G4UniformRand();
342	for(nCut = 0; nCut < nCutMax; nCut++) if(running[nCut] > random) break; // tNuc
343	delete [] running;
344	#ifdef debug
345	G4cout<<"G4QIonIonCollision::Construct: NLOOP-Soft Chosen nCut="<<nCut<<G4endl;
346	#endif
347	// @@ If nCut>0 interaction with a few nucleons is possible
348	// @@ nCut is found with big efforts and now nCut=0 ?? M.K. ?? !!
349	//nCut = 0; // @@ in original code ?? @@
350	aTarget->IncrementCollisionCount(nCut+1); // @@ What about multyNucleon target?
351	aProjectile->IncrementCollisionCount(nCut+1);
352	G4QInteraction* anInteraction = new G4QInteraction(aProjectile);
353	anInteraction->SetTarget(aTarget);
354	anInteraction->SetNumberOfSoftCollisions(nCut+1);
355	curInteractions.push_back(anInteraction); // Now curInteractions are not empty
356	totalCuts += nCut+1;
357	#ifdef debug
358	G4cout<<"G4QIonIonCollision::Constr:NLOOP SoftInteract,tC="<<totalCuts<<G4endl;
359	impactUsed=Distance2;
360	#endif
361	}
362	}// Probability selection
363	} // End of While over target nucleons
364	G4int nCurInt=curInteractions.size();
365	for(G4int ni=0; ni < nCurInt; ++ni) theInteractions.push_back(curInteractions[ni]);
366	curInteractions.clear(); // Probably, not necessary...
367	} // End of WHILE over projectile nucleons
368	} // End of WHILE over attempts to find at least one nucleus-nucleus interaction
369	theProjNucleus.DeleteNucleons();
370	theTargNucleus.DeleteNucleons();
371	theProjNucleus = curProjNucleus;
372	theTargNucleus = curTargNucleus;
373	curProjNucleus.DeleteNucleons();
374	curTargNucleus.DeleteNucleons();
375	G4int nInt=theInteractions.size();
376	#ifdef debug
377	G4cout<<"G4QIonIonCollision::Constr: #ofInteractions = "<<nInt<<", #OfDINR = "
378	<<theInteractions[0]->GetNumberOfDINRCollisions()<<G4endl;
379	#endif
380	if(!nInt \|\| (nInt==1 && theInteractions[0]->GetNumberOfDINRCollisions()==1)) // QFreeInel
381	{
382	G4QHadron* aTarget;
383	G4QHadron* aProjectile;
384	if(nInt) // Take Targ/Proj from the Interaction
385	{
386	aTarget=theInteractions[0]->GetTarget();
387	aProjectile=theInteractions[0]->GetProjectile();
388	delete theInteractions[0];
389	theInteractions.clear();
390	}
391	else // Create a new target nucleon (?)
392	{
393	theProjNucleus.StartLoop(); // To get the same nucleon from projec
394	pNucleon=theProjNucleus.GetNextNucleon(); // Get theNucleon to create projectNuc
395	aProjectile = new G4QHadron(*pNucleon); // Copy selected pNucleon for String
396	theProjNucleus.DoLorentzBoost(-theALVelocity); // Boost theResProjNucleus toRotatedLS
397	theProjNucleus.SubtractNucleon(pNucleon); // Pointer to SelProjNucleon to delete
398	theProjNucleus.DoLorentzBoost(theALVelocity); // Boost theResProNucleus back to CMS
399	theTargNucleus.StartLoop(); // To get the same nucleon from target
400	tNucleon=theTargNucleus.GetNextNucleon(); // Get theNucleon to create targetNucl
401	aTarget = new G4QHadron(*tNucleon); // Copy selected tNucleon for String
402	theTargNucleus.DoLorentzBoost(theCMVelocity); // Boost theResTargNucleus toRotatedLS
403	theTargNucleus.SubtractNucleon(tNucleon); // Pointer to SelTargNucleon to delete
404	theTargNucleus.DoLorentzBoost(-theCMVelocity); // Boost theResTargNucleus back to CMS
405	}
406	G4QContent QQC=aTarget->GetQC()+aProjectile->GetQC(); // QContent of the compound
407	G4LorentzVector Q4M=aTarget->Get4Momentum()+aProjectile->Get4Momentum(); // 4-mom of Q
408	delete aTarget;
409	delete aProjectile;
410	// @@ 4-Mom should be converted to LS for the TargQuasmon and to AL for the ProjQuasmon
411	Q4M.boost(theCMVelocity);
412	Q4M=toLab*Q4M;
413	G4Quasmon* stringQuasmon = new G4Quasmon(QQC, Q4M);
414	theQuasmons.push_back(stringQuasmon);
415	theTargNucleus.DoLorentzBoost(theCMVelocity); // BoostTheResidualNucleus toRotatedLS
416	theTargNucleus.DoLorentzRotation(toLab); // Recove Z-direction in LS ("LS"->LS)
417	return;
418	}
419	//
420	// ------------------ now build the parton pairs for the strings ------------------
421	//
422	for(G4int i=0; i<nInt; i++)
423	{
424	theInteractions[i]->SplitHadrons();
425	#ifdef edebug
426	G4QHadron* projH=theInteractions[i]->GetProjectile(); // Projectile of theInteraction
427	G4QHadron* targH=theInteractions[i]->GetTarget(); // Target of the Interaction
428	G4LorentzVector pSP(0.,0.,0.,0.); // Sum of parton's 4mom's for proj
429	G4LorentzVector tSP(0.,0.,0.,0.); // Sum of parton's 4mom's for proj
430	std::list<G4QParton*> projCP=projH->GetColor(); // Pointers to proj Color-partons
431	std::list<G4QParton*> projAC=projH->GetAntiColor();// PointersTo projAntiColorPartons
432	std::list<G4QParton*> targCP=targH->GetColor(); // Pointers to targ Color-partons
433	std::list<G4QParton*> targAC=targH->GetAntiColor();// PointersTo targAntiColorPartons
434	std::list<G4QParton*>::iterator picp = projCP.begin();
435	std::list<G4QParton*>::iterator pecp = projCP.end();
436	std::list<G4QParton*>::iterator piac = projAC.begin();
437	std::list<G4QParton*>::iterator peac = projAC.end();
438	std::list<G4QParton*>::iterator ticp = targCP.begin();
439	std::list<G4QParton*>::iterator tecp = targCP.end();
440	std::list<G4QParton*>::iterator tiac = targAC.begin();
441	std::list<G4QParton*>::iterator teac = targAC.end();
442	for(; picp!=pecp&& piac!=peac&& ticp!=tecp&& tiac!=teac; ++picp,++piac,++ticp,++tiac)
443	{
444	pSP+=(*picp)->Get4Momentum();
445	pSP+=(*piac)->Get4Momentum();
446	tSP+=(*ticp)->Get4Momentum();
447	tSP+=(*tiac)->Get4Momentum();
448	}
449	G4cout<<"-EMC-G4QIonIonCollision::Construct: Interaction#"<<i<<",dP4M="
450	<<projH->Get4Momentum()-pSP<<",dT4M="<<targH->Get4Momentum()-tSP<<G4endl;
451	#endif
452	}
453	//
454	// >>>>>>>> make soft collisions (ordering is vital)
455	//
456	G4QInteractionVector::iterator it;
457	#ifdef debug
458	G4cout<<"G4QIonIonCollision::Constr: Creation ofSoftCollisionPartonPair STARTS"<<G4endl;
459	#endif
460	G4bool rep=true;
461	while(rep && theInteractions.size())
462	{
463	for(it = theInteractions.begin(); it != theInteractions.end(); ++it)
464	{
465	G4QInteraction* anIniteraction = *it;
466	G4QPartonPair* aPair=0;
467	G4int nSoftCollisions = anIniteraction->GetNumberOfSoftCollisions();
468	#ifdef debug
469	G4cout<<"G4QIonIonCollision::Construct: #0f SoftCollisions ="<<nSoftCollisions<<G4endl;
470	#endif
471	if (nSoftCollisions)
472	{
473	G4QHadron* pProjectile = anIniteraction->GetProjectile();
474	G4QHadron* pTarget = anIniteraction->GetTarget();
475	for (G4int j = 0; j < nSoftCollisions; j++)
476	{
477	aPair = new G4QPartonPair(pTarget->GetNextParton(),
478	pProjectile->GetNextAntiParton(),
479	G4QPartonPair::SOFT, G4QPartonPair::TARGET);
480	thePartonPairs.push_back(aPair); // A target pair (Why TAGRET?)
481	aPair = new G4QPartonPair(pProjectile->GetNextParton(),
482	pTarget->GetNextAntiParton(),
483	G4QPartonPair::SOFT, G4QPartonPair::PROJECTILE);
484	thePartonPairs.push_back(aPair); // A projectile pair (Why Projectile?)
485	#ifdef debug
486	G4cout<<"--->G4QIonIonCollision::Constr: SOFT, 2 parton pairs are filled"<<G4endl;
487	#endif
488	}
489	delete *it;
490	it=theInteractions.erase(it); // SoftInteractions're converted&erased
491	if( it != theInteractions.begin() )// To avoid going below begin() (for Windows)
492	{
493	it--;
494	rep=false;
495	}
496	else
497	{
498	rep=true;
499	break;
500	}
501	}
502	else rep=false;
503	}
504	}
505	#ifdef debug
506	G4cout<<"G4QIonIonCollision::Constr: -> Parton pairs for SOFT strings are made"<<G4endl;
507	#endif
508	//
509	// >>>>>>>>>>>>>>> make the rest as the diffractive interactions
510	//
511	for(unsigned i = 0; i < theInteractions.size(); i++) // Interactions are reduced bySoft
512	{
513	// The double or single diffraction is defined by the presonce of proj/targ partons
514	G4QInteraction* anIniteraction = theInteractions[i];
515	G4QPartonPair* aPartonPair;
516	#ifdef debug
517	G4cout<<"G4QIonIonCollision::Constr: CreationOfDiffractivePartonPairs, i="<<i<<G4endl;
518	#endif
519	// the projectile diffraction parton pair is created first
520	G4QHadron* aProjectile = anIniteraction->GetProjectile();
521	G4QParton* aParton = aProjectile->GetNextParton();
522	if (aParton)
523	{
524	aPartonPair = new G4QPartonPair(aParton, aProjectile->GetNextAntiParton(),
525	G4QPartonPair::DIFFRACTIVE,
526	G4QPartonPair::PROJECTILE);
527	thePartonPairs.push_back(aPartonPair);
528	#ifdef debug
529	G4cout<<"G4QIonIonCollision::Constr: proj Diffractive PartonPair is filled"<<G4endl;
530	#endif
531	}
532	// then the target diffraction parton pair is created
533	G4QHadron* aTarget = anIniteraction->GetTarget();
534	aParton = aTarget->GetNextParton();
535	if (aParton)
536	{
537	aPartonPair = new G4QPartonPair(aParton, aTarget->GetNextAntiParton(),
538	G4QPartonPair::DIFFRACTIVE, G4QPartonPair::TARGET);
539	thePartonPairs.push_back(aPartonPair);
540	#ifdef debug
541	G4cout<<"G4QIonIonCollision::Constr: target Diffractive PartonPair's filled"<<G4endl;
542	#endif
543	}
544	}
545	#ifdef debug
546	G4cout<<"G4QIonIonCollision::Construct: DiffractivePartonPairs are created"<<G4endl;
547	#endif
548	//
549	// >>>>>>>>>>>>>> clean-up Interactions, if necessary
550	//
551	std::for_each(theInteractions.begin(),theInteractions.end(), DeleteQInteraction());
552	theInteractions.clear();
553	#ifdef debug
554	G4cout<<"G4QIonIonCollision::Construct: Temporary objects are cleaned up"<<G4endl;
555	#endif
556	// This function prepares theBoost for transformation of secondaries to LS (-ProjRot!)
557	theProjNucleus.DoLorentzBoost(theCMVelocity);// Boost theResidualProjNucleus to RotatedLS
558	theTargNucleus.DoLorentzBoost(theCMVelocity);// Boost theResidualTargNucleus to RotatedLS
559	// @@ Nucleus isn't completely in LS, it needs the toZ (-ProjRot) rotation to consE/M
560	#ifdef debug
561	G4cout<<">>>>>>>>>>>>G4QIonIonCollision::Construct: >>>>>> Strings are created "<<G4endl;
562	#endif
563	G4QPartonPair* aPair;
564	G4QString* aString=0;
565	while(thePartonPairs.size()) // @@ At present noDifference in stringBuild (? M.K.)
566	{
567	aPair = thePartonPairs.back(); // Get the parton pair
568	thePartonPairs.pop_back(); // Clean up thePartonPairPointer in the Vector
569	#ifdef debug
570	G4cout<<"G4QIonIonCollision::Construct:StringType="<<aPair->GetCollisionType()<<G4endl;
571	#endif
572	aString= new G4QString(aPair);
573	#ifdef debug
574	G4cout<<"G4QIonIonCollision::Construct:NewString4M="<<aString->Get4Momentum()<<G4endl;
575	#endif
576	aString->Boost(theCMVelocity); // ! Strings are moved to ZLS when pushed !
577	strings.push_back(aString);
578	stringsInitted=true;
579	delete aPair;
580	} // End of the String Creation LOOP
581	#ifdef edebug
582	G4LorentzVector sum=theProjNucleus.Get4Momentum()+theTargNucleus.Get4Momentum(); // in LS
583	G4int rChg=totChg-theProjNucleus.GetZ()-theTargNucleus.GetZ();
584	G4int rBaN=totBaN-theProjNucleus.GetA()-theTargNucleus.GetA();
585	G4int nStrings=strings.size();
586	G4cout<<"-EMC-G4QIonIonCollision::Constr:#ofString="<<nStrings<<",tNuc4M="<<sum<<G4endl;
587	for(G4int i=0; i<nStrings; i++)
588	{
589	G4QString* prString=strings[i];
590	G4LorentzVector strI4M=prString->Get4Momentum();
591	sum+=strI4M;
592	G4int sChg=prString->GetCharge();
593	G4int sBaN=prString->GetBaryonNumber();
594	G4int LPDG=prString->GetLeftParton()->GetPDGCode();
595	G4int RPDG=prString->GetRightParton()->GetPDGCode();
596	G4QContent LQC =prString->GetLeftParton()->GetQC();
597	G4QContent RQC =prString->GetRightParton()->GetQC();
598	rChg-=sChg;
599	rBaN-=sBaN;
600	G4cout<<"-EMC-G4QIonIonCollision::Construct: String#"<<i<<", 4M="<<strI4M<<", LPDG="
601	<<LPDG<<LQC<<",RPDG="<<RPDG<<RQC<<", Ch="<<sChg<<", BN="<<sBaN<<G4endl;
602	}
603	G4cout<<"-EMC-G4QInel::Constr: r4M="<<sum-totZLS4M<<", rC="<<rChg<<", rB="<<rBaN<<G4endl;
604	#endif
605	if(!stringsInitted)
606	{
607	G4cerr<<"***G4QIonIonCollision::Construct: No strings are created ***"<<G4endl;
608	G4Exception("G4QIonIonCollision::Construct:","72",FatalException,"No Strings created");
609	}
610	#ifdef debug
611	G4cout<<"G4QIonIonCollision::Constr:BeforeRotation, #0fStrings="<<strings.size()<<G4endl;
612	#endif
613	//
614	// ---------------- At this point the strings must be already created in "LS" -----------
615	//
616	for(unsigned astring=0; astring < strings.size(); astring++)
617	strings[astring]->LorentzRotate(toLab); // Recove Z-direction in LS ("LS"->LS)
618	theProjNucleus.DoLorentzRotation(toLab); // Recove Z-dir in LS ("LS"->LS) for ProjNucleus
619	theTargNucleus.DoLorentzRotation(toLab); // Recove Z-dir in LS ("LS"->LS) for TargNucleus
620	// Now everything is in LS system
621	#ifdef edebug
622	G4LorentzVector sm=theProjNucleus.Get4Momentum()+theTargNucleus.Get4Momentum();// NucInLS
623	G4int rCg=totChg-theProjNucleus.GetZ()-theTargNucleus.GetZ();
624	G4int rBC=totBaN-theProjNucleus.GetA()-theTargNucleus.GetA();
625	G4int nStrs=strings.size();
626	G4cout<<"-EMCLS-G4QIonIonCollision::Constr:#ofS="<<nStrings<<",Nuc4M(E=M)="<<sum<<G4endl;
627	for(G4int i=0; i<nStrs; i++)
628	{
629	G4LorentzVector strI4M=strings[i]->Get4Momentum();
630	sm+=strI4M;
631	G4int sChg=strings[i]->GetCharge();
632	rCg-=sChg;
633	G4int sBaN=strings[i]->GetBaryonNumber();
634	rBC-=sBaN;
635	G4cout<<"-EMCLS-G4QInel::Construct:String#"<<i<<",4M="<<strI4M<<strI4M.m()<<",Charge="
636	<<sChg<<",BaryN="<<sBaN<<G4endl;
637	}
638	G4cout<<"-EMCLS-...G4QInel::Constr: r4M="<<sm-totLS4M<<", rC="<<rCg<<",rB="<<rBC<<G4endl;
639	#endif
640	//
641	// --- Strings are created, but we should try to get rid of negative mass strings -----
642	//
643	SwapPartons();
644	//
645	// --- Strings are created, but we should get rid of too light strings (Mmin+MPi0) -----
646	//
647	G4int problem=0; // 0="no problem", incremented by ASIS
648	G4QStringVector::iterator ist;
649	G4bool con=true;
650	while(con && strings.size())
651	{
652	for(ist = strings.begin(); ist < strings.end(); ++ist)
653	{
654	G4bool bad=true;
655	G4LorentzVector cS4M=(*ist)->Get4Momentum();
656	G4double cSM2=cS4M.m2(); // Squared mass of the String
657	G4QParton* cLeft=(*ist)->GetLeftParton();
658	G4QParton* cRight=(*ist)->GetRightParton();
659	G4int cLT=cLeft->GetType();
660	G4int cRT=cRight->GetType();
661	G4int cLPDG=cLeft->GetPDGCode();
662	G4int cRPDG=cRight->GetPDGCode();
663	G4int aLPDG=0;
664	G4int aRPDG=0;
665	if (cLPDG > 7) aLPDG= cLPDG/100;
666	else if(cLPDG <-7) aLPDG=(-cLPDG)/100;
667	if (cRPDG > 7) aRPDG= cRPDG/100;
668	else if(cRPDG <-7) aRPDG=(-cRPDG)/100;
669	G4int L1=0;
670	G4int L2=0;
671	if(aLPDG)
672	{
673	L1=aLPDG/10;
674	L2=aLPDG%10;
675	}
676	G4int R1=0;
677	G4int R2=0;
678	if(aRPDG)
679	{
680	R1=aRPDG/10;
681	R2=aRPDG%10;
682	}
683	G4double cSM=cSM2;
684	if(cSM2>0.) cSM=std::sqrt(cSM2);
685	#ifdef debug
686	G4cout<<"G4QIonIonCollision::Construct: Needs Fusion? cLPDG="<<cLPDG<<",cRPDG="<<cRPDG
687	<<",cM(cM2 if neg)="<<cSM<<G4endl;
688	#endif
689	if(cSM>0.) // Mass can be calculated
690	{
691	G4bool single=true;
692	G4double miM=0.; // Proto of the Min HadronString Mass
693	if(cLT==2 && cRT==2)
694	{
695	if(L1!=R1 && L1!=R2 && L2!=R1 && L2!=R2) // Unreducable DiQ-aDiQ
696	{
697	single=false;
698	G4QPDGCode tmp;
699	std::pair<G4int,G4int> paB=tmp.MakeTwoBaryons(L1, L2, R1, R2);
700	miM=G4QPDGCode(paB.first).GetMass()+G4QPDGCode(paB.second).GetMass();
701	}
702	}
703	if(single) miM=G4QPDGCode((*ist)->GetQC().GetSPDGCode()).GetMass() + mPi0;//MinHSMass
704	//if(single) miM=G4QPDGCode((*ist)->GetQC().GetSPDGCode()).GetMass();//MinHSMass
705	#ifdef debug
706	G4cout<<"G4QInel::Const:*IsItGood? realM="<<std::sqrt(cSM2)<<" > GSM="<<miM<<G4endl;
707	#endif
708	if(std::sqrt(cSM2) > miM) bad=false; // String is OK
709	}
710	if(bad) // String should be merged with others
711	{
712	#ifdef debug
713	G4cout<<"G4QInel::Const:TryFuse,L1="<<L1<<",L2="<<L2<<",R1="<<R1<<",R2="<<R2<<G4endl;
714	#endif
715	G4int cST=cLT+cRT;
716	G4double excess=-DBL_MAX; // The value to be maximized excess M
717	G4double maxiM2=-DBL_MAX; // The value to be maximized M2
718	G4QStringVector::iterator sst; // Selected partner string
719	G4QStringVector::iterator pst;
720	G4int sLPDG=0; // selectedLeft (like inStringPartner)
721	G4int sRPDG=0; // selectedRight(like inStringPartner)
722	G4int sOrd=0; // selected Order of PartonFusion
723	G4bool minC=true; // for the case when M2<0
724	if(cSM2>0.) minC=false; // If M2>0 already don'tSearchFor M2>0
725	for(pst = strings.begin(); pst < strings.end(); pst++) if(pst != ist)
726	{
727	G4LorentzVector pS4M=(*pst)->Get4Momentum()+cS4M; // Summed 4-momentum
728	G4int nLPDG=0; // new Left (like in theStringPartner)
729	G4int nRPDG=0; // new Right(like in theStringPartner)
730	G4double pExcess=-DBL_MAX; // Prototype of the excess
731	G4double pSM2=pS4M.m2(); // Squared mass of the Fused Strings
732	#ifdef debug
733	G4cout<<"->G4QIonIonCollision::Construct: sum4M="<<pS4M<<", M2="<<pSM2<<G4endl;
734	#endif
735	//if(pSM2>0.) // The partner can be a candidate
736	//{
737	G4QParton* pLeft=(*pst)->GetLeftParton();
738	G4QParton* pRight=(*pst)->GetRightParton();
739	G4int pLT=pLeft->GetType();
740	G4int pRT=pRight->GetType();
741	G4int pLPDG=pLeft->GetPDGCode();
742	G4int pRPDG=pRight->GetPDGCode();
743	G4int LPDG=0;
744	G4int RPDG=0;
745	if (pLPDG > 7) LPDG= pLPDG/100;
746	else if(pLPDG <-7) LPDG=(-pLPDG)/100;
747	if (pRPDG > 7) RPDG= pRPDG/100;
748	else if(pRPDG <-7) RPDG=(-pRPDG)/100;
749	G4int pL1=0;
750	G4int pL2=0;
751	if(LPDG)
752	{
753	pL1=LPDG/10;
754	pL2=LPDG%10;
755	}
756	G4int pR1=0;
757	G4int pR2=0;
758	if(RPDG)
759	{
760	pR1=RPDG/10;
761	pR2=RPDG%10;
762	}
763	G4int pST=pLT+pRT;
764	#ifdef debug
765	G4cout<<"G4QInel::Construct: Partner/w pLPDG="<<pLPDG<<", pRPDG="<<pRPDG<<", pM2="
766	<<pSM2<<G4endl;
767	#endif
768	// Different fromCompactAlrorithm ofStringFusionAfterDecay (no DiQaDiQ reduction)
769	G4int tf=0; // Type combination flag
770	G4int af=0; // Annihilatio combination flag
771	if (cST==2 && pST==2) // QaQ + QaQ = DiQaDiQ (always)
772	{
773	tf=1;
774	af=1;
775	}
776	else if(cST==2 && pST==3) // QaQ + QDiQ/aQaDiQ = QDiQ/aQaDiQ (s)
777	{
778	tf=2;
779	if (pLPDG > 7 &&
780	( (cLPDG<0 && (-cLPDG==pL1 \|\| -cLPDG==pL2 \|\| -cLPDG==pRPDG) ) \|\|
781	(cRPDG<0 && (-cRPDG==pL1 \|\| -cRPDG==pL2 \|\| -cRPDG==pRPDG) )
782	)
783	) af=1;
784	else if(pRPDG > 7 &&
785	( (cLPDG<0 && (-cLPDG==pR1 \|\| -cLPDG==pR2 \|\| -cLPDG==pLPDG) ) \|\|
786	(cRPDG<0 && (-cRPDG==pR1 \|\| -cRPDG==pR2 \|\| -cRPDG==pLPDG) )
787	)
788	) af=2;
789	else if(pLPDG <-7 &&
790	( (cLPDG>0 && ( cLPDG==pL1 \|\| cLPDG==pL2 \|\| cLPDG==-pRPDG) ) \|\|
791	(cRPDG>0 && ( cRPDG==pL1 \|\| cRPDG==pL2 \|\| cRPDG==-pRPDG) )
792	)
793	) af=3;
794	else if(pRPDG <-7 &&
795	( (cLPDG>0 && ( cLPDG==pR1 \|\| cLPDG==pR2 \|\| cLPDG==-pLPDG) ) \|\|
796	(cRPDG>0 && ( cRPDG==pR1 \|\| cRPDG==pR2 \|\| cRPDG==-pLPDG) )
797	)
798	) af=4;
799	#ifdef debug
800	else G4cout<<"G4QIonIonCollision::Constr:2(QaQ+QDiQ/aQaDiQ) Can't fuse"<<G4endl;
801	#endif
802	}
803	else if(cST==3 && pST==2) // QDiQ/aQaDiQ + QaQ = QDiQ/aQaDiQ (s)
804	{
805	tf=3;
806	if (cLPDG > 7 &&
807	( (pLPDG<0 && (-pLPDG==L1 \|\| -pLPDG==L2) ) \|\|
808	(pRPDG<0 && (-pRPDG==L1 \|\| -pRPDG==L2) )
809	)
810	) af=1;
811	else if(cRPDG > 7 &&
812	( (pLPDG<0 && (-pLPDG==R1 \|\| -pLPDG==R2) ) \|\|
813	(pRPDG<0 && (-pRPDG==R1 \|\| -pRPDG==R2) )
814	)
815	) af=2;
816	else if(cLPDG <-7 &&
817	( (pLPDG>0 && ( pLPDG==L1 \|\| pLPDG==L2) ) \|\|
818	(pRPDG>0 && ( pRPDG==L1 \|\| pRPDG==L2) )
819	)
820	) af=3;
821	else if(cRPDG <-7 &&
822	( (pLPDG>0 && ( pLPDG==R1 \|\| pLPDG==R2) ) \|\|
823	(pRPDG>0 && ( pRPDG==R1 \|\| pRPDG==R2) )
824	)
825	) af=4;
826	#ifdef debug
827	else G4cout<<"G4QIonIonCollision::Constr:3(QDiQ/aQaDiQ+QaQ) Can't fuse"<<G4endl;
828	#endif
829	}
830	else if(cST==2 && pST==4) // QaQ + aDiQDiQ = QaQ (double)
831	{
832	tf=4;
833	if (pLPDG > 7) // pRPDG <-7
834	{
835	if ( (-cLPDG==pL1 \|\| -cLPDG==pL2) && (cRPDG==pR1 \|\| cRPDG==pR2) ) af=1;
836	else if( (-cRPDG==pL1 \|\| -cRPDG==pL2) && (cLPDG==pR1 \|\| cLPDG==pR2) ) af=2;
837	}
838	else if(pRPDG > 7) // pLPDG <-7
839	{
840	if ( (-cRPDG==pR1 \|\| -cRPDG==pR2) && (cLPDG==pL1 \|\| cLPDG==pL2) ) af=3;
841	else if( (-cLPDG==pR1 \|\| -cLPDG==pR2) && (cRPDG==pL1 \|\| cRPDG==pL2) ) af=4;
842	}
843	#ifdef debug
844	else G4cout<<"-G4QIonIonCollision::Constr: 4 (QaQ+aQDiQDiQ) Can't fuse"<<G4endl;
845	#endif
846	}
847	else if(cST==4 && pST==2) // aDiQDiQ + QaQ = QaQ (double)
848	{
849	tf=5;
850	if (cLPDG > 7) // cRPDG<-7
851	{
852	if ( (-pLPDG==L1 \|\| -pLPDG==L2) && (pRPDG==R1 \|\| pRPDG==R2) ) af=1;
853	else if( (-pRPDG==L1 \|\| -pRPDG==L2) && (pLPDG==R1 \|\| pLPDG==R2) ) af=2;
854	}
855	else if(cRPDG > 7) // cLPDG<-7
856	{
857	if ( (-pRPDG==R1 \|\| -pRPDG==R2) && (pLPDG==L1 \|\| pLPDG==L2) ) af=3;
858	else if( (-pLPDG==R1 \|\| -pLPDG==R2) && (pRPDG==L1 \|\| pRPDG==L2) ) af=4;
859	}
860	#ifdef debug
861	else G4cout<<"-G4QIonIonCollision::Constr: 5 (aQDiQDiQ+QaQ) Can't fuse"<<G4endl;
862	#endif
863	}
864	else if(cST==3 && pST==3) // QDiQ + aQaDiQ = QaQ (double)
865	{
866	tf=6;
867	if(pLPDG > 7)
868	{
869	if (cLPDG<-7 && (-cRPDG==pL1 \|\| -cRPDG==pL2) && (pRPDG==L1 \|\| pRPDG==L2))
870	af=1;
871	else if(cRPDG<-7 && (-cLPDG==pL1 \|\| -cLPDG==pL2) && (pRPDG==R1 \|\| pRPDG==R2))
872	af=2;
873	}
874	else if(pRPDG > 7)
875	{
876	if (cLPDG<-7 && (-cRPDG==pR1 \|\| -cRPDG==pR2) && (pLPDG==L1 \|\| pLPDG==L2))
877	af=3;
878	else if(cRPDG<-7 && (-cLPDG==pR1 \|\| -cLPDG==pR2) && (pLPDG==R1 \|\| pLPDG==R2))
879	af=4;
880	}
881	else if(cLPDG > 7)
882	{
883	if (pLPDG<-7 && (-pRPDG==L1 \|\| -pRPDG==L2) && (cRPDG==pL1 \|\| cRPDG==pL2))
884	af=5;
885	else if(pRPDG<-7 && (-pLPDG==L1 \|\| -pLPDG==L2) && (cRPDG==pR1 \|\| cRPDG==pR2))
886	af=6;
887	}
888	else if(cRPDG > 7)
889	{
890	if (pLPDG<-7 && (-pRPDG==R1 \|\| -pRPDG==R2) && (cLPDG==pL1 \|\| cLPDG==pL2))
891	af=7;
892	else if(pRPDG<-7 && (-pLPDG==R1 \|\| -pLPDG==R2) && (cLPDG==pR1 \|\| cLPDG==pR2))
893	af=8;
894	}
895	#ifdef debug
896	else G4cout<<"-G4QIonIonCollision::Constr: 6 (QDiQ+aQaDiQ) Can't fuse"<<G4endl;
897	#endif
898	}
899	#ifdef debug
900	G4cout<<"G4QIonIonCollision::Constr: Possibility, tf="<<tf<<",af="<<af<<G4endl;
901	#endif
902	if(tf && af)
903	{
904	// Strings can be fused, update the max excess and remember usefull switches
905	G4int order=0; // LL/RR (1) or LR/RL (2) PartonFusion
906	switch (tf)
907	{
908	case 1: // ------------------------------------> QaQ + QaQ = DiQaDiQ (always)
909	if (cLPDG > 0 && pLPDG > 0)
910	{
911	order= 1;
912	if (cLPDG > pLPDG) nLPDG=cLPDG1000+pLPDG100+1;
913	else if(cLPDG < pLPDG) nLPDG=pLPDG1000+cLPDG100+1;
914	else nLPDG=pLPDG1000+cLPDG100+3;
915	if (cRPDG < pRPDG) nRPDG=cRPDG1000+pRPDG100-1;
916	else if(cRPDG > pRPDG) nRPDG=pRPDG1000+cRPDG100-1;
917	else nRPDG=pRPDG1000+cRPDG100-3;
918	}
919	else if(cLPDG < 0 && pLPDG < 0)
920	{
921	order= 1;
922	if (cRPDG > pRPDG) nRPDG=cRPDG1000+pRPDG100+1;
923	else if(cRPDG < pRPDG) nRPDG=pRPDG1000+cRPDG100+1;
924	else nRPDG=pRPDG1000+cRPDG100+3;
925	if (cLPDG < pLPDG) nLPDG=cLPDG1000+pLPDG100-1;
926	else if(cLPDG > pLPDG) nLPDG=pLPDG1000+cLPDG100-1;
927	else nLPDG=pLPDG1000+cLPDG100-3;
928	}
929	else if(cRPDG > 0 && pLPDG > 0)
930	{
931	order=-1;
932	if (cRPDG > pLPDG) nLPDG=cRPDG1000+pLPDG100+1;
933	else if(cRPDG < pLPDG) nLPDG=pLPDG1000+cRPDG100+1;
934	else nLPDG=pLPDG1000+cRPDG100+3;
935	if (cLPDG < pRPDG) nRPDG=cLPDG1000+pRPDG100-1;
936	else if(cLPDG > pRPDG) nRPDG=pRPDG1000+cLPDG100-1;
937	else nRPDG=pRPDG1000+cLPDG100-3;
938	}
939	else if(cRPDG < 0 && pLPDG < 0)
940	{
941	order=-1;
942	if (cLPDG > pRPDG) nRPDG=cLPDG1000+pRPDG100+1;
943	else if(cLPDG < pRPDG) nRPDG=pRPDG1000+cLPDG100+1;
944	else nRPDG=pRPDG1000+cLPDG100+3;
945	if (cRPDG < pLPDG) nLPDG=cRPDG1000+pLPDG100-1;
946	else if(cRPDG > pLPDG) nLPDG=pLPDG1000+cRPDG100-1;
947	else nLPDG=pLPDG1000+cRPDG100-3;
948	}
949	break;
950	case 2: // ------------------------> QaQ + QDiQ/aQaDiQ = QDiQ/aQaDiQ (single)
951	switch (af)
952	{
953	case 1: // ....................... pLPDG > 7
954	if(cLPDG < 0)
955	{
956	order= 1;
957	if(-cLPDG==pRPDG)
958	{
959	nLPDG=pLPDG;
960	nRPDG=cRPDG;
961	}
962	else
963	{
964	if (cRPDG > pRPDG) nRPDG=cRPDG1000+pRPDG100+1;
965	else if(cRPDG < pRPDG) nRPDG=pRPDG1000+cRPDG100+1;
966	else nRPDG=pRPDG1000+cRPDG100+3;
967	if (-cLPDG == pL1) nLPDG=pL2;
968	else nLPDG=pL1; // -cLPDG == pL2
969	}
970	}
971	else // cRPDG < 0
972	{
973	order=-1;
974	if(-cRPDG==pRPDG)
975	{
976	nLPDG=pLPDG;
977	nRPDG=cLPDG;
978	}
979	else
980	{
981	if (cLPDG > pRPDG) nRPDG=cLPDG1000+pRPDG100+1;
982	else if(cLPDG < pRPDG) nRPDG=pRPDG1000+cLPDG100+1;
983	else nRPDG=pRPDG1000+cLPDG100+3;
984	if (-cRPDG == pL1) nLPDG=pL2;
985	else nLPDG=pL1; // -cRPDG == pL2
986	}
987	}
988	break;
989	case 2: // ....................... pRPDG > 7
990	if(cLPDG < 0)
991	{
992	order=-1;
993	if(-cLPDG==pLPDG)
994	{
995	nLPDG=cRPDG;
996	nRPDG=pRPDG;
997	}
998	else
999	{
1000	if (cRPDG > pLPDG) nLPDG=cRPDG1000+pLPDG100+1;
1001	else if(cRPDG < pLPDG) nLPDG=pLPDG1000+cRPDG100+1;
1002	else nLPDG=pLPDG1000+cRPDG100+3;
1003	if (-cLPDG == pR1) nRPDG=pR2;
1004	else nRPDG=pR1; // -cLPDG == pR2
1005	}
1006	}
1007	else // cRPDG < 0
1008	{
1009	order= 1;
1010	if(-cRPDG==pLPDG)
1011	{
1012	nLPDG=cLPDG;
1013	nRPDG=pRPDG;
1014	}
1015	else
1016	{
1017	if (cLPDG > pLPDG) nLPDG=cLPDG1000+pLPDG100+1;
1018	else if(cLPDG < pLPDG) nLPDG=pLPDG1000+cLPDG100+1;
1019	else nLPDG=pLPDG1000+cLPDG100+3;
1020	if (-cRPDG == pR1) nRPDG=pR2;
1021	else nRPDG=pR1; // -cRPDG == pR2
1022	}
1023	}
1024	break;
1025	case 3: // ....................... pLPDG <-7
1026	if(cLPDG > 0)
1027	{
1028	order= 1;
1029	if(cLPDG==-pRPDG)
1030	{
1031	nLPDG=pLPDG;
1032	nRPDG=cRPDG;
1033	}
1034	else
1035	{
1036	if (cRPDG < pRPDG) nRPDG=cRPDG1000+pRPDG100-1;
1037	else if(cRPDG > pRPDG) nRPDG=pRPDG1000+cRPDG100-1;
1038	else nRPDG=pRPDG1000+cRPDG100-3;
1039	if ( cLPDG == pL1) nLPDG=-pL2;
1040	else nLPDG=-pL1; // cLPDG == pL2
1041	}
1042	}
1043	else // cRPDG > 0
1044	{
1045	order=-1;
1046	if(cRPDG==-pRPDG)
1047	{
1048	nLPDG=pLPDG;
1049	nRPDG=cLPDG;
1050	}
1051	else
1052	{
1053	if (cLPDG < pRPDG) nRPDG=cLPDG1000+pRPDG100-1;
1054	else if(cLPDG > pRPDG) nRPDG=pRPDG1000+cLPDG100-1;
1055	else nRPDG=pRPDG1000+cLPDG100-3;
1056	if ( cRPDG == pL1) nLPDG=-pL2;
1057	else nLPDG=-pL1; // cRPDG == pL2
1058	}
1059	}
1060	break;
1061	case 4: // ....................... pRPDG <-7
1062	if(cLPDG > 0)
1063	{
1064	order=-1;
1065	if(cLPDG==-pLPDG)
1066	{
1067	nLPDG=cRPDG;
1068	nRPDG=pRPDG;
1069	}
1070	else
1071	{
1072	if (cRPDG < pLPDG) nLPDG=cRPDG1000+pLPDG100-1;
1073	else if(cRPDG > pLPDG) nLPDG=pLPDG1000+cRPDG100-1;
1074	else nLPDG=pLPDG1000+cRPDG100-3;
1075	if ( cLPDG == pR1) nRPDG=-pR2;
1076	else nRPDG=-pR1; // cLPDG == pR2
1077	}
1078	}
1079	else // cRPDG > 0
1080	{
1081	order= 1;
1082	if(cRPDG==-pLPDG)
1083	{
1084	nLPDG=cLPDG;
1085	nRPDG=pRPDG;
1086	}
1087	else
1088	{
1089	if (cLPDG < pLPDG) nLPDG=cLPDG1000+pLPDG100-1;
1090	else if(cLPDG > pLPDG) nLPDG=pLPDG1000+cLPDG100-1;
1091	else nLPDG=pLPDG1000+cLPDG100-3;
1092	if ( cRPDG == pR1) nRPDG=-pR2;
1093	else nRPDG=-pR1; // cRPDG == pR2
1094	}
1095	}
1096	break;
1097	}
1098	break;
1099	case 3: // ------------------------> QDiQ/aQaDiQ + QaQ = QDiQ/aQaDiQ (single)
1100	switch (af)
1101	{
1102	case 1: // ....................... cLPDG > 7
1103	if(pLPDG < 0)
1104	{
1105	order= 1;
1106	if (pRPDG > cRPDG) nRPDG=pRPDG1000+cRPDG100+1;
1107	else if(pRPDG < cRPDG) nRPDG=cRPDG1000+pRPDG100+1;
1108	else nRPDG=cRPDG1000+pRPDG100+3;
1109	if (-pLPDG == L1) nLPDG=L2;
1110	else nLPDG=L1; // -pLPDG == L2
1111	}
1112	else // pRPDG < 0
1113	{
1114	order=-1;
1115	if (pLPDG > cRPDG) nLPDG=pLPDG1000+cRPDG100+1;
1116	else if(pLPDG < cRPDG) nLPDG=cRPDG1000+pLPDG100+1;
1117	else nLPDG=cRPDG1000+pLPDG100+3;
1118	if (-pRPDG == L1) nRPDG=L2;
1119	else nRPDG=L1; // -pRPDG == L2
1120	}
1121	break;
1122	case 2: // ....................... cRPDG > 7
1123	if(pLPDG < 0)
1124	{
1125	order=-1;
1126	if (pRPDG > cLPDG) nRPDG=pRPDG1000+cLPDG100+1;
1127	else if(pRPDG < cLPDG) nRPDG=cLPDG1000+pRPDG100+1;
1128	else nRPDG=cLPDG1000+pRPDG100+3;
1129	if (-pLPDG == R1) nLPDG=R2;
1130	else nLPDG=R1; // -pLPDG == R2
1131	}
1132	else // pRPDG < 0
1133	{
1134	order= 1;
1135	if (pLPDG > cLPDG) nLPDG=pLPDG1000+cLPDG100+1;
1136	else if(pLPDG < cLPDG) nLPDG=cLPDG1000+pLPDG100+1;
1137	else nLPDG=cLPDG1000+pLPDG100+3;
1138	if (-pRPDG == R1) nRPDG=R2;
1139	else nRPDG=R1; // -pRPDG == R2
1140	}
1141	break;
1142	case 3: // ....................... cLPDG <-7 (cRPDG <0)
1143	if(pLPDG > 0)
1144	{
1145	order= 1;
1146	if (pRPDG < cRPDG) nRPDG=pRPDG1000+cRPDG100-1;
1147	else if(pRPDG > cRPDG) nRPDG=cRPDG1000+pRPDG100-1;
1148	else nRPDG=cRPDG1000+pRPDG100-3;
1149	if ( pLPDG == L1) nLPDG=-L2;
1150	else nLPDG=-L1; // pLPDG == L2
1151	}
1152	else // pRPDG > 0
1153	{
1154	order=-1;
1155	if (pLPDG < cRPDG) nLPDG=pLPDG1000+cRPDG100-1;
1156	else if(pLPDG > cRPDG) nLPDG=cRPDG1000+pLPDG100-1;
1157	else nLPDG=cRPDG1000+pLPDG100-3;
1158	if ( pRPDG == L1) nRPDG=-L2;
1159	else nRPDG=-L1; // pRPDG == L2
1160	}
1161	break;
1162	case 4: // ....................... cRPDG <-7 (cLPDG <0)
1163	if(pLPDG > 0) // pRPDG & cLPDG are anti-quarks
1164	{
1165	order=-1;
1166	if (pRPDG < cLPDG) nRPDG=pRPDG1000+cLPDG100-1;
1167	else if(pRPDG > cLPDG) nRPDG=cLPDG1000+pRPDG100-1;
1168	else nRPDG=cLPDG1000+pRPDG100-3;
1169	if ( pLPDG == R1) nLPDG=-R2;
1170	else nLPDG=-R1; // pLPDG == R2
1171	}
1172	else // pRPDG > 0
1173	{
1174	order= 1;
1175	if (pLPDG < cLPDG) nLPDG=pLPDG1000+cLPDG100-1;
1176	else if(pLPDG > cLPDG) nLPDG=cLPDG1000+pLPDG100-1;
1177	else nLPDG=cLPDG1000+pLPDG100-3;
1178	if ( pRPDG == R1) nRPDG=-R2;
1179	else nRPDG=-R1; // pRPDG == R2
1180	}
1181	break;
1182	}
1183	break;
1184	case 4: // ------------------------------------> QaQ + aDiQDiQ = QaQ (double)
1185	switch (af)
1186	{
1187	case 1: // ....................... pLPDG > 7 && pRPDG <-7 === LL/RR
1188	order= 1;
1189	if(-cLPDG == pL1) nLPDG= pL2;
1190	else nLPDG= pL1;
1191	if( cRPDG == pR1) nRPDG=-pR2;
1192	else nRPDG=-pR1;
1193	break;
1194	case 2: // ...................... pLPDG > 7 && pRPDG <-7 === LR/RL
1195	order=-1;
1196	if(-cRPDG == pL1) nLPDG= pL2;
1197	else nLPDG= pL1;
1198	if( cLPDG == pR1) nRPDG=-pR2;
1199	else nRPDG=-pR1;
1200	break;
1201	case 3: // ...................... pRPDG > 7 && pLPDG <-7 === LL/RR
1202	order= 1;
1203	if( cLPDG == pL1) nLPDG=-pL2;
1204	else nLPDG=-pL1;
1205	if(-cRPDG == pR1) nRPDG= pR2;
1206	else nRPDG= pR1;
1207	break;
1208	case 4: // ...................... pRPDG > 7 && pLPDG <-7 === LR/RL
1209	order=-1;
1210	if( cRPDG == pL1) nLPDG=-pL2;
1211	else nLPDG=-pL1;
1212	if(-cLPDG == pR1) nRPDG= pR2;
1213	else nRPDG= pR1;
1214	break;
1215	}
1216	break;
1217	case 5: // ------------------------------------> aDiQDiQ + QaQ = QaQ (double)
1218	switch (af)
1219	{
1220	case 1: // ...................... cLPDG > 7 && cRPDG <-7 === LL/RR
1221	order= 1;
1222	if(-pLPDG == L1) nLPDG= L2;
1223	else nLPDG= L1;
1224	if( pRPDG == R1) nRPDG=-R2;
1225	else nRPDG=-R1;
1226	break;
1227	case 2: // ...................... cLPDG > 7 && cRPDG <-7 === LR/RL
1228	order=-1;
1229	if(-pRPDG == L1) nRPDG= L2;
1230	else nRPDG= L1;
1231	if( pLPDG == R1) nLPDG=-R2;
1232	else nLPDG=-R1;
1233	break;
1234	case 3: // ...................... cRPDG > 7 && cLPDG <-7 === LL/RR
1235	order= 1;
1236	if( pLPDG == L1) nLPDG=-L2;
1237	else nLPDG=-L1;
1238	if(-pRPDG == R1) nRPDG= R2;
1239	else nRPDG= R1;
1240	break;
1241	case 4: // ...................... cRPDG > 7 && cLPDG <-7 === LR/RL
1242	order=-1;
1243	if( pRPDG == L1) nRPDG=-L2;
1244	else nRPDG=-L1;
1245	if(-pLPDG == R1) nLPDG= R2;
1246	else nLPDG= R1;
1247	break;
1248	}
1249	break;
1250	case 6: // ------------------------------------> QDiQ + aQaDiQ = QaQ (double)
1251	switch (af)
1252	{
1253	case 1:
1254	order=-1;
1255	if(-cRPDG == pL1) nLPDG= pL2;
1256	else nLPDG= pL1;
1257	if( pRPDG == L1) nRPDG= -L2;
1258	else nRPDG= -L1;
1259	break;
1260	case 2:
1261	order= 1;
1262	if(-cLPDG == pL1) nLPDG= pL2;
1263	else nLPDG= pL1;
1264	if( pRPDG == R1) nRPDG= -R2;
1265	else nRPDG= -R1;
1266	break;
1267	case 3:
1268	order= 1;
1269	if(-cRPDG == pR1) nRPDG= pR2;
1270	else nRPDG= pR1;
1271	if( pLPDG == L1) nLPDG= -L2;
1272	else nLPDG= -L1;
1273	break;
1274	case 4:
1275	order=-1;
1276	if(-cLPDG == pR1) nRPDG= pR2;
1277	else nRPDG= pR1;
1278	if( pLPDG == R1) nLPDG= -R2;
1279	else nLPDG= -R1;
1280	break;
1281	case 5:
1282	order=-1;
1283	if(-pRPDG == L1) nRPDG= L2;
1284	else nRPDG= L1;
1285	if( cRPDG == pL1) nLPDG=-pL2;
1286	else nLPDG=-pL1;
1287	break;
1288	case 6:
1289	order= 1;
1290	if(-pLPDG == L1) nLPDG= L2;
1291	else nLPDG= L1;
1292	if( cRPDG == pR1) nRPDG=-pR2;
1293	else nRPDG=-pR1;
1294	break;
1295	case 7:
1296	order= 1;
1297	if(-pRPDG == R1) nRPDG= R2;
1298	else nRPDG= R1;
1299	if( cLPDG == pL1) nLPDG=-pL2;
1300	else nLPDG=-pL1;
1301	break;
1302	case 8:
1303	order=-1;
1304	if(-pLPDG == R1) nLPDG= R2;
1305	else nLPDG= R1;
1306	if( cLPDG == pR1) nRPDG=-pR2;
1307	else nRPDG=-pR1;
1308	break;
1309	}
1310	break;
1311	}
1312	if(!order) G4cerr<<"-Warning-G4QInel::Constr: t="<<tf<<", a="<<af<<", cL="<<cLPDG
1313	<<", cR="<<cRPDG<<", pL="<<pLPDG<<", pR="<<pRPDG<<G4endl;
1314	else
1315	{
1316	// With theNewHypotheticalPartons the min mass must be calculated & compared
1317	G4int LT=1;
1318	if(std::abs(nLPDG) > 7) ++LT;
1319	G4int RT=1;
1320	if(std::abs(nRPDG) > 7) ++RT;
1321	G4double minM=0.;
1322	G4bool sing=true;
1323	if(cLT==2 && cRT==2)
1324	{
1325	G4int aLPDG=0;
1326	G4int aRPDG=0;
1327	if(cLPDG>0)
1328	{
1329	aLPDG=nLPDG/100;
1330	aRPDG=(-nRPDG)/100;
1331	}
1332	else //cRPDG>0
1333	{
1334	aRPDG=nRPDG/100;
1335	aLPDG=(-nLPDG)/100;
1336	}
1337	G4int nL1=aLPDG/10;
1338	G4int nL2=aLPDG%10;
1339	G4int nR1=aRPDG/10;
1340	G4int nR2=aRPDG%10;
1341	if(nL1!=nR1 && nL1!=nR2 && nL2!=nR1 && nL2!=nR2) // Unreducable DiQ-aDiQ
1342	{
1343	#ifdef debug
1344	G4cout<<"G4QIonIonCollis::Constr:aLPDG="<<aLPDG<<", aRPDG="<<aRPDG<<G4endl;
1345	#endif
1346	sing=false;
1347	G4QPDGCode tmp;
1348	std::pair<G4int,G4int> pB=tmp.MakeTwoBaryons(nL1, nL2, nR1, nR2);
1349	minM=G4QPDGCode(pB.first).GetMass()+G4QPDGCode(pB.second).GetMass();
1350	}
1351	}
1352	if(sing)
1353	{
1354	std::pair<G4int,G4int> newPair = std::make_pair(nLPDG,nRPDG);
1355	G4QContent newStQC(newPair); // NewString QuarkContent
1356	#ifdef debug
1357	G4cout<<"G4QIn::Con: LPDG="<<nLPDG<<",RPDG="<<nRPDG<<",QC="<<newStQC<<G4endl;
1358	#endif
1359	G4int minPDG=newStQC.GetSPDGCode(); // PDG of the Lightest Hadron=String
1360	minM=G4QPDGCode(minPDG).GetMass() + mPi0; // Min SingleHadron=String Mass
1361	}
1362	// Compare this mass
1363	G4bool win=false;
1364	G4double pSM=0.;
1365	if(pSM2>0.) pSM=std::sqrt(pSM2);
1366	if(minC && pSM2 > maxiM2) // Up to now any positive mass is good
1367	{
1368	maxiM2=pSM2;
1369	win=true;
1370	}
1371	else if(!minC \|\| pSM > minM)
1372	{
1373	pExcess=pSM-minM;
1374	if(minC \|\| pExcess > excess)
1375	{
1376	minC=false;
1377	excess=pExcess;
1378	win=true;
1379	}
1380	}
1381	if(win)
1382	{
1383	sst=pst;
1384	sLPDG=nLPDG;
1385	sRPDG=nRPDG;
1386	sOrd=order;
1387	}
1388	} // End of IF(new partons are created)
1389	} // End of IF(compatible partons)
1390	//} // End of positive squared mass of the fused string
1391	} // End of the LOOP over the possible partners (with the exclusive if for itself)
1392	if(sOrd) // The best pStringCandidate was found
1393	{
1394	G4LorentzVector cL4M=cLeft->Get4Momentum();
1395	G4LorentzVector cR4M=cRight->Get4Momentum();
1396	G4QParton* pLeft=(*sst)->GetLeftParton();
1397	G4QParton* pRight=(*sst)->GetRightParton();
1398	G4LorentzVector pL4M=pLeft->Get4Momentum();
1399	G4LorentzVector pR4M=pRight->Get4Momentum();
1400	#ifdef debug
1401	G4cout<<"G4QIonIonCollis::Constr:cS4M="<<cS4M<<" fused/w pS4M="<<pL4M+pR4M<<G4endl;
1402	#endif
1403	if(sOrd>0)
1404	{
1405	pL4M+=cL4M;
1406	pR4M+=cR4M;
1407	}
1408	else
1409	{
1410	pL4M+=cR4M;
1411	pR4M+=cL4M;
1412	}
1413	pLeft->SetPDGCode(sLPDG);
1414	pLeft->Set4Momentum(pL4M);
1415	pRight->SetPDGCode(sRPDG);
1416	pRight->Set4Momentum(pR4M);
1417	delete (*ist);
1418	strings.erase(ist);
1419	#ifdef debug
1420	G4LorentzVector ss4M=pL4M+pR4M;
1421	G4cout<<"G4QIonIonCollision::Constr:Created,4M="<<ss4M<<",m2="<<ss4M.m2()<<G4endl;
1422	#endif
1423	if( ist != strings.begin() ) // To avoid going below begin() (for Windows)
1424	{
1425	ist--;
1426	con=false;
1427	#ifdef debug
1428	G4cout<<"G4QIonIonCollision::Construct: * IST Decremented *"<<G4endl;
1429	#endif
1430	}
1431	else
1432	{
1433	con=true;
1434	#ifdef debug
1435	G4cout<<"G4QIonIonCollision::Construct: * IST Begin *"<<G4endl;
1436	#endif
1437	break;
1438	}
1439	} // End of the IF(the best partnerString candidate was found)
1440	else
1441	{
1442	#ifdef debug
1443	G4cout<<"-Warning-G4QInel::Const: S4M="<<cS4M<<",M2="<<cSM2<<" Leave ASIS"<<G4endl;
1444	#endif
1445	++problem;
1446	con=false;
1447	}
1448	}
1449	else con=false;
1450	} // End of loop over ist iterator
1451	#ifdef debug
1452	G4cout<<"G4QIonIonCollision::Construct: * IST While * , con="<<con<<G4endl;
1453	#endif
1454	} // End of "con" while
1455	#ifdef edebug
1456	// This print has meaning only if something appear between it and the StringFragmLOOP
1457	G4LorentzVector t4M=theProjNucleus.Get4Momentum()+theTargNucleus.Get4Momentum();//NucInLS
1458	G4int rC=totChg-theProjNucleus.GetZ()-theTargNucleus.GetZ();
1459	G4int rB=totBaN-theProjNucleus.GetA()-theTargNucleus.GetA();
1460	G4int nStr=strings.size();
1461	G4cout<<"-EMCLS-G4QIn::Const: AfterSUPPRESION #ofS="<<nStr<<",tNuc4M(E=M)="<<sum<<G4endl;
1462	for(G4int i=0; i<nStr; i++)
1463	{
1464	G4LorentzVector strI4M=strings[i]->Get4Momentum();
1465	t4M+=strI4M;
1466	G4int sChg=strings[i]->GetCharge();
1467	rC-=sChg;
1468	G4int sBaN=strings[i]->GetBaryonNumber();
1469	rB-=sBaN;
1470	G4cout<<"-EMCLS-G4QIonIonCollision::Construct: St#"<<i<<", 4M="<<strI4M<<", M="
1471	<<strI4M.m()<<", C="<<sChg<<", B="<<sBaN<<G4endl;
1472	}
1473	G4cout<<"-EMCLS-G4QInel::Construct: r4M="<<t4M-totLS4M<<", rC="<<rC<<", rB="<<rB<<G4endl;
1474	#endif
1475	//
1476	// --- If a problem is foreseen then the DiQaDiQ strings should be reduced if possible --
1477	//
1478	#ifdef debug
1479	G4cout<<"G4QIonIonCollision::Construct: problem="<<problem<<G4endl;
1480	#endif
1481	if(problem)
1482	{
1483	G4int nOfStr=strings.size();
1484	#ifdef debug
1485	G4cout<<"G4QIonIonCollision::Constr:SecurityDiQaDiQReduction, #OfStr="<<nOfStr<<G4endl;
1486	#endif
1487	for (G4int astring=0; astring < nOfStr; astring++)
1488	{
1489	G4QString* curString=strings[astring];
1490	G4QParton* cLeft=curString->GetLeftParton();
1491	G4QParton* cRight=curString->GetRightParton();
1492	G4int LT=cLeft->GetType();
1493	G4int RT=cRight->GetType();
1494	G4int sPDG=cLeft->GetPDGCode();
1495	G4int nPDG=cRight->GetPDGCode();
1496	if(LT==2 && RT==2)
1497	{
1498	#ifdef debug
1499	G4cout<<"G4QIonIonCollis::Constr:TrySelfReduString,L="<<sPDG<<",R="<<nPDG<<G4endl;
1500	#endif
1501	if( cLeft->ReduceDiQADiQ(cLeft, cRight) ) // DiQ-aDiQ pair was successfully reduced
1502	{
1503	sPDG=cLeft->GetPDGCode();
1504	nPDG=cRight->GetPDGCode();
1505	#ifdef debug
1506	G4cout<<"+G4QInel::Const:#"<<astring<<" Reduced, L="<<sPDG<<", R="<<nPDG<<G4endl;
1507	#endif
1508	}
1509	#ifdef debug
1510	else G4cout<<"--*--G4QInel::Const: #"<<astring<<" DQ-aDQ reduction Failed"<<G4endl;
1511	#endif
1512	} // End of the found DiQ/aDiQ pair
1513	else if(sPDG==3 && nPDG==-3)
1514	{
1515	sPDG= 1;
1516	nPDG=-1;
1517	cLeft->SetPDGCode(sPDG);
1518	cRight->SetPDGCode(nPDG);
1519	}
1520	else if(sPDG==-3 && nPDG==3)
1521	{
1522	sPDG=-1;
1523	nPDG= 1;
1524	cLeft->SetPDGCode(sPDG);
1525	cRight->SetPDGCode(nPDG);
1526	}
1527	}
1528	SwapPartons();
1529	} // End of IF(problem)
1530	#ifdef edebug
1531	G4LorentzVector u4M=theProjNucleus.Get4Momentum()+theTargNucleus.Get4Momentum();//NucInLS
1532	G4int rCh=totChg-theProjNucleus.GetZ()-theTargNucleus.GetZ();
1533	G4int rBa=totBaN-theProjNucleus.GetA()-theTargNucleus.GetA();
1534	G4int nStri=strings.size();
1535	G4cout<<"-EMCLS-G4QIn::Const: FinalConstruct, #ofSt="<<nStri<<",tN4M(E=M)="<<t4M<<G4endl;
1536	for(G4int i=0; i<nStri; i++)
1537	{
1538	G4LorentzVector strI4M=strings[i]->Get4Momentum();
1539	u4M+=strI4M;
1540	G4int sChg=strings[i]->GetCharge();
1541	rCh-=sChg;
1542	G4int sBaN=strings[i]->GetBaryonNumber();
1543	rBa-=sBaN;
1544	G4cout<<"-EMCLS-G4QIonIonCollision::Construct: St#"<<i<<", 4M="<<strI4M<<", M="
1545	<<strI4M.m()<<", C="<<sChg<<", B="<<sBaN<<G4endl;
1546	}
1547	G4cout<<"-EMCLS-G4QInel::Construct: r4M="<<u4M-totLS4M<<",rC="<<rCh<<",rB="<<rBa<<G4endl;
1548	#endif
1549	}
1550
1551	G4QIonIonCollision::~G4QIonIonCollision()
1552	{
1553	std::for_each(strings.begin(), strings.end(), DeleteQString() );
1554	}
1555
1556	G4QHadronVector* G4QIonIonCollision::Fragment()
1557	{ // This is the member function fragmenting Strings & Quasmons (in nuclear matter)
1558	#ifdef debug
1559	G4cout<<"*****>G4QIonIonCollision::Fragment: Called**, Res="<<theResult<<G4endl;
1560	#endif
1561	G4int striNum=strings.size(); // Find out if there're strings
1562	G4int hadrNum=theResult->size(); // Find out if there're hadrons
1563	#ifdef edebug
1564	G4int nQm=theQuasmons.size();
1565	G4LorentzVector totLS4M=theProjNucleus.Get4Momentum()+theTargNucleus.Get4Momentum(); //LS
1566	G4int totChg=theProjNucleus.GetZ()+theTargNucleus.GetZ();
1567	G4int totBaN=theTargNucleus.GetA()+theTargNucleus.GetA();
1568	G4cout<<"-EMCLS-G4QIn::Fragment: CHECKRecovery, #ofS="<<striNum<<", Nuc4M(E=M)="<<totLS4M
1569	<<",#Q="<<nQm<<",#H="<<hadrNum<<G4endl;
1570	for(G4int i=0; i < striNum; i++)
1571	{
1572	G4LorentzVector strI4M=strings[i]->Get4Momentum();
1573	totLS4M+=strI4M;
1574	G4int sChg=strings[i]->GetCharge();
1575	totChg+=sChg;
1576	G4int sBaN=strings[i]->GetBaryonNumber();
1577	totBaN+=sBaN;
1578	G4cout<<"-EMCLS-G4QIonIonCollision::Fragment:String#"<<i<<", 4M="<<strI4M<<", M="
1579	<<strI4M.m()<<", C="<<sChg<<", B="<<sBaN<<G4endl;
1580	}
1581	for(G4int i=0; i < nQm; i++)
1582	{
1583	G4LorentzVector hI4M=theQuasmons[i]->Get4Momentum();
1584	totLS4M+=hI4M;
1585	G4int hChg=theQuasmons[i]->GetCharge();
1586	totChg+=hChg;
1587	G4int hBaN=theQuasmons[i]->GetBaryonNumber();
1588	totBaN+=hBaN;
1589	G4cout<<"-EMCLS-G4QIonIonCollision::Fragment: Quasmon#"<<i<<", 4M="<<hI4M<<", C="<<hChg
1590	<<", B="<<hBaN<<G4endl;
1591	}
1592	for(G4int i=0; i < hadrNum; i++)
1593	{
1594	G4LorentzVector hI4M=(*theResult)[i]->Get4Momentum();
1595	totLS4M+=hI4M;
1596	G4int hChg=(*theResult)[i]->GetCharge();
1597	totChg+=hChg;
1598	G4int hBaN=(*theResult)[i]->GetBaryonNumber();
1599	totBaN+=hBaN;
1600	G4cout<<"-EMCLS-G4QIn::Fragment:H#"<<i<<",4M="<<hI4M<<",C="<<hChg<<",B="<<hBaN<<G4endl;
1601	}
1602	#endif
1603	#ifdef debug
1604	G4cout<<"***>G4QIonIonCollision::Fragm: #OfStr="<<striNum<<", #OfRes="<<hadrNum<<G4endl;
1605	#endif
1606	if(!striNum && hadrNum) // Quasi-elastic or decoupled p
1607	{
1608	#ifdef debug
1609	G4cout<<"*>G4QIonIonCollision::Fragm:Quasi-Elastic**, #OfResult="<<hadrNum<<G4endl;
1610	#endif
1611	return theResult;
1612	}
1613	else if(striNum) Breeder(); // Strings fragmentation
1614	else // No strings, make HadrNucleus
1615	{
1616	if(hadrNum)
1617	{
1618	for(G4int ih=0; ih<hadrNum; ih++) delete (*theResult)[ih];
1619	theResult->clear();
1620	}
1621	G4LorentzVector rp4M=theProjNucleus.Get4Momentum(); // Nucleus 4-momentum in LS
1622	G4int rpPDG=theProjNucleus.GetPDG(); // Nuclear PDG
1623	G4QHadron* resPNuc = new G4QHadron(rpPDG,rp4M); // Nucleus -> Hadron
1624	theResult->push_back(resPNuc); // Fill the residual nucleus
1625	G4LorentzVector rt4M=theTargNucleus.Get4Momentum(); // Nucleus 4-momentum in LS
1626	G4int rtPDG=theTargNucleus.GetPDG(); // Nuclear PDG
1627	G4QHadron* resTNuc = new G4QHadron(rtPDG,rt4M); // Nucleus -> Hadron
1628	theResult->push_back(resTNuc); // Fill the residual nucleus
1629	}
1630	G4int nQuas=theQuasmons.size(); // Size of the Quasmon OUTPUT
1631	G4int theRS=theResult->size(); // Size of Hadron Output by now
1632	#ifdef debug
1633	G4cout<<"***>G4QIonIonCollision::Fragm:beforeEnv, #OfQ="<<nQuas<<",#OfR="<<theRS<<G4endl;
1634	#endif
1635	if(nQuas && theRS)
1636	{
1637
1638	G4QHadron* resNuc = (*theResult)[theRS-1]; // Pointer to Residual Nucleus
1639	G4LorentzVector resNuc4M = resNuc->Get4Momentum(); // 4-Momentum of the Nucleuz
1640	G4int resNucPDG= resNuc->GetPDGCode(); // PDG Code of the Nucleus
1641	G4QNucleus theEnv(resNucPDG); // NucleusHadron->NucleusAtRest
1642	delete resNuc; // Delete resNucleus as aHadron
1643	theResult->pop_back(); // Exclude the nucleus from HV
1644	--theRS; // Reduce the OUTPUT by theNucl
1645	#ifdef pdebug
1646	G4cout<<"G4QIonIonCollision::Fragm:#OfRemainingHadron="<<theRS<<", A="<<theEnv<<G4endl;
1647	#endif
1648	// Now we need to be sure that the compound nucleus is heavier than the Ground State
1649	for(G4int j=theRS-1; j>-2; --j) // try to reach M_compound>M_GS
1650	{
1651	G4LorentzVector qsum4M=resNuc4M; // Proto compound 4-momentum
1652	G4QContent qsumQC=theEnv.GetQCZNS(); // Proto compound Quark Content
1653	#ifdef pdebug
1654	G4cout<<"G4QIonIonCollis::Fragm: rN4M"<<qsum4M<<qsum4M.m()<<",rNQC="<<qsumQC<<G4endl;
1655	#endif
1656	G4Quasmon* firstQ=0; // Prototype of theFirstQuasmon
1657	G4LorentzVector first4M; // Proto of the FirstQuasmon 4M
1658	G4QContent firstQC; // Proto of the FirstQuasmon QC
1659	for(G4int i=0; i<nQuas; ++i) // LOOP over Quasmons
1660	{
1661	G4Quasmon* curQuasm=theQuasmons[i]; // current Quasmon
1662	G4LorentzVector cur4M=curQuasm->Get4Momentum(); // 4-Mom of the Quasmon
1663	G4QContent curQC=curQuasm->GetQC(); // Quark Content of the Quasmon
1664	qsum4M+=cur4M; // Add quasmon's 4-momentum
1665	qsumQC+=curQC; // Add quasmon's Quark Content
1666	#ifdef pdebug
1667	G4cout<<"G4QIonIonCollis::Fragm: Q#"<<i<<", QQC="<<curQC<<", sQC="<<qsumQC<<G4endl;
1668	#endif
1669	if(!i) // Remember 1-st for correction
1670	{
1671	firstQ =curQuasm;
1672	first4M=cur4M;
1673	firstQC=curQC;
1674	}
1675	}
1676	G4int miPDG=qsumQC.GetSPDGCode(); // PDG of minM of hadron/fragm.
1677	G4double gsM=0.; // Proto minM of had/frag forQC
1678	if(miPDG == 10)
1679	{
1680	G4QChipolino QCh(qsumQC); // define TotNuc as a Chipolino
1681	gsM=QCh.GetQPDG1().GetMass()+QCh.GetQPDG2().GetMass(); // Sum of Hadron Masses
1682	//gsM=theWorld->GetQParticle(QCh.GetQPDG1())->MinMassOfFragm() +
1683	// theWorld->GetQParticle(QCh.GetQPDG2())->MinMassOfFragm();
1684	}
1685	// @@ it is not clear, why it does not work ?!
1686	//else if(miPDG>80000000) // Compound Nucleus
1687	//{
1688	// G4QNucleus rtN(qsumQC); // CreatePseudoNucl for totComp
1689	// gsM=rtN.GetGSMass(); // MinMass of residQ+(Env-ParC)
1690	//}
1691	else if(miPDG < 80000000 && std::abs(miPDG)%10 > 2)
1692	gsM=theWorld->GetQParticle(G4QPDGCode(miPDG))->MinMassOfFragm();
1693	else gsM=G4QPDGCode(miPDG).GetMass(); // minM of hadron/fragm. for QC
1694	G4double reM=qsum4M.m(); // real mass of the compound
1695	#ifdef pdebug
1696	G4cout<<"G4QIonIonCollision::Fragm: PDG="<<miPDG<<", rM="<<reM<<",GSM="<<gsM<<G4endl;
1697	#endif
1698	if(reM > gsM) break; // CHIPS can be called
1699	if(j > -1) // Can try to add hadrons to Q0
1700	{
1701	G4QHadron* cH = (*theResult)[j]; // Pointer to the last Hadron
1702	G4LorentzVector h4M = cH->Get4Momentum(); // 4-Momentum of the Hadron
1703	G4QContent hQC = cH->GetQC(); // QC of the Hadron
1704	firstQ->Set4Momentum(first4M+h4M); // Update the Quasmon's 4-Mom
1705	firstQ->SetQC(firstQC+hQC); // Update the Quasmon's QCont
1706	delete cH; // Delete the Hadron
1707	theResult->pop_back(); // Exclude the hadron from HV
1708	#ifdef pdebug
1709	G4cout<<"G4QInel::Fragm: H#"<<j<<", hQC="<<hQC<<",hPDG="<<cH->GetPDGCode()<<G4endl;
1710	#endif
1711	}
1712	else
1713	{
1714	G4cerr<<"***G4QIonIonCollision::Fra:PDG="<<miPDG<<",M="<<reM<<",GSM="<<gsM<<G4endl;
1715	G4Exception("G4QIonIonCollision::Fragment:","27",FatalException,"Can'tRecoverGSM");
1716	}
1717	}
1718	G4double nucE=resNuc4M.e(); // Total energy of the nuclEnv
1719	if(nucE<1.E-12) nucE=0.; // Computer accuracy safety
1720	G4ThreeVector nucVel(0.,0.,0.); // Proto of the NucleusVelocity
1721	G4QHadronVector* output=0; // NucleusFragmentation Hadrons
1722	G4QEnvironment* pan= new G4QEnvironment(theEnv); // ---> DELETED --->----------+
1723	#ifdef pdebug
1724	G4cout<<"G4QIonIonCollision::Fragm: nucE="<<nucE<<", nQ="<<nQuas<<G4endl; // \|
1725	#endif
1726	if(nucE) nucVel=resNuc4M.vect()/nucE; // The NucleusVelocity \|
1727	for(G4int i=0; i<nQuas; ++i) // LOOP over Quasmons \|
1728	{ // \|
1729	G4Quasmon* curQuasm=theQuasmons[i]; // current Quasmon \|
1730	if(nucE) curQuasm->Boost(-nucVel); // Boost it to CMS of Nucleus \|
1731	pan->AddQuasmon(curQuasm); // Fill the predefined Quasmon\|
1732	#ifdef pdebug
1733	G4LorentzVector cQ4M=curQuasm->Get4Momentum(); // Just for printing \|
1734	G4cout<<"G4QIonIonCollision::Fragment: Quasmon# "<<i<<" added, 4M="<<cQ4M<<G4endl;//\|
1735	#endif
1736	} // \|
1737	try // \|
1738	{ // \|
1739	delete output; // \|
1740	output = pan->Fragment();// DESTROYED after theHadrons are transferred to theResult \|
1741	} // \| \|
1742	catch (G4QException& error) // \| \|
1743	{ // \| \|
1744	G4cerr<<"***G4QIonIonCollision::Fragment: G4QE Exception is catched"<<G4endl; // \| \|
1745	G4Exception("G4QIonIonCollision::Fragment:","27",FatalException,"CHIPSCrash");// \| \|
1746	} // \| \|
1747	delete pan; // Delete the Nuclear Environment <-----<--+-+
1748	if(output) // Output exists \|
1749	{ // \|
1750	G4int nOut=output->size(); // #ofHadrons in the Nuclear Fragmentation \|
1751	for(G4int j=0; j<nOut; j++) // LOOP over Hadrons transferring to LS \|
1752	{ // \|
1753	G4QHadron* curHadron=(*output)[j]; // Hadron from the nucleus fragmentation \|
1754	if(nucE) curHadron->Boost(nucVel); // Boost it back to Laboratory System \|
1755	theResult->push_back(curHadron); // Transfer it to the result \|
1756	} // \|
1757	delete output; // Delete the OUTPUT <-----<-----<-----<---+
1758	}
1759	}
1760	else if(!striNum) G4cout<<"-Warning-G4QIonIonCollision::Fragment:NothingWasDone"<<G4endl;
1761	#ifdef debug
1762	G4cout<<"====>G4QIonIonCollision::Fragment: Final #OfResult="<<theResult->size()<<G4endl;
1763	#endif
1764	G4int nQ =theQuasmons.size();
1765	if(nQ) theQuasmons.clear(); // @@ Not necesary ?
1766	#ifdef edebug
1767	G4LorentzVector f4M(0.,0.,0.,0.); // Sum of the Result in LS
1768	G4int fCh=totChg;
1769	G4int fBN=totBaN;
1770	G4int nHd=theResult->size();
1771	G4cout<<"-EMCLS-G4QIonIonCollision::Fragment: #ofHadr="<<nHd<<",#OfQuasm="<<nQ<<G4endl;
1772	for(G4int i=0; i<nHd; i++)
1773	{
1774	G4LorentzVector hI4M=(*theResult)[i]->Get4Momentum();
1775	f4M+=hI4M;
1776	G4int hChg=(*theResult)[i]->GetCharge();
1777	fCh-=hChg;
1778	G4int hBaN=(*theResult)[i]->GetBaryonNumber();
1779	fBN-=hBaN;
1780	G4cout<<"-EMCLS-G4QIonIonCollision::Fragment: Hadron#"<<i<<", 4M="<<hI4M<<", PDG="
1781	<<(*theResult)[i]->GetPDGCode()<<", C="<<hChg<<", B="<<hBaN<<G4endl;
1782	}
1783	G4cout<<"-EMCLS-G4QInel::Fragm:, r4M="<<f4M-totLS4M<<", rC="<<fCh<<",rB="<<fBN<<G4endl;
1784	#endif
1785	return theResult;
1786	} // End of fragmentation
1787
1788	void G4QIonIonCollision::Breeder()
1789	{ // This is the member function, which returns the resulting vector of Hadrons & Quasmons
1790	static const G4double eps = 0.001; // Tolerance in MeV
1791	//
1792	// ------------ At this point the strings are fragmenting to hadrons in LS -------------
1793	//
1794	#ifdef edebug
1795	G4LorentzVector totLS4M=theProjNucleus.Get4Momentum()+theTargNucleus.Get4Momentum(); //LS
1796	G4int totChg=theProjNucleus.GetZ()+theTargNucleus.GetZ();
1797	G4int totBaN=theProjNucleus.GetA()+theTargNucleus.GetA();
1798	G4int nStri=strings.size();
1799	G4cout<<"-EMCLS-G4QIn::Breed: CHECKRecovery #ofS="<<nStri<<",N4M(E=M)="<<totLS4M<<G4endl;
1800	for(G4int i=0; i<nStri; i++)
1801	{
1802	G4LorentzVector strI4M=strings[i]->Get4Momentum();
1803	totLS4M+=strI4M;
1804	G4int sChg=strings[i]->GetCharge();
1805	totChg+=sChg;
1806	G4int sBaN=strings[i]->GetBaryonNumber();
1807	totBaN+=sBaN;
1808	G4cout<<"-EMCLS-G4QIonIonCollision::Breeder: St#"<<i<<", 4M="<<strI4M<<", M="
1809	<<strI4M.m()<<", C="<<sChg<<", B="<<sBaN<<G4endl;
1810	}
1811	#endif
1812	G4int nOfStr=strings.size();
1813	#ifdef debug
1814	G4cout<<"G4QIonIonCollision::Breeder: BeforeFragmentation, #OfStr="<<nOfStr<<G4endl;
1815	#endif
1816	G4LorentzVector ft4M(0.,0.,0.,0.);
1817	G4QContent ftQC(0,0,0,0,0,0);
1818	G4bool ftBad=false;
1819	for(G4int i=0; i < nOfStr; ++i)
1820	{
1821	G4LorentzVector pS4M=strings[i]->Get4Momentum(); // String 4-momentum
1822	ft4M+=pS4M;
1823	G4QContent pSQC=strings[i]->GetQC(); // String Quark Content
1824	ftQC+=pSQC;
1825	if(pS4M.m2() < 0.) ftBad=true;
1826	#ifdef debug
1827	G4cout<<">G4QIonIonCollision::Breed:1stTest,S#"<<i<<",4M="<<pS4M<<",QC="<<pSQC<<G4endl;
1828	#endif
1829	}
1830	if(ftBad)
1831	{
1832	G4Quasmon* stringQuasmon = new G4Quasmon(ftQC, ft4M);
1833	#ifdef debug
1834	G4cout<<"->G4QIonIonCollision::Breed:TotQ,QC="<<ftQC<<",4M="<<ft4M<<ft4M.m()<<G4endl;
1835	#endif
1836	theQuasmons.push_back(stringQuasmon);
1837	G4LorentzVector rp4M=theProjNucleus.Get4Momentum(); // Nucleus 4-momentum in LS
1838	G4int rpPDG=theProjNucleus.GetPDG();
1839	G4QHadron* resPNuc = new G4QHadron(rpPDG,rp4M);
1840	theResult->push_back(resPNuc); // Fill the residual projectile nucleus
1841	G4LorentzVector rt4M=theTargNucleus.Get4Momentum(); // Nucleus 4-momentum in LS
1842	G4int rtPDG=theTargNucleus.GetPDG();
1843	G4QHadron* resTNuc = new G4QHadron(rtPDG,rt4M);
1844	theResult->push_back(resTNuc); // Fill the residual target nucleus
1845	return;
1846	}
1847	for (G4int astring=0; astring < nOfStr; astring++)
1848	{
1849	#ifdef edebug
1850	G4LorentzVector sum=theProjNucleus.Get4Momentum()+theTargNucleus.Get4Momentum(); //inLS
1851	G4int rChg=totChg-theProjNucleus.GetZ()-theTargNucleus.GetZ();
1852	G4int rBaN=totBaN-theProjNucleus.GetA()-theTargNucleus.GetA();
1853	G4int nOfHadr=theResult->size();
1854	G4cout<<"-EMCLS-G4QIonIonCollision::Breed:#ofSt="<<nOfStr<<",#ofHad="<<nOfHadr<<G4endl;
1855	for(G4int i=astring; i<nOfStr; i++)
1856	{
1857	G4LorentzVector strI4M=strings[i]->Get4Momentum();
1858	sum+=strI4M;
1859	G4int sChg=strings[i]->GetCharge();
1860	rChg-=sChg;
1861	G4int sBaN=strings[i]->GetBaryonNumber();
1862	rBaN-=sBaN;
1863	G4cout<<"-EMCLS-G4QI::Breed:S#"<<i<<",4M="<<strI4M<<",C="<<sChg<<",B="<<sBaN<<G4endl;
1864	}
1865	for(G4int i=0; i<nOfHadr; i++)
1866	{
1867	G4LorentzVector hI4M=(*theResult)[i]->Get4Momentum();
1868	sum+=hI4M;
1869	G4int hChg=(*theResult)[i]->GetCharge();
1870	rChg-=hChg;
1871	G4int hBaN=(*theResult)[i]->GetBaryonNumber();
1872	rBaN-=hBaN;
1873	G4cout<<"-EMCLS-G4QIn::Breed: H#"<<i<<",4M="<<hI4M<<",C="<<hChg<<",B="<<hBaN<<G4endl;
1874	}
1875	G4cout<<"....-EMCLS-G4QInel::Br:r4M="<<sum-totLS4M<<",rC="<<rChg<<",rB="<<rBaN<<G4endl;
1876	#endif
1877	G4QString* curString=strings[astring];
1878	if(!curString->GetDirection()) continue; // Historic for the dead strings: DoesNotWork
1879	#ifdef edebug
1880	G4int curStrChg = curString->GetCharge();
1881	G4int curStrBaN = curString->GetBaryonNumber();
1882	#endif
1883	G4LorentzVector curString4M = curString->Get4Momentum();
1884	#ifdef debug
1885	G4cout<<"====>G4QIonIonCollision::Breeder: String#"<<astring<<",s4M/m="<<curString4M
1886	<<curString4M.m()<<", LPart="<<curString->GetLeftParton()->GetPDGCode()
1887	<<", RPart="<<curString->GetRightParton()->GetPDGCode()<<G4endl;
1888	#endif
1889	G4QHadronVector* theHadrons = 0; // Prototype of theStringFragmentationOUTPUT
1890	theHadrons=curString->FragmentString(true);// !! Fragmenting the String !!
1891	if (!theHadrons) // The string can not be fragmented
1892	{
1893	// First try to correct the diQ-antiDiQ strings, converting them to Q-antiQ
1894	G4QParton* cLeft=curString->GetLeftParton();
1895	G4QParton* cRight=curString->GetRightParton();
1896	G4int sPDG=cLeft->GetPDGCode();
1897	G4int nPDG=cRight->GetPDGCode();
1898	G4int LT=cLeft->GetType();
1899	G4int RT=cRight->GetType();
1900	G4int LS=LT+RT;
1901	if(LT==2 && RT==2)
1902	{
1903	#ifdef debug
1904	G4cout<<"G4QIonIonCollision::Breed:TryReduceString, L="<<sPDG<<",R="<<nPDG<<G4endl;
1905	#endif
1906	if( cLeft->ReduceDiQADiQ(cLeft, cRight) ) // DiQ-aDiQ pair was successfully reduced
1907	{
1908	LT=1;
1909	RT=1;
1910	LS=2;
1911	sPDG=cLeft->GetPDGCode();
1912	nPDG=cRight->GetPDGCode();
1913	#ifdef debug
1914	G4cout<<"G4QIonIonCollision::Breed:AfterReduction,L="<<sPDG<<",R="<<nPDG<<G4endl;
1915	#endif
1916	theHadrons=curString->FragmentString(true);//!! Try to fragment the new String !!
1917	cLeft=curString->GetLeftParton();
1918	cRight=curString->GetRightParton();
1919	#ifdef debug
1920	G4cout<<"G4QInel::Breed:L="<<cLeft->Get4Momentum()<<",R="<<cRight->Get4Momentum()
1921	<<G4endl;
1922	#endif
1923	}
1924	#ifdef debug
1925	else G4cout<<"^G4QIonIonCollision::Breed: DQ-aDQ reduction to Q-aQ Failed"<<G4endl;
1926	#endif
1927	} // End of the SelfReduction
1928	#ifdef debug
1929	G4cout<<"G4QIonIonCollision::Breed:AfterRedAttempt, theH="<<theHadrons<<", L4M="
1930	<<cLeft->Get4Momentum()<<", R4M="<<cRight->Get4Momentum()<<G4endl;
1931	#endif
1932	unsigned next=astring+1; // The next string position
1933	if (!theHadrons) // The string can not be fragmented
1934	{
1935	G4int fusionDONE=0; // StringFusion didn't happen (1=Fuse L+L/R+R, -1=Fuse L+R/R+L)
1936	if(next < strings.size()) // TheString isn't theLastString can fuse
1937	{
1938	G4int fustr=0; // The found partner index (never can be 0)
1939	G4int swap=0; // working interger for swapping parton PDG
1940	G4double Vmin=DBL_MAX; // Prototype of the found Velocity Distance
1941	G4int dPDG=nPDG;
1942	G4int qPDG=sPDG;
1943	if(dPDG<-99 \|\| (dPDG>0&&dPDG<7) \|\| qPDG>99 \|\| (qPDG<0 && qPDG>-7))
1944	{
1945	swap=qPDG;
1946	qPDG=dPDG;
1947	dPDG=swap;
1948	}
1949	if(dPDG>99) dPDG/=100;
1950	if(qPDG<-99) qPDG=-(-qPDG)/100;
1951	#ifdef debug
1952	G4cout<<"G4QIonIonCollision::Breed:TryFuseStringS, q="<<qPDG<<", a="<<dPDG
1953	<<", n="<<next<<G4endl;
1954	#endif
1955	G4ThreeVector curV=curString4M.vect()/curString4M.e();
1956	G4int reduce=0; // a#of reduced Q-aQ pairs
1957	G4int restr=0; // To use beyon the LOOP for printing
1958	G4int MPS=0; // PLS for the selected string
1959	for (restr=next; restr < nOfStr; restr++)
1960	{
1961	reduce=0;
1962	G4QString* reString=strings[restr];
1963	G4QParton* Left=reString->GetLeftParton();
1964	G4QParton* Right=reString->GetRightParton();
1965	G4int uPDG=Left->GetPDGCode();
1966	G4int mPDG=Right->GetPDGCode();
1967	G4int PLT =Left->GetType();
1968	G4int PRT =Right->GetType();
1969	G4int aPDG=mPDG;
1970	G4int rPDG=uPDG;
1971	if(aPDG<-99 \|\| (aPDG>0 && aPDG<7) \|\| rPDG>99 \|\| (rPDG<0 && rPDG>-7))
1972	{
1973	swap=rPDG;
1974	rPDG=aPDG;
1975	aPDG=swap;
1976	}
1977	if(aPDG > 99) aPDG/=100;
1978	if(rPDG <-99) rPDG=-(-rPDG)/100;
1979	// Try to reduce two DQ-aDQ strings
1980	#ifdef debug
1981	G4cout<<"G4Qnel::Breed: TryReduce #"<<restr<<", q="<<rPDG<<",a="<<aPDG<<G4endl;
1982	#endif
1983	if(LT==2 && RT==2 && PLT==2 && PRT==2) // Have a chance for the reduction
1984	{
1985	G4int cQ1=(-qPDG)/10;
1986	G4int cQ2=(-qPDG)%10;
1987	G4int cA1=dPDG/10;
1988	G4int cA2=dPDG%10;
1989	G4int pQ1=(-rPDG)/10;
1990	G4int pQ2=(-rPDG)%10;
1991	G4int pA1=aPDG/10;
1992	G4int pA2=aPDG%10;
1993	#ifdef debug
1994	G4cout<<"G4QIonIonCollision::Breeder: cQ="<<cQ1<<","<<cQ2<<", cA="<<cA1<<","
1995	<<cA2<<", pQ="<<pQ1<<","<<pQ2<<", pA="<<pA1<<","<<pA2<<G4endl;
1996	#endif
1997	G4bool iQA = (cA1==pQ1 \|\| cA1==pQ2 \|\| cA2==pQ1 \|\| cA2==pQ2);
1998	G4bool iAQ = (cQ1==pA1 \|\| cQ1==pA2 \|\| cQ2==pA1 \|\| cQ2==pA2);
1999	if(iQA) reduce++;
2000	if(iAQ) reduce++;
2001	if (reduce==2) // Two quark pairs can be reduced
2002	{
2003	if(sPDG>0 && uPDG<0) // LL/RR Reduction
2004	{
2005	std::pair<G4int,G4int> resLL=ReducePair(sPDG/100, (-uPDG)/100);
2006	G4int newCL=resLL.first;
2007	G4int newPL=resLL.second;
2008	if(!newCL \|\| !newPL)
2009	{
2010	G4cerr<<"*G4QIonIonCollision::Breed:CL="<<newCL<<",PL="<<newPL<<G4endl;
2011	G4Exception("G4QIonIonCollision::Breeder:","72",FatalException,"2LL-");
2012	}
2013	std::pair<G4int,G4int> resRR=ReducePair((-nPDG)/100, mPDG/100);
2014	G4int newCR=resRR.first;
2015	G4int newPR=resRR.second;
2016	if(!newCR \|\| !newPR)
2017	{
2018	G4cerr<<"*G4QIonIonCollision::Breed:CR="<<newCR<<",PR="<<newPR<<G4endl;
2019	G4Exception("G4QIonIonCollision::Breeder:","72",FatalException,"2RR-");
2020	}
2021	cLeft->SetPDGCode(newCL); // Reset the left quark of curString
2022	cRight->SetPDGCode(-newCR); // Reset the right quark of curString
2023	Left->SetPDGCode(-newPL); // Reset the left quark of reString
2024	Right->SetPDGCode(newPR); // Reset the right quark of reString
2025	break; // Break out of the reString internal LOOP
2026	}
2027	else if(sPDG<0 && uPDG>0) // LL/RR Reduction
2028	{
2029	std::pair<G4int,G4int> resLL=ReducePair((-sPDG)/100, uPDG/100);
2030	G4int newCL=resLL.first;
2031	G4int newPL=resLL.second;
2032	if(!newCL \|\| !newPL)
2033	{
2034	G4cerr<<"*G4QIonIonCollision::Breed:CL="<<newCL<<",PL="<<newPL<<G4endl;
2035	G4Exception("G4QIonIonCollision::Breeder:","72",FatalException,"2LL+");
2036	}
2037	std::pair<G4int,G4int> resRR=ReducePair(nPDG/100, (-mPDG)/100);
2038	G4int newCR=resRR.first;
2039	G4int newPR=resRR.second;
2040	if(!newCR \|\| !newPR)
2041	{
2042	G4cerr<<"*G4QIonIonCollision::Breed:CR="<<newCR<<",PR="<<newPR<<G4endl;
2043	G4Exception("G4QIonIonCollision::Breeder:","72",FatalException,"2RR+");
2044	}
2045	cLeft->SetPDGCode(-newCL); // Reset the left quark of curString
2046	cRight->SetPDGCode(newCR); // Reset the right quark of curString
2047	Left->SetPDGCode(newPL); // Reset the left quark of reString
2048	Right->SetPDGCode(-newPR); // Reset the right quark of reString
2049	break; // Break out of the reString internal LOOP
2050	}
2051	else if(sPDG>0 && mPDG<0) // LR Reduction
2052	{
2053	std::pair<G4int,G4int> resLL=ReducePair(sPDG/100, (-mPDG)/100);
2054	G4int newCL=resLL.first;
2055	G4int newPR=resLL.second;
2056	if(!newCL \|\| !newPR)
2057	{
2058	G4cerr<<"*G4QIonIonCollision::Breed:CL="<<newCL<<",PR="<<newPR<<G4endl;
2059	G4Exception("G4QIonIonCollision::Breeder:","72",FatalException,"2-LR");
2060	}
2061	std::pair<G4int,G4int> resRR=ReducePair((-nPDG)/100, uPDG/100);
2062	G4int newCR=resRR.first;
2063	G4int newPL=resRR.second;
2064	if(!newCR \|\| !newPR)
2065	{
2066	G4cerr<<"*G4QIonIonCollision::Breed:CR="<<newCR<<",PL="<<newPL<<G4endl;
2067	G4Exception("G4QIonIonCollision::Breeder:","72",FatalException,"2-LR");
2068	}
2069	cLeft->SetPDGCode(newCL); // Reset the left quark of curString
2070	cRight->SetPDGCode(-newCR); // Reset the right quark of curString
2071	Left->SetPDGCode(newPL); // Reset the left quark of reString
2072	Right->SetPDGCode(-newPR); // Reset the right quark of reString
2073	break; // Break out of the reString internal LOOP
2074	}
2075	else // RL Reduction
2076	{
2077	std::pair<G4int,G4int> resLL=ReducePair((-sPDG)/100, mPDG/100);
2078	G4int newCL=resLL.first;
2079	G4int newPR=resLL.second;
2080	if(!newCL \|\| !newPR)
2081	{
2082	G4cerr<<"*G4QIonIonCollision::Breed:CL="<<newCL<<",PR="<<newPR<<G4endl;
2083	G4Exception("G4QIonIonCollision::Breeder:","72",FatalException,"2-RL");
2084	}
2085	std::pair<G4int,G4int> resRR=ReducePair(nPDG/100, (-uPDG)/100);
2086	G4int newCR=resRR.first;
2087	G4int newPL=resRR.second;
2088	if(!newCR \|\| !newPR)
2089	{
2090	G4cerr<<"*G4QIonIonCollision::Breed:CR="<<newCR<<",PL="<<newPL<<G4endl;
2091	G4Exception("G4QIonIonCollision::Breeder:","72",FatalException,"2-RL");
2092	}
2093	cLeft->SetPDGCode(-newCL); // Reset the left quark of curString
2094	cRight->SetPDGCode(newCR); // Reset the right quark of curString
2095	Left->SetPDGCode(-newPL); // Reset the left quark of reString
2096	Right->SetPDGCode(newPR); // Reset the right quark of reString
2097	break; // Break out of the reString internal LOOP
2098	}
2099	} // End of IF(possible reduction)
2100	}
2101	// Check StringsCanBeFused: all QaQ+QaQ(22), single QaQ+QDiQ/aQaDtQ(23/32),
2102	// double QaQ+DiQaDiQ(24/42), QDiQ+aDiQaQ(34/43)
2103	#ifdef debug
2104	G4cout<<"G4QInel::Breed: TryFuse/w #"<<restr<<",q="<<rPDG<<",a="<<aPDG<<G4endl;
2105	#endif
2106	G4int PLS=PLT+PRT;
2107	if( (LS==2 && PLS==2) \|\| // QaQ+QaQ always to DiQaDiQ
2108	( ( (LS==2 && PLS==3) \|\| (LS==3 && PLS==2) ) &&// QaQ w QDiQ/aQaDiQ(single)
2109	( (aPDG> 7 && (-dPDG==aPDG/10 \|\| -dPDG==aPDG%10) ) \|\| // cAQ w DQ
2110	(dPDG> 7 && (-aPDG==dPDG/10 \|\| -aPDG==dPDG%10) ) \|\| // AQ w cDQ
2111	(rPDG<-7 && (qPDG==(-rPDG)/10 \|\| qPDG==(-rPDG)%10) ) \|\| // cQ w ADQ
2112	(qPDG<-7 && (rPDG==(-qPDG)/10 \|\| rPDG==(-qPDG)%10) ) // Q w cADQ
2113	//\|\| (aPDG< 0 && -aPDG==qPDG) \|\| (dPDG< 0 && -dPDG==rPDG) // aQ w Q
2114	)
2115	) \|\| // -----------------------
2116	( ( LS==2 && PLS==4 &&// QaQ w DiQaDiQ (double)
2117	(aPDG> 7 && (-dPDG == aPDG/10 \|\| -dPDG == aPDG%10) ) &&
2118	(rPDG<-7 && (qPDG==(-rPDG)/10 \|\| qPDG==(-rPDG)%10) )
2119	) \|\|
2120	( LS==4 && PLS==2 &&// DiQaDiQ w QaQ (double)
2121	(dPDG> 7 && (-aPDG == dPDG/10 \|\| -aPDG == dPDG%10) ) &&
2122	(qPDG<-7 && (rPDG==(-qPDG)/10 \|\| rPDG==(-qPDG)%10) )
2123	)
2124	) \|\| // -----------------------
2125	( LS==3 && PLS==3 &&// QDiQ w aDiQaQ (double)
2126	( (aPDG> 7 && (-dPDG == aPDG/10 \|\| -dPDG == aPDG%10) &&
2127	qPDG<-7 && (rPDG==(-qPDG)/10 \|\| rPDG==(-qPDG)%10)
2128	) \|\|
2129	(dPDG> 7 && (-aPDG == dPDG/10 \|\| -aPDG == dPDG%10) &&
2130	rPDG<-7 && (dPDG==(-rPDG)/10 \|\| dPDG==(-rPDG)%10)
2131	)
2132	)
2133	)
2134	)
2135	{
2136	G4LorentzVector reString4M = reString->Get4Momentum();
2137	G4ThreeVector reV = reString4M.vect()/reString4M.e();
2138	G4double dV=(curV-reV).mag2(); // Squared difference between relVelocities
2139	#ifdef debug
2140	G4cout<<"G4QIonIonCollision::Breeder: StringCand#"<<restr<<", q="<<rPDG
2141	<<", a="<<aPDG<<", L="<<uPDG<<", R="<<mPDG<<",dV="<<dV<<G4endl;
2142	#endif
2143	if(dV < Vmin)
2144	{
2145	Vmin=dV;
2146	fustr=restr;
2147	MPS=PLS;
2148	}
2149	}
2150	}
2151	if(reduce==2) // String mutual reduction happened
2152	{
2153	#ifdef debug
2154	G4cout<<"-G4QIonIonCollision::Breed:Reduced #"<<astring<<" & #"<<restr<<G4endl;
2155	#endif
2156	astring--; // String was QCreduced using another String
2157	continue; // Repeat fragmentation of the same string
2158	}
2159	if(fustr) // The partner was found -> fuse strings
2160	{
2161	#ifdef debug
2162	G4cout<<"G4QInel::Breeder: StPartner#"<<fustr<<", LT="<<LT<<",RT="<<RT<<G4endl;
2163	#endif
2164	G4QString* fuString=strings[fustr];
2165	G4QParton* Left=fuString->GetLeftParton();
2166	G4QParton* Right=fuString->GetRightParton();
2167	G4int uPDG=Left->GetPDGCode();
2168	G4int mPDG=Right->GetPDGCode();
2169	G4int rPDG=uPDG;
2170	G4int aPDG=mPDG;
2171	if(aPDG<-99 \|\| (aPDG>0 && aPDG<7) \|\| rPDG>99 \|\| (rPDG<0 && rPDG>-7))
2172	{
2173	swap=rPDG;
2174	rPDG=aPDG;
2175	aPDG=swap;
2176	}
2177	if(aPDG > 99) aPDG/=100;
2178	if(rPDG <-99) rPDG=-(-rPDG)/100;
2179	// Check that the strings can fuse (no anti-diquarks assumed)
2180	G4LorentzVector resL4M(0.,0.,0.,0.);
2181	G4LorentzVector resR4M(0.,0.,0.,0.);
2182	G4LorentzVector L4M=cLeft->Get4Momentum();
2183	G4LorentzVector R4M=cRight->Get4Momentum();
2184	#ifdef debug
2185	G4cout<<"G4QIonIonCollision::Breeder:BeforeFuDir,sL="<<sPDG<<",nR="<<nPDG
2186	<<",uL="<<uPDG<<",mR="<<mPDG<<",L4M="<<L4M<<",R4M="<<R4M<<G4endl;
2187	#endif
2188	G4LorentzVector PL4M=Left->Get4Momentum();
2189	G4LorentzVector PR4M=Right->Get4Momentum();
2190	fusionDONE=AnnihilationOrder(LS, MPS, uPDG, mPDG, sPDG, nPDG);
2191	if (fusionDONE>0)
2192	{
2193	if(fusionDONE>1) // Anihilation algorithm
2194	{
2195	if ( (uPDG<0 \|\| nPDG<0) && -uPDG==nPDG ) Left->SetPDGCode(sPDG);
2196	else if( (mPDG<0 \|\| sPDG<0) && -mPDG==sPDG ) Right->SetPDGCode(nPDG);
2197	else if( (uPDG<0 \|\| sPDG<0) && -uPDG==sPDG ) Left->SetPDGCode(nPDG);
2198	else if( (mPDG<0 \|\| nPDG<0) && -mPDG==nPDG ) Right->SetPDGCode(sPDG);
2199	}
2200	{
2201	Left->SetPDGCode(SumPartonPDG(uPDG, sPDG));
2202	Right->SetPDGCode(SumPartonPDG(mPDG, nPDG));
2203	}
2204	Left->Set4Momentum(L4M+PL4M);
2205	Right->Set4Momentum(R4M+PR4M);
2206	#ifdef debug
2207	G4cout<<"G4QIonIonCollision::Breeder:LL/RR s4M="<<fuString->Get4Momentum()
2208	<<",S="<<L4M+PL4M+R4M+PR4M<<", L="<<Left->Get4Momentum()<<", R="
2209	<<Right->Get4Momentum()<<G4endl;
2210	#endif
2211	}
2212	else if(fusionDONE<0)
2213	{
2214	Left->SetPDGCode(SumPartonPDG(uPDG, nPDG));
2215	Left->Set4Momentum(L4M+PR4M);
2216	Right->SetPDGCode(SumPartonPDG(mPDG, sPDG));
2217	Right->Set4Momentum(R4M+PL4M);
2218	#ifdef debug
2219	G4cout<<"G4QIonIonCollision::Breeder:LR/RL s4M="<<fuString->Get4Momentum()
2220	<<",S="<<L4M+PL4M+R4M+PR4M<<", L="<<Left->Get4Momentum()<<", R="
2221	<<Right->Get4Momentum()<<G4endl;
2222	#endif
2223	}
2224	#ifdef debug
2225	else G4cout<<"-Warning-G4QIonIonCollision::Breeder: WrongStringFusion"<<G4endl;
2226	#endif
2227	#ifdef edebug
2228	G4cout<<"#EMC#G4QIonIonCollision::Breed:StringFused,F="<<fusionDONE<<",L="<<L4M
2229	<<",R="<<R4M<<",pL="<<PL4M<<",pR="<<PR4M<<",nL="<<Left->Get4Momentum()
2230	<<",nR="<<Right->Get4Momentum()<<",S="<<fuString->Get4Momentum()<<G4endl;
2231	#endif
2232	if(fusionDONE)
2233	{
2234	#ifdef debug
2235	G4cout<<"###G4QIonIonCollision::Breed: Str#"<<astring<<" fused/w Str#"<<fustr
2236	<<"->"<<fuString->Get4Momentum()<<fuString->Get4Momentum().m()<<G4endl;
2237	#endif
2238	continue; // @@ killing of the curString is excluded
2239	}
2240	}
2241	#ifdef debug
2242	else
2243	{
2244
2245	G4cerr<<"*G4QIonIonCollision::Breed:NoPart*M="<<curString->Get4Momentum().m()
2246	<<", F="<<fusionDONE<<", LPDG="<<curString->GetLeftParton()->GetPDGCode()
2247	<<", RPDG="<<curString->GetRightParton()->GetPDGCode()<<G4endl;
2248	}
2249	#endif
2250	}
2251	if(!fusionDONE) // Fusion of theString failed, try hadrons
2252	{
2253	G4int nHadr=theResult->size(); // #of collected resulting hadrons upToNow
2254	#ifdef debug
2255	G4cout<<"+++4QFragmentation::Breeder:TryHadr M="<<curString->Get4Momentum().m()
2256	<<", nH="<<nHadr<<", LPDG="<<curString->GetLeftParton()->GetPDGCode()
2257	<<", RPDG="<<curString->GetRightParton()->GetPDGCode()<<G4endl;
2258	#endif
2259	// The only solution now is to try fusion with the secondary hadron
2260	while( (nHadr=theResult->size()) && !theHadrons)
2261	{
2262	#ifdef edebug
2263	for(G4int i=0; i<nHadr; i++)
2264	{
2265	G4LorentzVector h4M=(*theResult)[i]->Get4Momentum();
2266	G4int hPDG=(*theResult)[i]->GetPDGCode();
2267	G4QContent hQC=(*theResult)[i]->GetQC();
2268	G4cout<<"-EMC-G4QInel::Breed:H#"<<i<<",4M="<<h4M<<hQC<<",PDG="<<hPDG<<G4endl;
2269	}
2270	#endif
2271	G4int fusDONE=0; // String+Hadron fusion didn't happen
2272	G4int fuhad=-1; // The found hadron index
2273	G4int newPDG=0; // PDG ofTheParton afterMerging with Hadr
2274	G4int secPDG=0; // Second PDG oParton afterMerging w/Hadr
2275	G4double maM2=-DBL_MAX; // Prototype of the max ResultingStringM2
2276	G4LorentzVector selH4M(0.,0.,0.,0.); // 4-mom of the selected hadron
2277	G4QHadron* selHP=0; // Pointer to the used hadron for erasing
2278	G4QString* cString=strings[astring]; // Must be the last string by definition
2279	G4LorentzVector cString4M = cString->Get4Momentum();
2280	G4QParton* cLeft=cString->GetLeftParton();
2281	G4QParton* cRight=cString->GetRightParton();
2282	G4int sumT=cLeft->GetType()+cRight->GetType();
2283	G4int sPDG=cLeft->GetPDGCode();
2284	G4int nPDG=cRight->GetPDGCode();
2285	G4int cMax=0; // Both or only one cuark can merge
2286	for (G4int reh=0; reh < nHadr; reh++)
2287	{
2288	G4QHadron* curHadr=(*theResult)[reh];
2289	G4QContent curQC=curHadr->GetQC(); // Quark content of the hadron
2290	G4int partPDG1=curQC.AddParton(sPDG);
2291	G4int partPDG2=curQC.AddParton(nPDG);
2292	#ifdef debug
2293	G4cout<<"G4QIonIonCollision::Breeder: Hadron # "<<reh<<", QC="<<curQC
2294	<<", P1="<<partPDG1<<", P2="<<partPDG2<<G4endl;
2295	#endif
2296	if(partPDG1 \|\| partPDG2) // Hadron can merge at least w/one parton
2297	{
2298	G4int cCur=1;
2299	if(sumT>3 && partPDG1 && partPDG2) cCur=2;
2300	G4LorentzVector curHadr4M = curHadr->Get4Momentum();
2301	G4double M2=(cString4M+curHadr4M).m2();// SquaredMass of theResultingString
2302	#ifdef debug
2303	G4cout<<"G4QIonIonCollision::Breeder:INHadron#"<<reh<<",M2="<<M2<<G4endl;
2304	#endif
2305	if( (sumT<4 \|\| cCur>=cMax) && M2 > maM2) // DeepAnnihilation for DiQ-aDiQ
2306	{
2307	maM2=M2;
2308	fuhad=reh;
2309	if(partPDG1)
2310	{
2311	fusDONE= 1;
2312	newPDG=partPDG1;
2313	if(partPDG2)
2314	{
2315	secPDG=partPDG2;
2316	cMax=cCur;
2317	}
2318	}
2319	else
2320	{
2321	fusDONE=-1;
2322	newPDG=partPDG2;
2323	}
2324	#ifdef debug
2325	G4cout<<"G4QInel::Br:Selected,P1="<<partPDG1<<",P2="<<partPDG2<<G4endl;
2326	#endif
2327	selH4M=curHadr4M;
2328	selHP=curHadr;
2329	} // End of IF(update selection)
2330	} // End of IF(HadronCanMergeWithTheString)
2331	} // End of the LOOP over Hadrons
2332	#ifdef debug
2333	G4cout<<"G4QIonIonCollision::Breeder: fuh="<<fuhad<<",fus="<<fusDONE<<G4endl;
2334	#endif
2335	if(fuhad>-1) // The partner was found -> fuse H&S
2336	{
2337	if (fusDONE>0) // Fuse hadron with the left parton
2338	{
2339	cLeft->SetPDGCode(newPDG);
2340	G4LorentzVector newLeft=cLeft->Get4Momentum()+selH4M;
2341	cLeft->Set4Momentum(newLeft);
2342	if(secPDG && cMax>1)
2343	{
2344	#ifdef debug
2345	G4cout<<"G4QInel::Br:TryReduce, nPDG="<<newPDG<<",sPDG="<<secPDG<<G4endl;
2346	#endif
2347	cLeft->ReduceDiQADiQ(cLeft, cRight);
2348	}
2349	#ifdef debug
2350	G4cout<<"G4QIonIonCollision::Breed: Left, s4M="<<curString->Get4Momentum()
2351	<<", L4M="<<newLeft<<", R4M="<<cRight->Get4Momentum()<<", h4M="
2352	<<selH4M<<", newL="<<newPDG<<", oldR="<<cRight->GetPDGCode()<<G4endl;
2353	#endif
2354	}
2355	else if(fusDONE<0) // Fuse hadron with the right parton
2356	{
2357	cRight->SetPDGCode(newPDG);
2358	G4LorentzVector newRight=cRight->Get4Momentum()+selH4M;
2359	cRight->Set4Momentum(newRight);
2360	#ifdef debug
2361	G4cout<<"G4QIonIonCollision::Breed: Right, s4M="<<curString->Get4Momentum()
2362	<<", L4M="<<cLeft->Get4Momentum()<<", R4M="<<newRight<<", h4M="
2363	<<selH4M<<", newR="<<newPDG<<", oldL="<<cLeft->GetPDGCode()<<G4endl;
2364	#endif
2365	}
2366	#ifdef debug
2367	else G4cout<<"-G4QIonIonCollision::Breed: Wrong String+HadronFusion"<<G4endl;
2368	#endif
2369	#ifdef debug
2370	if(fusDONE) G4cout<<"####G4QIonIonCollision::Breeder: String #"<<astring
2371	<<" is fused with Hadron #"<<fuhad
2372	<<", new4Mom="<<curString->Get4Momentum()
2373	<<", M2="<<curString->Get4Momentum().m2()
2374	<<", QC="<<curString->GetQC()<<G4endl;
2375	#endif
2376	}
2377	else
2378	{
2379	#ifdef debug
2380	G4cerr<<"*G4QIonIonCollision::Breed:NoH*,M="<<curString->Get4Momentum().m()
2381	<<", LPDG="<<curString->GetLeftParton()->GetPDGCode()
2382	<<", RPDG="<<curString->GetRightParton()->GetPDGCode()<<G4endl;
2383	// @@ Temporary exception for the future solution
2384	//G4Exception("G4QIonIonCollision::Breed:","72",FatalException,"SHNotFused");
2385	#endif
2386	break; // Breake the While LOOP
2387	} // End of the namespace where both Fusion and reduction have failed
2388	// The fused hadron must be excluded from theResult
2389	#ifdef debug
2390	G4cout<<"G4QIonIonCollision::Breed: before HR, nH="<<theResult->size()<<G4endl;
2391	G4int icon=0; // Loop counter
2392	#endif
2393	G4QHadronVector::iterator ih;
2394	G4bool found=false;
2395	for(ih = theResult->begin(); ih != theResult->end(); ih++)
2396	{
2397	#ifdef debug
2398	G4cout<<"G4QInelast::Breeder:#"<<icon<<", i="<<(*ih)<<", sH="<<selHP<<G4endl;
2399	icon++;
2400	#endif
2401	if((*ih)==selHP)
2402	{
2403	#ifdef debug
2404	G4cout<<"G4QInel::Breed: HadronFound, PDG="<<selHP->GetPDGCode()<<G4endl;
2405	#endif
2406	G4LorentzVector p4M=selHP->Get4Momentum();
2407	curString4M+=p4M;
2408	#ifdef edebug
2409	G4int Chg=selHP->GetCharge();
2410	G4int BaN=selHP->GetBaryonNumber();
2411	curStrChg+=Chg;
2412	curStrBaN+=BaN;
2413	G4cout<<"-EMC->>>>G4QIonIonCollision::Breed: S+=H, 4M="<<curString4M<<",M="
2414	<<curString4M.m()<<", Charge="<<curStrChg<<", B="<<curStrBaN<<G4endl;
2415	#endif
2416	delete selHP; // delete the Hadron
2417	theResult->erase(ih); // erase the Hadron from theResult
2418	found=true;
2419	break; // beak the LOOP over hadrons
2420	}
2421	} // End of the LOOP over hadrons
2422	#ifdef debug
2423	if(!found) G4cout<<"*G4QIonIonCollision::Breed:nH="<<theResult->size()<<G4endl;
2424	#endif
2425	// New attempt of the string decay
2426	theHadrons=curString->FragmentString(true);//! Try to fragment the new String !
2427	#ifdef debug
2428	G4cout<<"G4QInel::Breeder: tH="<<theHadrons<<",nH="<<theResult->size()<<G4endl;
2429	#endif
2430	} // End of the while LOOP over the fusion with hadrons
2431	#ifdef debug
2432	G4cout<<"G4QIonIonCollision::Breed: CanTryToDecay?* nH="<<theHadrons<<", next="
2433	<<next<<" =? nS="<<strings.size()<<", nR="<<theResult->size()<<G4endl;
2434	#endif
2435	if(!theHadrons && next == strings.size() && !(theResult->size()))// TryToDecay
2436	{
2437	G4QContent miQC=curString->GetQC(); // QContent of the Lightest Hadron
2438	G4int miPDG=miQC.GetSPDGCode(); // PDG of the Lightest Hadron
2439	if(miPDG == 10) // ==> Decay Hadron-Chipolino
2440	{
2441	G4QChipolino QCh(miQC); // define theTotalResidual as aChipolino
2442	G4QPDGCode h1QPDG=QCh.GetQPDG1(); // QPDG of the first hadron
2443	G4double h1M =h1QPDG.GetMass();// Mass of the first hadron
2444	G4QPDGCode h2QPDG=QCh.GetQPDG2(); // QPDG of the second hadron
2445	G4double h2M =h2QPDG.GetMass();// Mass of the second hadron
2446	G4double ttM =curString4M.m(); // Real Mass of the Chipolino
2447	if(h1M+h2M<ttM+eps) // Two particles decay of Chipolino
2448	{
2449	G4LorentzVector h14M(0.,0.,0.,h1M);
2450	G4LorentzVector h24M(0.,0.,0.,h2M);
2451	if(std::fabs(ttM-h1M-h2M)<=eps)
2452	{
2453	G4double part1=h1M/(h1M+h2M);
2454	h14M=part1*curString4M;
2455	h24M=curString4M-h14M;
2456	}
2457	else
2458	{
2459	if(!G4QHadron(curString4M).DecayIn2(h14M,h24M))
2460	{
2461	G4cerr<<"***G4QIonIonCollision::Breeder: tM="<<ttM<<"->h1="<<h1QPDG
2462	<<"("<<h1M<<")+h2="<<h1QPDG<<"("<<h2M<<")="<<h1M+h2M<<G4endl;
2463	G4Exception("G4QIonIonCollision::Breeder:","72",FatalException,"QDec");
2464	}
2465	}
2466	G4QHadron* h1H = new G4QHadron(h1QPDG.GetPDGCode(),h14M);
2467	theResult->push_back(h1H); // (delete equivalent)
2468	#ifdef debug
2469	G4LorentzVector f4M=h1H->Get4Momentum();
2470	G4int fPD=h1H->GetPDGCode();
2471	G4int fCg=h1H->GetCharge();
2472	G4int fBN=h1H->GetBaryonNumber();
2473	G4cout<<"-EMC->>G4QIonIonCollision::Breed:String=HadrChiPro1's filled,f4M="
2474	<<f4M<<", fPDG="<<fPD<<", fCg="<<fCg<<", fBN="<<fBN<<G4endl;
2475	#endif
2476	G4QHadron* h2H = new G4QHadron(h2QPDG.GetPDGCode(),h24M);
2477	theResult->push_back(h2H); // (delete equivalent)
2478	#ifdef debug
2479	G4LorentzVector s4M=h2H->Get4Momentum();
2480	G4int sPD=h2H->GetPDGCode();
2481	G4int sCg=h2H->GetCharge();
2482	G4int sBN=h2H->GetBaryonNumber();
2483	G4cout<<"-EMC->>G4QIonIonCollision::Breed:String=HadrChiPro2's filled,s4M="
2484	<<s4M<<", sPDG="<<sPD<<", sCg="<<sCg<<", sBN="<<sBN<<G4endl;
2485	#endif
2486	#ifdef edebug
2487	G4cout<<"-EMC-..Chi..G4QIonIonCollision::Breeder: DecayCHECK, Ch4M="
2488	<<curString4M<<", d4M="<<curString4M-h14M-h24M<<G4endl;
2489	#endif
2490	break; // Go out of the main StringDecay LOOP
2491	}
2492	else
2493	{
2494	G4Quasmon* stringQuasmon = new G4Quasmon(miQC, curString4M);// String->Quas
2495	theQuasmons.push_back(stringQuasmon);
2496	break; // Go out of the main StringDecay LOOP
2497	}
2498	}
2499	else if(miPDG) // Decay Hadron as a Resonans
2500	{
2501	if (miPDG>0 && miPDG%10 < 3) miPDG+=2; // Convert to Resonans
2502	else if(miPDG<0 && (-miPDG)%10< 3) miPDG-=2; // Convert to antiResonans
2503	G4Quasmon Quasm;
2504	G4QHadron* sHad = new G4QHadron(miPDG,curString4M);
2505	G4QHadronVector* tmpQHadVec=Quasm.DecayQHadron(sHad); // It deletes sHad
2506	G4int tmpN=tmpQHadVec->size();
2507	#ifdef debug
2508	G4cout<<"G4QIonIonCollision::Breeder: Decay the Last, Res#H="<<tmpN<<G4endl;
2509	#endif
2510	if(tmpN>1)
2511	{
2512	for(G4int aH=0; aH < tmpN; aH++)
2513	{
2514	theResult->push_back((*tmpQHadVec)[aH]);//TheDecayProdOfHadronDecIsFilled
2515	#ifdef debug
2516	G4QHadron* prodH =(*tmpQHadVec)[aH];
2517	G4LorentzVector p4M=prodH->Get4Momentum();
2518	G4int PDG=prodH->GetPDGCode();
2519	G4int Chg=prodH->GetCharge();
2520	G4int BaN=prodH->GetBaryonNumber();
2521	G4cout<<"-EMC->>G4QIonIonCollis::Breed:String=Hadr,H#"<<aH<<" filled,4M="
2522	<<p4M<<", PDG="<<PDG<<", Chg="<<Chg<<", BaN="<<BaN<<G4endl;
2523	#endif
2524	}
2525	}
2526	else
2527	{
2528	G4Quasmon* stringQuasmon = new G4Quasmon(miQC, curString4M);// String->Quas
2529	#ifdef debug
2530	G4cout<<"G4QIonIonCollision::Breeder:==> to Quasm="<<miQC<<curString4M
2531	<<", pNuc="<<theProjNucleus<<theProjNucleus.Get4Momentum()<<", tNuc="
2532	<<theTargNucleus<<theTargNucleus.Get4Momentum()<<", NString="
2533	<<strings.size()<<", nR="<<theResult->size()<<", nQ="
2534	<<theQuasmons.size()<<G4endl;
2535	#endif
2536	theQuasmons.push_back(stringQuasmon);
2537	delete sHad;
2538	tmpQHadVec->clear();
2539	delete tmpQHadVec; // WhoCallsDecayQHadron is responsible for clear&delete
2540	break; // Go out of the main StringDecay LOOP
2541	}
2542	tmpQHadVec->clear();
2543	delete tmpQHadVec; // Who calls DecayQHadron is responsible for clear&delete
2544	break; // Go out of the main String Decay LOOP
2545	}
2546	} // End of the DecayOfTheLast
2547	} // End of IF(String-Hadron fusion)
2548	} // End of IF(NO_Hadrons) for String-String and String-Hadron fusion
2549	// The last hope is to CORREC the string, using other strings (ForwardInLOOP)
2550	#ifdef debug
2551	G4cout<<"G4QIonIonCollision::Breeder: theH="<<theHadrons<<"?=0, next="<<next<<G4endl;
2552	#endif
2553	if(!theHadrons && next < strings.size()) // ForwardInLOOP strings exist
2554	{
2555	// @@ string can be not convertable to one hadron (2,0.0,0,2,0) - To be improved
2556	G4QContent miQC=curString->GetQC(); // QContent of the Lightest Hadron
2557	G4int miPDG=miQC.GetSPDGCode();// PDG of the Lightest Hadron
2558	#ifdef debug
2559	G4cout<<">>>G4QIonIonCollision::Breeder: SQC="<<miQC<<", miSPDG="<<miPDG<<G4endl;
2560	#endif
2561	G4double miM=0.; // Prototype of the Mass of the Cur LightestHadron
2562	if(miPDG!=10) miM=G4QPDGCode(miPDG).GetMass(); // Mass of theCurLightestOneHadron
2563	else
2564	{
2565	G4QChipolino QCh(miQC); // define the TotalString as a Chipolino
2566	miM=QCh.GetQPDG1().GetMass()+QCh.GetQPDG2().GetMass();//MinMass of theStringChipo
2567	}
2568	G4double cM2=curString4M.m2(); // Actual squared mass of the Cur String
2569	#ifdef debug
2570	G4cout<<">>>G4QIonIonCollision::Breeder: minMass="<<miM<<", realM2="<<cM2<<G4endl;
2571	#endif
2572	G4double cM=0.;
2573	if(cM2>0.)
2574	{
2575	cM=std::sqrt(cM2);
2576	if(std::fabs(cM-miM) < eps) // Convert to hadron(2 hadrons) w/o calculations
2577	{
2578	if(miPDG!=10)
2579	{
2580	G4QHadron* sHad = new G4QHadron(miPDG,curString4M);
2581	theResult->push_back(sHad);// Fill the curString as a hadron
2582	#ifdef debug
2583	G4cout<<">>>G4QIonIonCollision::Breeder:S->H="<<miPDG<<curString4M<<G4endl;
2584	#endif
2585	}
2586	else
2587	{
2588	G4QChipolino QCh(miQC); // define TotalResidual as a Chipolino
2589	G4QPDGCode h1QPDG=QCh.GetQPDG1(); // QPDG of the first hadron
2590	G4double h1M =h1QPDG.GetMass(); // Mass of the first hadron
2591	G4QPDGCode h2QPDG=QCh.GetQPDG2(); // QPDG of the second hadron
2592	G4double h2M =h2QPDG.GetMass(); // Mass of the second hadron
2593	G4double pt1 =h1M/(h1M+h2M); // 4-mom part of the first hadron
2594	G4LorentzVector h14M=pt1*curString4M; // 4-mom of the first hadron
2595	G4LorentzVector h24M=curString4M-h14M;// 4-mom of the second hadron
2596	G4QHadron* h1H = new G4QHadron(h1QPDG.GetPDGCode(),h14M);
2597	theResult->push_back(h1H); // (delete equivalent)
2598	#ifdef debug
2599	G4LorentzVector f4M=h1H->Get4Momentum();
2600	G4int fPD=h1H->GetPDGCode();
2601	G4int fCg=h1H->GetCharge();
2602	G4int fBN=h1H->GetBaryonNumber();
2603	G4cout<<"-EMC->>G4QIonIonCollision::Breed:Str=2HadrAR Prod-F is filled, f4M="
2604	<<f4M<<", fPDG="<<fPD<<", fCg="<<fCg<<", fBN="<<fBN<<G4endl;
2605	#endif
2606	G4QHadron* h2H = new G4QHadron(h2QPDG.GetPDGCode(),h24M);
2607	theResult->push_back(h2H); // (delete equivalent)
2608	#ifdef debug
2609	G4LorentzVector s4M=h2H->Get4Momentum();
2610	G4int sPD=h2H->GetPDGCode();
2611	G4int sCg=h2H->GetCharge();
2612	G4int sBN=h2H->GetBaryonNumber();
2613	G4cout<<"-EMC->>G4QIonIonCollision::Breed:Str=2HadrAR Prod-S is filled, s4M="
2614	<<s4M<<", sPDG="<<sPD<<", sCg="<<sCg<<", sBN="<<sBN<<G4endl;
2615	#endif
2616	}
2617	continue; // Continue the LOOP over the curStrings
2618	}
2619	else // Try to recover (+/-) to min Mass
2620	{
2621	G4ThreeVector cP=curString4M.vect(); // Momentum of the curString
2622	G4double cE=curString4M.e(); // Energy of the curString
2623	G4ThreeVector curV=cP/cE; // curRelativeVelocity
2624	G4double miM2=miM*miM;
2625	G4int restr=0; // To use beyon the LOOP for printing
2626	G4int fustr=0; // Selected String-partner (0 = NotFound)
2627	G4double selX=0.; // Selected value of x
2628	G4double maD=-DBL_MAX; // Maximum Free Mass
2629	G4double Vmin=DBL_MAX; // min Velocity Distance
2630	G4LorentzVector s4M(0.,0.,0.,0.); // Selected 4-mom of the hadron
2631	#ifdef debug
2632	G4cout<<"G4QIonIonCollision::Breeder: TryRecover, cV="<<curV<<G4endl;
2633	#endif
2634	nOfStr=strings.size();
2635	for(restr=next; restr < nOfStr; ++restr) if(restr != astring)
2636	{
2637	G4QString* pString=strings[restr];
2638	G4LorentzVector p4M=pString->Get4Momentum();
2639	G4ThreeVector pP=p4M.vect(); // Momentum of the partnerString
2640	G4double pE=p4M.e(); // Energy of the partnerString
2641	G4double D2=cE*pE-cP.dot(pP);
2642	G4double pM2=p4M.m2();
2643	G4double dM4=pM2*(miM2-cM2);
2644	G4double D=D2*D2+dM4;
2645	G4double x=-1.; // Bad preexpectation
2646	if(D >= 0. && pM2>.01) x=(std::sqrt(D)-D2)/pM2; // what we should get from p
2647	#ifdef debug
2648	else G4cout<<"G4QIonIonCollis::Breed:D="<<D<<",D2="<<D2<<",dM="<<dM4<<G4endl;
2649	G4cout<<"G4QIonIonCollision::Breed: pM2="<<pM2<<",D2="<<D2<<",x="<<x<<G4endl;
2650	#endif
2651	if(x > 0. && x < 1.) // We are getting x part of p4M
2652	{
2653	G4QContent pQC=pString->GetQC(); // Quark Content of The Partner
2654	G4int pPDG=pQC.GetSPDGCode();// PDG of The Lightest Hadron for the Partner
2655	G4double pM=0.; // Mass of the LightestHadron
2656	if(pPDG==10)
2657	{
2658	G4QChipolino QCh(pQC); // define the TotalString as a Chipolino
2659	pM=QCh.GetQPDG1().GetMass()+QCh.GetQPDG2().GetMass();// Mass of Chipolino
2660	}
2661	else pM=G4QPDGCode(pPDG).GetMass();// Mass of theLightestHadron for Partner
2662	G4double rM=std::sqrt(pM2); // Real mass of the string-partner
2663	G4double delta=(1.-x)*rM-pM;// @@ Minimum CM disterbance measure
2664	if(delta > 0. && delta > maD)
2665	{
2666	maD=delta;
2667	#ifdef debug
2668	G4cout<<"G4QIonIonCollision::Breed: Subtr,S#"<<restr<<",d="<<maD<<G4endl;
2669	#endif
2670	fustr=restr;
2671	selX=x;
2672	s4M=p4M;
2673	}
2674	}
2675	else if(x <= 0.) // We are adding to p4M, so use RelVelocity
2676	{
2677	G4ThreeVector pV=pP/pE; // curRelativeVelocity
2678	G4double dV=(curV-pV).mag2();// SquaredDifferenceBetweenRelVelocities
2679	if(dV < Vmin)
2680	{
2681	#ifdef debug
2682	G4cout<<"G4QIonIonCollision::Breed: FreeAdd,S#"<<restr<<",x="<<x<<G4endl;
2683	#endif
2684	Vmin=dV;
2685	fustr=restr;
2686	selX=x;
2687	s4M=p4M;
2688	}
2689	}
2690	#ifdef debug
2691	G4cout<<"G4QIonIonCollision::Breed:EndOfLOOP r="<<restr<<"<"<<nOfStr<<G4endl;
2692	#endif
2693	} // End of the LOOP over string-partners for Correction
2694	#ifdef debug
2695	G4cout<<"G4QIonIonCollision::Breeder: AfterLOOP fustr="<<fustr<<G4endl;
2696	#endif
2697	if(fustr)
2698	{
2699	#ifdef edebug
2700	G4LorentzVector sum4M=s4M+curString4M;
2701	G4cout<<"G4QIonIonCollision::Breeder: Found Sum4M="<<sum4M<<G4endl;
2702	#endif
2703	G4QString* pString=strings[fustr];
2704	curString4M+=selX*s4M;
2705	if(std::abs(miPDG)%10 > 2) // Decay String-Hadron-Resonance
2706	{
2707	G4Quasmon Quasm;
2708	G4QHadron* sHad = new G4QHadron(miPDG,curString4M);
2709	G4QHadronVector* tmpQHadVec=Quasm.DecayQHadron(sHad); // It deletes sHad
2710	#ifdef debug
2711	G4cout<<"G4QIonIonCollision::Breed:DecStH,nH="<<tmpQHadVec->size()<<G4endl;
2712	#endif
2713	for(unsigned aH=0; aH < tmpQHadVec->size(); aH++)
2714	{
2715	theResult->push_back((*tmpQHadVec)[aH]);//TheDecayProdOfHadron is filled
2716	#ifdef debug
2717	G4QHadron* prodH =(*tmpQHadVec)[aH];
2718	G4LorentzVector p4M=prodH->Get4Momentum();
2719	G4int PDG=prodH->GetPDGCode();
2720	G4int Chg=prodH->GetCharge();
2721	G4int BaN=prodH->GetBaryonNumber();
2722	G4cout<<"-EMC->>G4QIonIonCollision::Breed:St=Had,pH#"<<aH<<" filled, 4M="
2723	<<p4M<<", PDG="<<PDG<<", Chg="<<Chg<<", BaN="<<BaN<<G4endl;
2724	#endif
2725	}
2726	tmpQHadVec->clear();
2727	delete tmpQHadVec; // Who calls DecayQHadron is responsibleRorClear&Delete
2728	}
2729	else if(miPDG == 10) // ==> Decay Hadron-Chipolino
2730	{
2731	G4QChipolino QCh(miQC); // define theTotalResid as aChipolino
2732	G4QPDGCode h1QPDG=QCh.GetQPDG1(); // QPDG of the first hadron
2733	G4double h1M =h1QPDG.GetMass();// Mass of the first hadron
2734	G4QPDGCode h2QPDG=QCh.GetQPDG2(); // QPDG of the second hadron
2735	G4double h2M =h2QPDG.GetMass();// Mass of the second hadron
2736	G4double ttM =curString4M.m(); // Real Mass of the Chipolino
2737	if(h1M+h2M<ttM+eps) // Two particles decay of Chipolino
2738	{
2739	G4LorentzVector h14M(0.,0.,0.,h1M);
2740	G4LorentzVector h24M(0.,0.,0.,h2M);
2741	if(std::fabs(ttM-h1M-h2M)<=eps)
2742	{
2743	G4double part1=h1M/(h1M+h2M);
2744	h14M=part1*curString4M;
2745	h24M=curString4M-h14M;
2746	}
2747	else
2748	{
2749	if(!G4QHadron(curString4M).DecayIn2(h14M,h24M))
2750	{
2751	G4cerr<<"***G4QIonIonCollision::Breeder: tM="<<ttM<<"->h1="<<h1QPDG
2752	<<"("<<h1M<<")+h2="<<h1QPDG<<"("<<h2M<<")="<<h1M+h2M<<G4endl;
2753	G4Exception("G4QIonIonCollision::Breeder:","72",FatalException,"CD");
2754	}
2755	}
2756	G4QHadron* h1H = new G4QHadron(h1QPDG.GetPDGCode(),h14M);
2757	theResult->push_back(h1H); // (delete equivalent)
2758	#ifdef debug
2759	G4LorentzVector f4M=h1H->Get4Momentum();
2760	G4int fPD=h1H->GetPDGCode();
2761	G4int fCg=h1H->GetCharge();
2762	G4int fBN=h1H->GetBaryonNumber();
2763	G4cout<<"-EMC->>>G4QIonIonCollision::Breed:Str=Hadr Prod-F's filled,f4M="
2764	<<f4M<<", fPDG="<<fPD<<", fCg="<<fCg<<", fBN="<<fBN<<G4endl;
2765	#endif
2766	G4QHadron* h2H = new G4QHadron(h2QPDG.GetPDGCode(),h24M);
2767	theResult->push_back(h2H); // (delete equivalent)
2768	#ifdef debug
2769	G4LorentzVector s4M=h2H->Get4Momentum();
2770	G4int sPD=h2H->GetPDGCode();
2771	G4int sCg=h2H->GetCharge();
2772	G4int sBN=h2H->GetBaryonNumber();
2773	G4cout<<"-EMC->>>G4QIonIonCollision::Breed:Str=Hadr Prod-S's filled,s4M="
2774	<<s4M<<", sPDG="<<sPD<<", sCg="<<sCg<<", sBN="<<sBN<<G4endl;
2775	#endif
2776	#ifdef edebug
2777	G4cout<<"-EMC-Chipo.G4QIonIonCollision::Breed:DecCHECK,c4M="<<curString4M
2778	<<", ChQC="<<miQC<<", d4M="<<curString4M-h14M-h24M<<G4endl;
2779	#endif
2780	}
2781	else
2782	{
2783	G4cerr<<"***G4QInel::Breeder: tM="<<ttM<<miQC<<"->h1="<<h1QPDG<<"(" <<h1M
2784	<<")+h2="<<h1QPDG<<"("<<h2M<<") = "<<h1M+h2M<<G4endl;
2785	G4Exception("G4QIonIonCollision::Breeder:","72",FatalException,"ChiDec");
2786	}
2787	}
2788	else
2789	{
2790	G4QHadron* sHad = new G4QHadron(miPDG,curString4M);
2791	theResult->push_back(sHad); // The original string-hadron is filled
2792	#ifdef debug
2793	G4cout<<"-EMC->>>G4QIonIonCollision::Breeder:Str=Hadr Filled, 4M="
2794	<<curString4M<<", PDG="<<miPDG<<G4endl;
2795	#endif
2796	}
2797	G4double corF=1-selX;
2798	G4QParton* Left=pString->GetLeftParton();
2799	G4QParton* Right=pString->GetRightParton();
2800	Left->Set4Momentum(corF*Left->Get4Momentum());
2801	Right->Set4Momentum(corF*Right->Get4Momentum());
2802	#ifdef edebug
2803	G4cout<<"-EMC-...Cor...G4QIonIonCollision::Breeder:CorCHECK Sum="<<sum4M
2804	<<" =? "<<curString4M+pString->Get4Momentum()<<", M="<<miM<<" ?= "
2805	<<curString4M.m()<<G4endl;
2806	#endif
2807	#ifdef debug
2808	G4cout<<">>>G4QIonIonCollision::Breeder:Corrected String->Hadr="<<miPDG
2809	<<curString4M<<" by String #"<<fustr<<G4endl;
2810	#endif
2811	continue; // Continue the LOOP over the curStrings
2812	} // End of Found combination for the String on string correction
2813	} // End of the Try-to-recover String+String correction algorithm
2814	} // End of IF(CM2>0.)
2815	} // End of IF(Can try to correct by String-String)
2816	#ifdef debug
2817	else G4cerr<<"*G4QIonIonCollision::Breed:No SSCorrection**, next="<<next<<G4endl;
2818	#endif
2819	// ------------ At this point we can reduce the 3/-3 meson to 1/-1 meson ------------
2820	G4QParton* lpcS=curString->GetLeftParton();
2821	G4QParton* rpcS=curString->GetRightParton();
2822	G4int lPDGcS=lpcS->GetPDGCode();
2823	G4int rPDGcS=rpcS->GetPDGCode();
2824	if (lPDGcS==3 && rPDGcS==-3)
2825	{
2826	lpcS->SetPDGCode( 1);
2827	rpcS->SetPDGCode(-1);
2828	}
2829	else if(lPDGcS==-3 && rPDGcS==3)
2830	{
2831	lpcS->SetPDGCode(-1);
2832	rpcS->SetPDGCode( 1);
2833	}
2834	// -------- Now the only hope is Correction, using the already prodused Hadrons -----
2835	G4int nofRH=theResult->size(); // #of resulting Hadrons
2836	#ifdef debug
2837	G4cout<<"G4QIonIonCollision::Breeder: theH="<<theHadrons<<", #OfH="<<nofRH<<G4endl;
2838	#endif
2839	if(!theHadrons && nofRH) // Hadrons are existing for SH Correction
2840	{
2841	#ifdef debug
2842	G4cerr<<"!G4QIonIonCollision::Breeder:CanTrySHCor, nH="<<theResult->size()<<G4endl;
2843	#endif
2844	// @@ string can be not convertable to one hadron (2,0.0,0,2,0) - To be improved
2845	G4QContent miQC=curString->GetQC(); // QContent of the Lightest Hadron
2846	G4int miPDG=miQC.GetSPDGCode(); // PDG of the Lightest Hadron
2847	G4double miM=0.; // Prototype ofMass of theCurLightestHadron
2848	if(miPDG==10) // Mass of the Cur Lightest ChipolinoHadron
2849	{
2850	G4QChipolino QCh(miQC); // define the TotalString as a Chipolino
2851	miM=QCh.GetQPDG1().GetMass()+QCh.GetQPDG2().GetMass();//MinMass of theStringChipo
2852	}
2853	else miM=G4QPDGCode(miPDG).GetMass(); // Mass of the Cur Lightest Hadron
2854	G4double spM=0.; // Mass of the selected Hadron-Partner
2855	G4ThreeVector cP=curString4M.vect(); // Momentum of the curString
2856	G4double cE=curString4M.e(); // Energy of the curString
2857	G4ThreeVector curV=cP/cE; // curRelativeVelocity
2858	G4int reha=0; // Hadron # to use beyon the LOOP
2859	G4int fuha=0; // Selected Hadron-partner (0 = NotFound)
2860	G4double dMmin=DBL_MAX; // min Excess of the mass
2861	G4LorentzVector s4M(0.,0.,0.,0.); // Selected 4-mom of the Hadron+String
2862	G4double sM=0.; // Selected Mass of the Hadron+String
2863	for (reha=next; reha < nofRH; reha++) // LOOP over already collected hadrons
2864	{
2865	G4QHadron* pHadron=(*theResult)[reha];// Pointer to the current Partner-Hadron
2866	G4LorentzVector p4M=pHadron->Get4Momentum();
2867	G4double pM=p4M.m(); // Mass of the Partner-Hadron
2868	G4LorentzVector t4M=p4M+curString4M; // Total momentum of the compound
2869	G4double tM2=t4M.m2(); // Squared total mass of the compound
2870	if(tM2 >= sqr(pM+miM+eps)) // Condition of possible correction
2871	{
2872	G4double tM=std::sqrt(tM2); // Mass of the Hadron+String compound
2873	G4double dM=tM-pM-miM; // Excess of the compound mass
2874	if(dM < dMmin)
2875	{
2876	dMmin=dM;
2877	fuha=reha;
2878	spM=pM;
2879	s4M=t4M;
2880	sM=tM;
2881	}
2882	}
2883	} // End of the LOOP over string-partners for Correction
2884	if(fuha) // The hadron-partner was found
2885	{
2886	G4QHadron* pHadron=(*theResult)[fuha];// Necessary for update
2887	G4LorentzVector mi4M(0.,0.,0.,miM); // Prototype of the new String=Hadron
2888	if(miM+spM<sM+eps) // Decay into CorrectedString+theSameHadron
2889	{
2890	G4LorentzVector sp4M(0.,0.,0.,spM);
2891	if(std::fabs(sM-miM-spM)<=eps)
2892	{
2893	G4double part1=miM/(miM+spM);
2894	mi4M=part1*s4M;
2895	sp4M=s4M-mi4M;
2896	}
2897	else
2898	{
2899	if(!G4QHadron(s4M).DecayIn2(mi4M,sp4M))
2900	{
2901	G4cerr<<"**G4QIonIonCollision::Breeder: SH*, tM="<<sM<<"->h1=("<<miPDG
2902	<<")"<<miM<<" + h2="<<spM<<" = "<<miM+spM<<G4endl;
2903	G4Exception("G4QIonIonCollision::Breeder:","72",FatalException,"SHChiDec");
2904	}
2905	}
2906	pHadron->Set4Momentum(sp4M);
2907	#ifdef debug
2908	G4cout<<"-EMC->...>G4QIonIonCollision::Breed: H# "<<fuha<<" is updated, new4M="
2909	<<sp4M<<G4endl;
2910	#endif
2911	}
2912	else
2913	{
2914	G4cerr<<"***G4QInel::Breeder: HS Failed, tM="<<sM<<"->h1M=("<<miPDG<<")"<<miM
2915	<<"+h2M="<<spM<<" = "<<miM+spM<<G4endl;
2916	G4Exception("G4QIonIonCollision::Breeder:","72",FatalException,"HSChipoliDec");
2917	}
2918	if(std::abs(miPDG)%10 > 2) // Decay Hadron-Resonans
2919	{
2920	G4Quasmon Quasm;
2921	G4QHadron* sHad = new G4QHadron(miPDG,mi4M);
2922	G4QHadronVector* tmpQHadVec=Quasm.DecayQHadron(sHad); // It deletes sHad
2923	#ifdef debug
2924	G4cout<<"G4QInelast::Breeder: HS DecStrHad, nH="<<tmpQHadVec->size()<<G4endl;
2925	#endif
2926	for(unsigned aH=0; aH < tmpQHadVec->size(); aH++)
2927	{
2928	theResult->push_back((*tmpQHadVec)[aH]);// TheDecayProductOfTheHadronIsFilled
2929	#ifdef debug
2930	G4QHadron* prodH =(*tmpQHadVec)[aH];
2931	G4LorentzVector p4M=prodH->Get4Momentum();
2932	G4int PDG=prodH->GetPDGCode();
2933	G4int Chg=prodH->GetCharge();
2934	G4int BaN=prodH->GetBaryonNumber();
2935	G4cout<<"-EMC->>>G4QIonIonCollision::Breed:Str+Hadr PrH#"<<aH<<" filled, 4M="
2936	<<p4M<<", PDG="<<PDG<<", Chg="<<Chg<<", BaN="<<BaN<<G4endl;
2937	#endif
2938	}
2939	tmpQHadVec->clear();
2940	delete tmpQHadVec; // Who calls DecayQHadron is responsible for clear & delete
2941	}
2942	else if(miPDG == 10) // ==> Decay Hadron-Chipolino
2943	{
2944	G4QChipolino QCh(miQC); // define the TotalResidual as a Chipolino
2945	G4QPDGCode h1QPDG=QCh.GetQPDG1(); // QPDG of the first hadron
2946	G4double h1M =h1QPDG.GetMass();// Mass of the first hadron
2947	G4QPDGCode h2QPDG=QCh.GetQPDG2(); // QPDG of the second hadron
2948	G4double h2M =h2QPDG.GetMass();// Mass of the second hadron
2949	G4double ttM =curString4M.m(); // Real Mass of the Chipolino
2950	if(h1M+h2M<miM+eps) // Two particles decay of Chipolino
2951	{
2952	G4LorentzVector h14M(0.,0.,0.,h1M);
2953	G4LorentzVector h24M(0.,0.,0.,h2M);
2954	if(std::fabs(ttM-h1M-h2M)<=eps)
2955	{
2956	G4double part1=h1M/(h1M+h2M);
2957	h14M=part1*mi4M;
2958	h24M=mi4M-h14M;
2959	}
2960	else
2961	{
2962	if(!G4QHadron(mi4M).DecayIn2(h14M,h24M))
2963	{
2964	G4cerr<<"***G4QIonIonCollision::Breeder: HS tM="<<ttM<<"->h1="<<h1QPDG
2965	<<"("<<h1M<<")+h2="<<h1QPDG<<"("<<h2M<<")="<<h1M+h2M<<G4endl;
2966	G4Exception("G4QIonIonCollision::Breeder:","72",FatalException,"ChiDec");
2967	}
2968	}
2969	G4QHadron* h1H = new G4QHadron(h1QPDG.GetPDGCode(),h14M);
2970	theResult->push_back(h1H); // (delete equivalent)
2971	#ifdef debug
2972	G4LorentzVector f4M=h1H->Get4Momentum();
2973	G4int fPD=h1H->GetPDGCode();
2974	G4int fCg=h1H->GetCharge();
2975	G4int fBN=h1H->GetBaryonNumber();
2976	G4cout<<"-EMC->>G4QIonIonCollision::Breed: CorStrHadr Prod-1 is filled, f4M="
2977	<<f4M<<", fPDG="<<fPD<<", fCg="<<fCg<<", fBN="<<fBN<<G4endl;
2978	#endif
2979	G4QHadron* h2H = new G4QHadron(h2QPDG.GetPDGCode(),h24M);
2980	theResult->push_back(h2H); // (delete equivalent)
2981	#ifdef debug
2982	G4LorentzVector n4M=h2H->Get4Momentum();
2983	G4int nPD=h2H->GetPDGCode();
2984	G4int nCg=h2H->GetCharge();
2985	G4int nBN=h2H->GetBaryonNumber();
2986	G4cout<<"-EMC->>>G4QIonIonCollision::Breed:CorStrHadr Prod-2 is filled, n4M="
2987	<<n4M<<", nPDG="<<nPD<<", nCg="<<nCg<<", nBN="<<nBN<<G4endl;
2988	#endif
2989	#ifdef edebug
2990	G4cout<<"-EMC-...HS-Chipo...G4QIonIonCollision::Breeder:DecCHECK, Ch4M="
2991	<<mi4M<<", ChQC="<<miQC<<", d4M="<<mi4M-h14M-h24M<<G4endl;
2992	#endif
2993	}
2994	}
2995	else
2996	{
2997	G4QHadron* sHad = new G4QHadron(miPDG, mi4M);
2998	theResult->push_back(sHad); // The original string=hadron is filled
2999	#ifdef debug
3000	G4cout<<">>>>>>G4QIonIonCollision::Breeder: CorStr=Hadr is Filled, 4M="
3001	<<curString4M<<", StPDG="<<miPDG<<G4endl;
3002	#endif
3003	}
3004	#ifdef edebug
3005	G4cout<<"-EMC-...Cor...G4QIonIonCollision::Breeder:StHadCor CHECK Sum="<<s4M
3006	<<" =? "<<mi4M+pHadron->Get4Momentum()<<G4endl;
3007	#endif
3008	#ifdef debug
3009	G4cout<<">>>G4QIonIonCollision::Breeder:Corrected String+Hadr="<<miPDG
3010	<<mi4M<<" by Hadron #"<<reha<<G4endl;
3011	#endif
3012	continue; // Continue the LOOP over the curStrings
3013	}
3014	else
3015	{
3016	#ifdef debug
3017	G4cout<<"-EMC->>>G4QIonIonCollision::Breeder: Str+Hadr Failed, 4M="<<curString4M
3018	<<", PDG="<<miPDG<<G4endl;
3019	#endif
3020	}
3021	// @@@ convert string to Quasmon with curString4M
3022	G4QContent curStringQC=curString->GetQC();
3023	G4Quasmon* stringQuasmon = new G4Quasmon(curStringQC, curString4M);
3024	theQuasmons.push_back(stringQuasmon);
3025	continue; // Continue the LOOP over the curStrings
3026	} // End of IF(Can try the String-Hadron correction
3027	} // End of IF(NO_Hadrons) = Problem solution namespace
3028	G4Quasmon tmpQ; // @@ an issue of Q to decay resonances
3029	G4int nHfin=0;
3030	if(theHadrons) nHfin=theHadrons->size();
3031	else // !! Sum Up all strings and convert them in a Quasmon (Exception for development)
3032	{
3033	G4LorentzVector ss4M(0.,0.,0.,0.);
3034	G4QContent ssQC(0,0,0,0,0,0);
3035	G4int tnSt=strings.size();
3036	for(G4int i=astring; i < tnSt; ++i)
3037	{
3038	G4LorentzVector pS4M=strings[i]->Get4Momentum(); // String 4-momentum
3039	ss4M+=pS4M;
3040	G4QContent pSQC=strings[i]->GetQC(); // String Quark Content
3041	ssQC+=pSQC;
3042	#ifdef debug
3043	G4cout<<"====>G4QIonIonCollision::Breed:S#"<<i<<",4M="<<pS4M<<",QC="<<pSQC<<G4endl;
3044	#endif
3045	}
3046	#ifdef debug
3047	G4cout<<"==>G4QIonIonCollision::Breed:AllStrings are summed up in a Quasmon"<<G4endl;
3048	#endif
3049	G4Quasmon* stringQuasmon = new G4Quasmon(ssQC, ss4M);
3050	theQuasmons.push_back(stringQuasmon);
3051	break; // break the LOOP over Strings
3052	}
3053	#ifdef debug
3054	G4cout<<"G4QIonIonCollision::Breeder: Trying to decay hadrons #ofHRes="<<nHfin<<G4endl;
3055	#endif
3056	for(G4int aTrack=0; aTrack<nHfin; aTrack++)
3057	{
3058	G4QHadron* curHadron=(*theHadrons)[aTrack];
3059	G4int hPDG=curHadron->GetPDGCode();
3060	#ifdef edebug
3061	G4LorentzVector curH4M=curHadron->Get4Momentum();
3062	G4int curHCh=curHadron->GetCharge();
3063	G4int curHBN=curHadron->GetBaryonNumber();
3064	#endif
3065	#ifdef debug
3066	G4cout<<">>>>>>>>G4QIonIonCollision::Breeder:S#"<<astring<<",H#"<<aTrack<<",PDG="
3067	<<hPDG<<",4M="<<curHadron->Get4Momentum()<<G4endl;
3068	#endif
3069	if(std::abs(hPDG)%10 > 2)
3070	{
3071	G4QHadronVector* tmpQHadVec=tmpQ.DecayQHadron(curHadron); // It deletes curHadron
3072	#ifdef debug
3073	G4cout<<"G4QIonIonCollision::Breed:-DECAY'S DONE-,nH="<<tmpQHadVec->size()<<G4endl;
3074	#endif
3075	//G4int tmpS=tmpQHadVec->size(); // "The elegant method" (tested) is commented
3076	//theResult->resize(tmpS+theResult->size()); // Resize theQHadrons length
3077	//copy(tmpQHadVec->begin(), tmpQHadVec->end(), theResult->end()-tmpS);
3078	for(unsigned aH=0; aH < tmpQHadVec->size(); aH++)
3079	{
3080	theResult->push_back((*tmpQHadVec)[aH]);// TheDecayProduct of TheHadron is filled
3081	#ifdef edebug
3082	G4QHadron* prodH =(*tmpQHadVec)[aH];
3083	G4LorentzVector p4M=prodH->Get4Momentum();
3084	G4int PDG=prodH->GetPDGCode();
3085	G4int Chg=prodH->GetCharge();
3086	G4int BaN=prodH->GetBaryonNumber();
3087	curH4M-=p4M;
3088	curString4M-=p4M;
3089	curStrChg-=Chg;
3090	curStrBaN-=BaN;
3091	curHCh-=Chg;
3092	curHBN-=BaN;
3093	G4cout<<"-EMC->>>>G4QIonIonCollision::Breed:Str*Filled, 4M="<<p4M<<", PDG="<<PDG
3094	<<", Chg="<<Chg<<", BaN="<<BaN<<G4endl;
3095	#endif
3096	}
3097	#ifdef edebug
3098	G4cout<<"-EMC-.G4QIn::Br:Dec,r4M="<<curH4M<<",rC="<<curHCh<<",rB="<<curHBN<<G4endl;
3099	#endif
3100	tmpQHadVec->clear();
3101	delete tmpQHadVec; // Who calls DecayQHadron is responsible for clear & delete
3102	}
3103	else // Chipolino is not checked here
3104	{
3105	theResult->push_back(curHadron); // The original hadron is filled
3106	#ifdef edebug
3107	curString4M-=curH4M;
3108	G4int curCh=curHadron->GetCharge();
3109	G4int curBN=curHadron->GetBaryonNumber();
3110	curStrChg-=curCh;
3111	curStrBaN-=curBN;
3112	G4cout<<"-EMC->>>>>>G4QIonIonCollision::Breeder: curH filled 4M="<<curH4M<<",PDG="
3113	<<curHadron->GetPDGCode()<<", Chg="<<curCh<<", BaN="<<curBN<<G4endl;
3114	#endif
3115	}
3116	}
3117	// clean up (the issues are filled to theResult)
3118	if(theHadrons) delete theHadrons;
3119	#ifdef edebug
3120	G4cout<<"-EMC-.......G4QIonIonCollision::Breeder: StringDecay CHECK, r4M="<<curString4M
3121	<<", rChg="<<curStrChg<<", rBaN="<<curStrBaN<<G4endl;
3122	#endif
3123	} // End of the main LOOP over decaying strings
3124	G4LorentzVector rp4M=theProjNucleus.Get4Momentum(); // projNucleus 4-momentum in LS
3125	G4int rpPDG=theProjNucleus.GetPDG();
3126	G4QHadron* resPNuc = new G4QHadron(rpPDG,rp4M);
3127	theResult->push_back(resPNuc); // Fill the residual nucleus
3128	G4LorentzVector rt4M=theTargNucleus.Get4Momentum(); // Nucleus 4-momentum in LS
3129	G4int rtPDG=theTargNucleus.GetPDG();
3130	G4QHadron* resTNuc = new G4QHadron(rtPDG,rt4M);
3131	theResult->push_back(resTNuc); // Fill the residual nucleus
3132	#ifdef edebug
3133	G4LorentzVector s4M(0.,0.,0.,0.); // Sum of the Result in LS
3134	G4int rCh=totChg;
3135	G4int rBN=totBaN;
3136	G4int nHadr=theResult->size();
3137	G4int nQuasm=theQuasmons.size();
3138	G4cout<<"-EMCLS-G4QInel::Breeder: #ofHadr="<<nHadr<<", #OfQ="<<nQuasm<<", rPA="<<rp4M.m()
3139	<<"="<<G4QNucleus(rpPDG).GetGSMass()<<", rTA="<<rt4M.m()<<"="
3140	<<G4QNucleus(rtPDG).GetGSMass()<<G4endl;
3141	for(G4int i=0; i<nHadr; i++)
3142	{
3143	G4LorentzVector hI4M=(*theResult)[i]->Get4Momentum();
3144	s4M+=hI4M;
3145	G4int hChg=(*theResult)[i]->GetCharge();
3146	rCh-=hChg;
3147	G4int hBaN=(*theResult)[i]->GetBaryonNumber();
3148	rBN-=hBaN;
3149	G4cout<<"-EMCLS-G4QIonIonCollision::Breeder: Hadron#"<<i<<", 4M="<<hI4M<<", PDG="
3150	<<(*theResult)[i]->GetPDGCode()<<", C="<<hChg<<", B="<<hBaN<<G4endl;
3151	}
3152	for(G4int i=0; i<nQuasm; i++)
3153	{
3154	G4LorentzVector hI4M=theQuasmons[i]->Get4Momentum();
3155	s4M+=hI4M;
3156	G4int hChg=theQuasmons[i]->GetCharge();
3157	rCh-=hChg;
3158	G4int hBaN=theQuasmons[i]->GetBaryonNumber();
3159	rBN-=hBaN;
3160	G4cout<<"-EMCLS-G4QIonIonCollision::Breeder: Quasmon#"<<i<<", 4M="<<hI4M<<", C="<<hChg
3161	<<", B="<<hBaN<<G4endl;
3162	}
3163	G4cout<<"-EMCLS-G4QInel::Breed: LS r4M="<<s4M-totLS4M<<",rC="<<rCh<<",rB="<<rBN<<G4endl;
3164	#endif
3165	// Now we need to coolect particles for creation of a Quasmon @@ improve !!
3166	G4int nRes=theResult->size();
3167	if(nRes>2)
3168	{
3169	G4QHadronVector::iterator ih;
3170	G4QHadronVector::iterator nih;
3171	G4QHadronVector::iterator mih;
3172	G4QHadronVector::iterator lst=theResult->end();
3173	lst--;
3174	G4double minMesEn=DBL_MAX;
3175	G4double minBarEn=DBL_MAX;
3176	G4bool nfound=false;
3177	G4bool mfound=false;
3178	for(ih = theResult->begin(); ih < theResult->end(); ++ih) if(ih != lst)
3179	{
3180	G4LorentzVector h4M=(*ih)->Get4Momentum();
3181	G4int hPDG=(*ih)->GetPDGCode();
3182	#ifdef debug
3183	G4cout<<"%->G4QIonIonCollision::Breeder: TRY hPDG="<<hPDG<<", h4M="<<h4M<<G4endl;
3184	#endif
3185	if(hPDG>1111 && hPDG<3333)
3186	{
3187	G4double bE=h4M.e()-(*ih)->GetMass();
3188	if(bE < minBarEn)
3189	{
3190	minBarEn=bE;
3191	nih=ih;
3192	nfound=true;
3193	}
3194	}
3195	else if(hPDG>-1111)
3196	{
3197	G4double mE=h4M.e();
3198	if(mE < minMesEn)
3199	{
3200	minMesEn=mE;
3201	mih=ih;
3202	mfound=true;
3203	}
3204	}
3205	}
3206	if(nfound && mfound && minMesEn+minBarEn < maxEn)
3207	{
3208	G4QHadron* resNuc = (*theResult)[nRes-1]; // ResidualNucleus is theLastH
3209	theResult->pop_back(); // Fill the QHad-nucleus later
3210	G4LorentzVector q4M=(nih)->Get4Momentum()+(mih)->Get4Momentum();
3211	G4QContent qQC=(nih)->GetQC()+(mih)->GetQC();
3212	maxEn -= minBarEn+minMesEn; // Reduce the energy limit
3213	#ifdef debug
3214	G4cout<<"%->G4QIonIonCollision::Breeder:Exclude,4M="<<(*nih)->Get4Momentum()
3215	<<", & 4m="<<(*mih)->Get4Momentum()<<G4endl;
3216	#endif
3217	delete (*nih);
3218	delete (*mih);
3219	if(nih > mih) // Exclude nucleon first
3220	{
3221	theResult->erase(nih); // erase Baryon from theResult
3222	theResult->erase(mih); // erase Meson from theResult
3223	}
3224	else // Exclude meson first
3225	{
3226	theResult->erase(mih); // erase Baryon from theResult
3227	theResult->erase(nih); // erase Meson from theResult
3228	}
3229	#ifdef debug
3230	G4cout<<"%->G4QI::Breed: BeforeLOOP, dE="<<maxEn<<", nR="<<theResult->size()<<G4endl;
3231	#endif
3232	if(maxEn > 0.) // More hadrons to absorbe
3233	{
3234	for(ih = theResult->begin(); ih < theResult->end(); ih++)
3235	{
3236	G4LorentzVector h4M=(*ih)->Get4Momentum();
3237	G4int hPDG=(*ih)->GetPDGCode();
3238	G4double dE=0.;
3239	if (hPDG> 1111 && hPDG<3333) dE=h4M.e()-(*ih)->GetMass(); // Baryons
3240	else if(hPDG>-1111) dE=h4M.e(); // Mesons
3241	else continue; // Don't consider anti-baryons
3242	if(dE < maxEn)
3243	{
3244	maxEn-=dE;
3245	q4M+=h4M;
3246	qQC+=(*ih)->GetQC();
3247	#ifdef debug
3248	G4cout<<"%->G4QIonIonCollision::Breed:Exclude,4M="<<h4M<<",dE="<<maxEn<<G4endl;
3249	#endif
3250	delete (*ih);
3251	theResult->erase(ih); // erase Hadron from theResult
3252	--ih;
3253	}
3254	}
3255	}
3256	G4Quasmon* softQuasmon = new G4Quasmon(qQC, q4M); // SoftPart -> Quasmon
3257	#ifdef debug
3258	G4cout<<"%->G4QIonIonCollision::Breed:QuasmonIsFilled,4M="<<q4M<<",QC="<<qQC<<G4endl;
3259	#endif
3260	theQuasmons.push_back(softQuasmon);
3261	theResult->push_back(resNuc);
3262	}
3263	#ifdef edebug
3264	G4LorentzVector f4M(0.,0.,0.,0.); // Sum of the Result in LS
3265	G4int fCh=totChg;
3266	G4int fBN=totBaN;
3267	G4int nHd=theResult->size();
3268	G4int nQm=theQuasmons.size();
3269	G4cout<<"-EMCLS-G4QIonIonCollision::Breeder:#ofHadr="<<nHd<<", #OfQ="<<nQm
3270	<<",rPA="<<rt4M.m()<<"="<<G4QNucleus(rtPDG).GetGSMass()
3271	<<",rTA="<<rt4M.m()<<"="<<G4QNucleus(rtPDG).GetGSMass()<<G4endl;
3272	for(G4int i=0; i<nHd; i++)
3273	{
3274	G4LorentzVector hI4M=(*theResult)[i]->Get4Momentum();
3275	f4M+=hI4M;
3276	G4int hChg=(*theResult)[i]->GetCharge();
3277	fCh-=hChg;
3278	G4int hBaN=(*theResult)[i]->GetBaryonNumber();
3279	fBN-=hBaN;
3280	G4cout<<"-EMCLS-G4QIonIonCollision::Breeder: Hadron#"<<i<<", 4M="<<hI4M<<", PDG="
3281	<<(*theResult)[i]->GetPDGCode()<<", C="<<hChg<<", B="<<hBaN<<G4endl;
3282	}
3283	for(G4int i=0; i<nQm; i++)
3284	{
3285	G4LorentzVector hI4M=theQuasmons[i]->Get4Momentum();
3286	f4M+=hI4M;
3287	G4int hChg=theQuasmons[i]->GetCharge();
3288	fCh-=hChg;
3289	G4int hBaN=theQuasmons[i]->GetBaryonNumber();
3290	fBN-=hBaN;
3291	G4cout<<"-EMCLS-G4QIonIonCollision::Breeder: Quasmon#"<<i<<", 4M="<<hI4M<<", C="
3292	<<hChg<<", B="<<hBaN<<G4endl;
3293	}
3294	G4cout<<"-EMCLS-G4QInel::Breed:, r4M="<<f4M-totLS4M<<", rC="<<fCh<<",rB="<<fBN<<G4endl;
3295	#endif
3296	} // End of the soft Quasmon Creation
3297	return;
3298	} // End of Breeder
3299
3300	// Excite double diffractive string
3301	G4bool G4QIonIonCollision::ExciteDiffParticipants(G4QHadron* projectile,
3302	G4QHadron* target) const
3303	{
3304	G4LorentzVector Pprojectile=projectile->Get4Momentum();
3305	G4double Mprojectile=projectile->GetMass();
3306	G4double Mprojectile2=Mprojectile*Mprojectile;
3307	G4LorentzVector Ptarget=target->Get4Momentum();
3308	G4double Mtarget=target->GetMass();
3309	G4double Mtarget2=Mtarget*Mtarget;
3310	#ifdef debug
3311	G4cout<<"G4QInel::ExciteDiffPartici: Ep="<<Pprojectile.e()<<", Et="<<Ptarget.e()<<G4endl;
3312	#endif
3313	// Transform momenta to cms and then rotate parallel to z axis;
3314	G4LorentzVector Psum=Pprojectile+Ptarget;
3315	G4LorentzRotation toCms(-Psum.boostVector()); // Boost Rotation to CMS
3316	G4LorentzVector Ptmp=toCms*Pprojectile;
3317	if(Ptmp.pz()<=0.) // "String" moving backwards in CMS, abort collision !! ? M.K.
3318	{
3319	#ifdef debug
3320	G4cout<<"G4QIonIonCollision::ExciteDiffParticipants:1 abort Collision!! 1"<<G4endl;
3321	#endif
3322	return false;
3323	}
3324	toCms.rotateZ(-Ptmp.phi());
3325	toCms.rotateY(-Ptmp.theta());
3326	#ifdef debug
3327	G4cout<<"G4QIonIonCollision::ExciteDiffParticipantts:BeforBoost Pproj="<<Pprojectile
3328	<<",Ptarg="<<Ptarget<<G4endl;
3329	#endif
3330	G4LorentzRotation toLab(toCms.inverse()); // Boost Rotation to LabSys (LS)
3331	Pprojectile.transform(toCms);
3332	Ptarget.transform(toCms);
3333	#ifdef debug
3334	G4cout<< "G4QInelast::ExciteDiffParticipantts: AfterBoost Pproj="<<Pprojectile<<",Ptarg="
3335	<<Ptarget<<", cms4M="<<Pprojectile+Ptarget<<G4endl;
3336	G4cout<<"G4QIonIonCollis::ExciteDiffParticipants:ProjX+="<<Pprojectile.plus()<<",ProjX-="
3337	<<Pprojectile.minus()<<", tX+="<< Ptarget.plus()<<",tX-="<<Ptarget.minus()<<G4endl;
3338	#endif
3339	G4LorentzVector Qmomentum(0.,0.,0.,0.);
3340	G4int whilecount=0;
3341	#ifdef debug
3342	G4cout<<"G4QIonIonCollision::ExciteDiffParticipants: Before DO"<<G4endl;
3343	#endif
3344	do
3345	{
3346	// Generate pt
3347	G4double maxPtSquare=sqr(Ptarget.pz());
3348	#ifdef debug
3349	G4cout<<"G4QIonIonCollision::ExciteDiffParticipants: maxPtSq="<<maxPtSquare<<G4endl;
3350	if(whilecount++>=500 && whilecount%100==0) // @@ M.K. Hardwired limits
3351	G4cout<<"G4QIonIonCollision::ExciteDiffParticipantts: can loop, loopCount="<<whilecount
3352	<<", maxPtSquare="<<maxPtSquare<<G4endl;
3353	#endif
3354	if(whilecount>1000) // @@ M.K. Hardwired limits
3355	{
3356	#ifdef debug
3357	G4cout<<"G4QIonIonCollision::ExciteDiffParticipants: 2 abort Loop!! 2"<<G4endl;
3358	#endif
3359	return false; // Ignore this interaction
3360	}
3361	Qmomentum=G4LorentzVector(GaussianPt(widthOfPtSquare,maxPtSquare),0);
3362	#ifdef debug
3363	G4cout<<"G4QIonIonCollision::ExciteDiffParticipants: generated Pt="<<Qmomentum
3364	<<", ProjPt="<<Pprojectile+Qmomentum<<", TargPt="<<Ptarget-Qmomentum<<G4endl;
3365	#endif
3366	// Momentum transfer
3367	G4double Xmin=0.;
3368	G4double Xmax=1.;
3369	G4double Xplus =ChooseX(Xmin,Xmax);
3370	G4double Xminus=ChooseX(Xmin,Xmax);
3371	#ifdef debug
3372	G4cout<<"G4QIonIonCollision::ExciteDiffParticipant:X+="<<Xplus<<",X-="<<Xminus<<G4endl;
3373	#endif
3374	G4double pt2=Qmomentum.vect().mag2();
3375	G4double Qplus =-pt2/Xminus/Ptarget.minus();
3376	G4double Qminus= pt2/Xplus /Pprojectile.plus();
3377	Qmomentum.setPz((Qplus-Qminus)/2);
3378	Qmomentum.setE( (Qplus+Qminus)/2);
3379	#ifdef debug
3380	G4cout<<"G4QInelast::ExciteDiffParticip: Qplus="<<Qplus<<", Qminus="<<Qminus<<", pt2="
3381	<<pt2<<", Qmomentum="<<Qmomentum<<", ProjM="<<(Pprojectile+Qmomentum).mag()
3382	<<", TargM="<<(Ptarget-Qmomentum).mag()<<G4endl;
3383	#endif
3384	} while((Pprojectile+Qmomentum).mag2()<=Mprojectile2 \|\|
3385	(Ptarget-Qmomentum).mag2()<=Mtarget2);
3386	Pprojectile += Qmomentum;
3387	Ptarget -= Qmomentum;
3388	#ifdef debug
3389	G4cout<<"G4QInelast::ExciteDiffParticipan: Proj(Q)="<<Pprojectile<<", Targ(Q)="<<Ptarget
3390	<<", Proj(back)="<<toLabPprojectile<<", Targ(bac)="<< toLabPtarget << G4endl;
3391	#endif
3392	// Transform back and update SplitableHadron Participant.
3393	Pprojectile.transform(toLab);
3394	Ptarget.transform(toLab);
3395	#ifdef debug
3396	G4cout<< "G4QIonIonCollision::ExciteDiffParticipants:TargetMass="<<Ptarget.mag()<<G4endl;
3397	#endif
3398	target->Set4Momentum(Ptarget);
3399	#ifdef debug
3400	G4cout<<"G4QIonIonCollision::ExciteDiffParticipant:ProjMass="<<Pprojectile.mag()<<G4endl;
3401	#endif
3402	projectile->Set4Momentum(Pprojectile);
3403	return true;
3404	} // End of ExciteDiffParticipants
3405
3406
3407	// Excite single diffractive string
3408	G4bool G4QIonIonCollision::ExciteSingDiffParticipants(G4QHadron* projectile,
3409	G4QHadron* target) const
3410	{
3411	G4LorentzVector Pprojectile=projectile->Get4Momentum();
3412	G4double Mprojectile=projectile->GetMass();
3413	G4double Mprojectile2=Mprojectile*Mprojectile;
3414	G4LorentzVector Ptarget=target->Get4Momentum();
3415	G4double Mtarget=target->GetMass();
3416	G4double Mtarget2=Mtarget*Mtarget;
3417	#ifdef debug
3418	G4cout<<"G4QInel::ExSingDiffPartici: Ep="<<Pprojectile.e()<<", Et="<<Ptarget.e()<<G4endl;
3419	#endif
3420	G4bool KeepProjectile= G4UniformRand() > 0.5;
3421	// Reset minMass of the non diffractive particle to its value, (minus for rounding...)
3422	if(KeepProjectile )
3423	{
3424	#ifdef debug
3425	G4cout<<"-1/2-G4QIonIonCollision::ExSingDiffParticipants: Projectile is fixed"<<G4endl;
3426	#endif
3427	Mprojectile2 = projectile->GetMass2()*(1.-perCent); // Isn't it too big reduction? M.K.
3428	}
3429	else
3430	{
3431	#ifdef debug
3432	G4cout<<"---1/2---G4QIonIonCollision::ExSingDiffParticipants: Target is fixed"<<G4endl;
3433	#endif
3434	Mtarget2 = target->GetMass2()*(1.-perCent); // @@ Isn't it too big reduction? M.K.
3435	}
3436	// @@ From this point it repeats the Diffractional excitation (? Use flag ?)
3437	// Transform momenta to cms and then rotate parallel to z axis;
3438	G4LorentzVector Psum=Pprojectile+Ptarget;
3439	G4LorentzRotation toCms(-Psum.boostVector()); // Boost Rotation to CMS
3440	G4LorentzVector Ptmp=toCms*Pprojectile;
3441	if(Ptmp.pz()<=0.) // "String" moving backwards in CMS, abort collision !! ? M.K.
3442	{
3443	#ifdef debug
3444	G4cout<<"G4QIonIonCollision::ExciteSingDiffParticipants: 1abortCollision1"<<G4endl;
3445	#endif
3446	return false;
3447	}
3448	toCms.rotateZ(-Ptmp.phi());
3449	toCms.rotateY(-Ptmp.theta());
3450	#ifdef debug
3451	G4cout<<"G4QInel::ExciteSingDiffParticipantts: Be4Boost Pproj="<<Pprojectile<<", Ptarg="
3452	<<Ptarget<<G4endl;
3453	#endif
3454	G4LorentzRotation toLab(toCms.inverse()); // Boost Rotation to LabSys (LS)
3455	Pprojectile.transform(toCms);
3456	Ptarget.transform(toCms);
3457	#ifdef debug
3458	G4cout<< "G4QInelast::ExciteDiffParticipantts: AfterBoost Pproj="<<Pprojectile<<"Ptarg="
3459	<<Ptarget<<", cms4M="<<Pprojectile+Ptarget<<G4endl;
3460
3461	G4cout<<"G4QInelast::ExciteDiffParticipantts: ProjX+="<<Pprojectile.plus()<<", ProjX-="
3462	<<Pprojectile.minus()<<", TargX+="<< Ptarget.plus()<<", TargX-="<<Ptarget.minus()
3463	<<G4endl;
3464	#endif
3465	G4LorentzVector Qmomentum(0.,0.,0.,0.);
3466	G4int whilecount=0;
3467	do
3468	{
3469	// Generate pt
3470	G4double maxPtSquare=sqr(Ptarget.pz());
3471	if(whilecount++>=500 && whilecount%100==0) // @@ M.K. Hardwired limits
3472	#ifdef debug
3473	G4cout<<"G4QIonIonCollision::ExciteSingDiffParticipantts: can loop, loopCount="
3474	<<whilecount<<", maxPtSquare="<<maxPtSquare<<G4endl;
3475	#endif
3476	if(whilecount>1000) // @@ M.K. Hardwired limits
3477	{
3478	#ifdef debug
3479	G4cout<<"G4QIonIonCollision::ExciteSingDiffParticipants:2 abortLoop!! 2"<<G4endl;
3480	#endif
3481	return false; // Ignore this interaction
3482	}
3483	Qmomentum=G4LorentzVector(GaussianPt(widthOfPtSquare,maxPtSquare),0);
3484	#ifdef debug
3485	G4cout<<"G4QInelast::ExciteSingDiffParticipants: generated Pt="<<Qmomentum
3486	<<", ProjPt="<<Pprojectile+Qmomentum<<", TargPt="<<Ptarget-Qmomentum<<G4endl;
3487	#endif
3488	// Momentum transfer
3489	G4double Xmin=0.;
3490	G4double Xmax=1.;
3491	G4double Xplus =ChooseX(Xmin,Xmax);
3492	G4double Xminus=ChooseX(Xmin,Xmax);
3493	#ifdef debug
3494	G4cout<<"G4QInel::ExciteSingDiffPartici: X-plus="<<Xplus<<", X-minus="<<Xminus<<G4endl;
3495	#endif
3496	G4double pt2=G4ThreeVector(Qmomentum.vect()).mag2();
3497	G4double Qplus =-pt2/Xminus/Ptarget.minus();
3498	G4double Qminus= pt2/Xplus /Pprojectile.plus();
3499	if (KeepProjectile)
3500	Qminus=(projectile->GetMass2()+pt2)/(Pprojectile.plus()+Qplus) - Pprojectile.minus();
3501	else Qplus=Ptarget.plus() - (target->GetMass2()+pt2)/(Ptarget.minus()-Qminus);
3502	Qmomentum.setPz((Qplus-Qminus)/2);
3503	Qmomentum.setE( (Qplus+Qminus)/2);
3504	#ifdef debug
3505	G4cout<<"G4QInel::ExciteDiffParticip: Qplus="<<Qplus<<", Qminus="<<Qminus<<", pt2="
3506	<<pt2<<", Qmomentum="<<Qmomentum<<", ProjM="<<(Pprojectile+Qmomentum).mag()
3507	<<", TargM="<<(Ptarget-Qmomentum).mag()<<G4endl;
3508	#endif
3509	// while is different from the Double Diffractive Excitation (@@ !)
3510	//} while((Pprojectile+Qmomentum).mag2()<= Mprojectile2 \|\|
3511	// (Ptarget-Qmomentum).mag2()<=Mtarget2);
3512	} while((Ptarget-Qmomentum).mag2()<=Mtarget2 \|\|
3513	(Pprojectile+Qmomentum).mag2()<=Mprojectile2 \|\|
3514	(Ptarget-Qmomentum).e() < 0. \|\| (Pprojectile+Qmomentum).e() < 0.);
3515	Pprojectile += Qmomentum;
3516	Ptarget -= Qmomentum;
3517	#ifdef debug
3518	G4cout<<"G4QIonIonCollision::ExciteSingDiffParticipan: Proj(Q)="<<Pprojectile<<"(E="
3519	<<Pprojectile.e()<<"), Targ(Q)="<<Ptarget<<"(E="<<Ptarget.e()
3520	<<"), Proj(back)="<<toLabPprojectile<<", Targ(bac)="<< toLabPtarget << G4endl;
3521	#endif
3522	// Transform back and update SplitableHadron Participant.
3523	Pprojectile.transform(toLab);
3524	Ptarget.transform(toLab);
3525	#ifdef debug
3526	G4cout<< "G4QIonIonCollision::ExciteSingleDiffParticipants: TgM="<<Ptarget.mag()<<G4endl;
3527	#endif
3528	target->Set4Momentum(Ptarget);
3529	#ifdef debug
3530	G4cout<<"G4QIonIonCollision::ExciteSingleParticipants:ProjM="<<Pprojectile.mag()<<G4endl;
3531	#endif
3532	projectile->Set4Momentum(Pprojectile);
3533	return true;
3534	} // End of ExciteSingleDiffParticipants
3535
3536	void G4QIonIonCollision::SetParameters(G4int nC,G4double stT, G4double tbD, G4double SigPt)
3537	{// ======================================================================================
3538	nCutMax = nC; // max number of pomeron cuts
3539	stringTension = stT; // string tension for absorbed energy
3540	tubeDensity = tbD; // Flux Tube Density of nuclear nucleons
3541	widthOfPtSquare = -2SigPtSigPt; // width^2 of pt for string excitation
3542	}
3543
3544	G4double G4QIonIonCollision::ChooseX(G4double Xmin, G4double Xmax) const
3545	{
3546	// choose an x between Xmin and Xmax with P(x) ~ 1/x @@ M.K. -> 1/sqrt(x)
3547	//G4double range=Xmax-Xmin;
3548	if(Xmax == Xmin) return Xmin;
3549	if( Xmin < 0. \|\| Xmax < Xmin)
3550	{
3551	G4cerr<<"***G4QIonIonCollision::ChooseX: Xmin="<<Xmin<<", Xmax="<<Xmax<< G4endl;
3552	G4Exception("G4QIonIonCollision::ChooseX:","72",FatalException,"Bad X or X-Range");
3553	}
3554	//G4double x;
3555	//do {x=Xmin+G4UniformRand()*range;} while ( Xmin/x < G4UniformRand() );
3556	G4double sxi=std::sqrt(Xmin);
3557	G4double x=sqr(sxi+G4UniformRand()*(std::sqrt(Xmax)-sxi));
3558	#ifdef debug
3559	G4cout<<"G4QIonIonCollision::ChooseX: DiffractiveX="<<x<<G4endl;
3560	#endif
3561	return x;
3562	} // End of ChooseX
3563
3564	// Pt distribution @@ one can use 1/(1+A*Pt^2)^B
3565	G4ThreeVector G4QIonIonCollision::GaussianPt(G4double widthSq, G4double maxPtSquare) const
3566	{
3567	#ifdef debug
3568	G4cout<<"G4QIonIonCollision::GaussianPt: w^2="<<widthSq<<",maxPt2="<<maxPtSquare<<G4endl;
3569	#endif
3570	G4double pt2=0.;
3571	G4double rm=maxPtSquare/widthSq; // Negative
3572	if(rm>-.01) pt2=widthSq(std::sqrt(1.-G4UniformRand()(1.-sqr(1.+rm)))-1.);
3573	else pt2=widthSqstd::log(1.-G4UniformRand()(1.-std::exp(rm)));
3574	pt2=std::sqrt(pt2); // It is not pt2, but pt now
3575	G4double phi=G4UniformRand()*twopi;
3576	return G4ThreeVector(pt2std::cos(phi),pt2std::sin(phi),0.);
3577	} // End of GaussianPt
3578
3579	G4int G4QIonIonCollision::SumPartonPDG(G4int PDG1, G4int PDG2) const
3580	{
3581	if (PDG1 < 7 && PDG1 > 0 && PDG2 < 7 && PDG2 > 0) // Sum up two Q in DiQ (S=0)
3582	{
3583	if(PDG1 > PDG2) return PDG11000+PDG2100+1;
3584	else return PDG21000+PDG1100+1;
3585	}
3586	else if (PDG1 >-7 && PDG1 < 0 && PDG2 >-7 && PDG2 < 0) // Sum up two AQ in ADiQ (S=0)
3587	{
3588	if(-PDG1 > -PDG2) return PDG11000+PDG2100-1;
3589	else return PDG21000+PDG1100-1;
3590	}
3591	else if (PDG1 <-99 && PDG2 < 7 && PDG2 > 0) // Sum up Q and ADiQ in AQ
3592	{
3593	G4int PDG=-PDG1/100;
3594	if(PDG2==PDG/10) return -PDG%10;
3595	if(PDG2==PDG%10) return -PDG/10;
3596	else
3597	{
3598	G4cerr<<"***4QFragmentation::SumPartonPDG: PDG1="<<PDG1<<", PDG2="<<PDG2<<G4endl;
3599	G4Exception("G4QIonIonCollision::SumPartonPDG:","72",FatalException,"Q&ADiQnoMatch");
3600	}
3601	}
3602	else if (PDG2 <-99 && PDG1 < 7 && PDG1 > 0) // Sum up ADiQ and Q in AQ
3603	{
3604	G4int PDG=-PDG2/100;
3605	if(PDG1==PDG/10) return -PDG%10;
3606	if(PDG1==PDG%10) return -PDG/10;
3607	else
3608	{
3609	G4cerr<<"***4QFragmentation::SumPartonPDG: PDG1="<<PDG1<<", PDG2="<<PDG2<<G4endl;
3610	G4Exception("G4QIonIonCollision::SumPartonPDG:","72",FatalException,"ADiQ&QnoMatch");
3611	}
3612	}
3613	else if (PDG1 > 99 && PDG2 >-7 && PDG2 < 0) // Sum up DiQ and AQ in Q
3614	{
3615	G4int PDG=PDG1/100;
3616	if(PDG2==-PDG/10) return PDG%10;
3617	if(PDG2==-PDG%10) return PDG/10;
3618	else
3619	{
3620	G4cerr<<"***4QFragmentation::SumPartonPDG: PDG1="<<PDG1<<", PDG2="<<PDG2<<G4endl;
3621	G4Exception("G4QIonIonCollision::SumPartonPDG:","72",FatalException,"DiQ&AQnoMatch");
3622	}
3623	}
3624	else if (PDG2 > 99 && PDG1 >-7 && PDG1 < 0) // Sum up AQ and DiQ in Q
3625	{
3626	G4int PDG=PDG2/100;
3627	if(PDG1==-PDG/10) return PDG%10;
3628	if(PDG1==-PDG%10) return PDG/10;
3629	else
3630	{
3631	G4cerr<<"***G4QIonIonCollision::SumPartonPDG: PDG1="<<PDG1<<", PDG2="<<PDG2<<G4endl;
3632	G4Exception("G4QIonIonCollision::SumPartonPDG:","72",FatalException,"AQ&DiQnoMatch");
3633	}
3634	}
3635	else
3636	{
3637	G4cerr<<"***G4QIonIonCollision::SumPartonPDG: PDG1="<<PDG1<<", PDG2="<<PDG2<<G4endl;
3638	G4Exception("G4QIonIonCollision::SumPartonPDG:","72",FatalException,"Can'tSumUpParts");
3639	}
3640	return 0;
3641	} // End of SumPartonPDG
3642
3643	// Reduces quark pairs (unsigned 2 digits) to quark singles (unsigned)
3644	std::pair<G4int,G4int> G4QIonIonCollision::ReducePair(G4int P1, G4int P2) const
3645	{
3646	#ifdef debug
3647	G4cout<<"G4QIonIonCollision::ReducePair: Called P1="<<P1<<", P2="<<P2<<G4endl;
3648	#endif
3649	G4int P11=P1/10;
3650	G4int P12=P1%10;
3651	G4int P21=P2/10;
3652	G4int P22=P2%10;
3653	if (P11==P21) return std::make_pair(P12,P22);
3654	else if(P11==P22) return std::make_pair(P12,P21);
3655	else if(P12==P21) return std::make_pair(P11,P22);
3656	else if(P12==P22) return std::make_pair(P11,P21);
3657	//#ifdef debug
3658	G4cout<<"-Warning-G4QIonIonCollision::ReducePair:Failed,P1="<<P1<<",P2="<<P2<<G4endl;
3659	//#endif
3660	return std::make_pair(0,0); // Reduction failed
3661	}
3662
3663	// Select LL/RR (1) or LR/RL (-1) annihilation order (0, if the annihilation is impossible)
3664	G4int G4QIonIonCollision::AnnihilationOrder(G4int LS, G4int MPS, G4int uPDG, G4int mPDG,
3665	G4int sPDG, G4int nPDG) // ^L^^ Partner^^R^
3666	{// ^L^^ Curent ^^R^
3667	G4int Ord=0;
3668	// Curent Partner
3669	if (LS==2 && MPS==2 ) // Fuse 2 Q-aQ strings to DiQ-aDiQ
3670	{
3671	#ifdef debug
3672	G4cout<<"G4QIonIonCollision::AnnihilationOrder:QaQ/QaQ->DiQ-aDiQ, uPDG="<<uPDG<<G4endl;
3673	#endif
3674	if ( (uPDG>0 && sPDG>0 && mPDG<0 && nPDG<0) \|\|
3675	(uPDG<0 && sPDG<0 && mPDG>0 && nPDG>0) ) Ord= 1; // LL/RR
3676	else if( (uPDG>0 && nPDG>0 && mPDG<0 && sPDG<0) \|\|
3677	(uPDG<0 && nPDG<0 && mPDG>0 && sPDG>0) ) Ord=-1; // LR/RL
3678	else G4cerr<<"-Warning-G4QIonIonCollision::AnnihilationOrder: Wrong 22 fusion, L="
3679	<<uPDG<<",R="<<mPDG<<",cL="<<sPDG<<",cR="<<nPDG<<G4endl;
3680	}
3681	else if ( LS==2 && MPS==3 )
3682	{
3683	if (uPDG > 7) // Fuse pLDiQ
3684	{
3685	#ifdef debug
3686	G4cout<<"G4QIonIonCollision::AnnihOrder:pLDiQ, sPDG="<<sPDG<<", nPDG="<<nPDG<<G4endl;
3687	#endif
3688	if ( sPDG<0 && (-sPDG==uPDG/1000 \|\| -sPDG==(uPDG/100)%10) ) Ord= 1; // LL/RR
3689	else if( nPDG<0 && (-nPDG==uPDG/1000 \|\| -nPDG==(uPDG/100)%10) ) Ord=-1; // LR/RL
3690	else G4cerr<<"-Warning-G4QIonIonCollision::AnnihilationOrder: Wrong pLDiQ, L="<<uPDG
3691	<<",R="<<mPDG<<",cL="<<sPDG<<",cR="<<nPDG<<G4endl;
3692	}
3693	else if (mPDG > 7) // Fuse pRDiQ
3694	{
3695	#ifdef debug
3696	G4cout<<"G4QIonIonCollision::AnnihOrder:pRDiQ, sPDG="<<sPDG<<", nPDG="<<nPDG<<G4endl;
3697	#endif
3698	if ( sPDG<0 && (-sPDG==mPDG/1000 \|\| -sPDG==(mPDG/100)%10) ) Ord=-1; // LR/RL
3699	else if( nPDG<0 && (-nPDG==mPDG/1000 \|\| -nPDG==(mPDG/100)%10) ) Ord= 1; // LL/RR
3700	else G4cerr<<"-Warning-G4QIonIonCollision::AnnihilationOrder: Wrong pRDiQ, L="<<uPDG
3701	<<",R="<<mPDG<<",cL="<<sPDG<<",cR="<<nPDG<<G4endl;
3702	}
3703	else if (uPDG <-7) // Fuse pLaDiQ
3704	{
3705	#ifdef debug
3706	G4cout<<"G4QIonIonCollision::AnnihOrder:pLaDiQ, sPDG="<<sPDG<<",nPDG="<<nPDG<<G4endl;
3707	#endif
3708	if ( sPDG>0 && (sPDG==(-uPDG)/1000 \|\| sPDG==((-uPDG)/100)%10) ) Ord= 1; // LL/RR
3709	else if( nPDG>0 && (nPDG==(-uPDG)/1000 \|\| nPDG==((-uPDG)/100)%10) ) Ord=-1; // LR/RL
3710	else G4cerr<<"-Warning-G4QIonIonCollision::AnnihilationOrder: Wrong pLaDiQ, L="<<uPDG
3711	<<", R="<<mPDG<<", cL="<<sPDG<<", cR="<<nPDG<<G4endl;
3712	}
3713	else if (mPDG <-7) // Fuse pRaDiQ
3714	{
3715	#ifdef debug
3716	G4cout<<"G4QIonIonCollision::AnnihOrder:pRaDiQ, sPDG="<<sPDG<<",nPDG="<<nPDG<<G4endl;
3717	#endif
3718	if ( sPDG>0 && (sPDG==(-mPDG)/1000 \|\| sPDG==((-mPDG)/100)%10) ) Ord=-1; // LR/RL
3719	else if( nPDG>0 && (nPDG==(-mPDG)/1000 \|\| nPDG==((-mPDG)/100)%10) ) Ord= 1; // LL/RR
3720	else G4cerr<<"-Warning-G4QIonIonCollision::AnnihilationOrder: Wrong pRaDiQ, L="<<uPDG
3721	<<", R="<<mPDG<<", cL="<<sPDG<<", cR="<<nPDG<<G4endl;
3722	}
3723	else if( (sPDG<0 && (-sPDG==mPDG \|\| -sPDG==uPDG) ) \|\|
3724	(nPDG<0 && (-nPDG==mPDG \|\| -nPDG==uPDG) ) ) Ord= 2; // @@ Annihilation fusion
3725	#ifdef debug
3726	else G4cout<<"-Warning-G4QIonIonCollision::AnnihilatOrder:Wrong23StringFusion"<<G4endl;
3727	G4cout<<"G4QIonIonCollision::AnnihilationOrder: Ord="<<Ord<<",sPDG="<<sPDG<<",nPDG="
3728	<<nPDG<<", uPDG="<<uPDG<<",mPDG="<<mPDG<<G4endl;
3729	#endif
3730	}
3731	else if ( LS==3 && MPS==2 )
3732	{
3733	if (sPDG > 7) // Fuse cLDiQ
3734	{
3735	#ifdef debug
3736	G4cout<<"G4QIonIonCollision::AnnihOrder:cLDiQ, uPDG="<<uPDG<<", mPDG="<<mPDG<<G4endl;
3737	#endif
3738	if ( uPDG<0 && (-uPDG==sPDG/1000 \|\| -uPDG==(sPDG/100)%10) ) Ord= 1; // LL/RR
3739	else if( mPDG<0 && (-mPDG==sPDG/1000 \|\| -mPDG==(sPDG/100)%10) ) Ord=-1; // LR/RL
3740	else G4cerr<<"-Warning-G4QIonIonCollision::AnnihilationOrder: Wrong cLDiQ, L="<<uPDG
3741	<<",R="<<mPDG<<",cL="<<sPDG<<",cR="<<nPDG<<G4endl;
3742	}
3743	else if (nPDG > 7) // Fuse cRDiQ
3744	{
3745	#ifdef debug
3746	G4cout<<"G4QIonIonCollision::AnnihOrder:cRDiQ, uPDG="<<uPDG<<", mPDG="<<mPDG<<G4endl;
3747	#endif
3748	if ( uPDG<0 && (-uPDG==nPDG/1000 \|\| -uPDG==(nPDG/100)%10) ) Ord=-1; // LR/RL
3749	else if( mPDG<0 && (-mPDG==nPDG/1000 \|\| -mPDG==(nPDG/100)%10) ) Ord= 1; // LL/RR
3750	else G4cerr<<"-Warning-G4QIonIonCollision::AnnihilationOrder: Wrong cRDiQ, L="<<uPDG
3751	<<",R="<<mPDG<<",cL="<<sPDG<<",cR="<<nPDG<<G4endl;
3752	}
3753	else if (sPDG <-7) // Fuse cLaDiQ
3754	{
3755	#ifdef debug
3756	G4cout<<"G4QIonIonCollision::AnnihOrder:cLaDiQ, uPDG="<<uPDG<<",mPDG="<<mPDG<<G4endl;
3757	#endif
3758	if ( uPDG>0 && (uPDG==(-sPDG)/1000 \|\| uPDG==((-sPDG)/100)%10) ) Ord= 1; // LL/RR
3759	else if( mPDG>0 && (mPDG==(-sPDG)/1000 \|\| mPDG==((-sPDG)/100)%10) ) Ord=-1; // LR/RL
3760	else G4cerr<<"-Warning-G4QIonIonCollision::AnnihilationOrder: Wrong cLaDiQ, L="<<uPDG
3761	<<", R="<<mPDG<<", cL="<<sPDG<<", cR="<<nPDG<<G4endl;
3762	}
3763	else if (nPDG <-7) // Fuse cRaDiQ
3764	{
3765	#ifdef debug
3766	G4cout<<"G4QIonIonCollision::AnnihOrder:cRaDiQ, uPDG="<<uPDG<<",mPDG="<<mPDG<<G4endl;
3767	#endif
3768	if ( uPDG>0 && (uPDG==(-nPDG)/1000 \|\| uPDG==((-nPDG)/100)%10) ) Ord=-1; // LR/RL
3769	else if( mPDG>0 && (mPDG==(-nPDG)/1000 \|\| mPDG==((-nPDG)/100)%10) ) Ord= 1; // LL/RR
3770	else G4cerr<<"-Warning-G4QIonIonCollision::AnnihilationOrder: Wrong cRaDiQ, L="<<uPDG
3771	<<", R="<<mPDG<<", cL="<<sPDG<<", cR="<<nPDG<<G4endl;
3772	}
3773	else if( (uPDG<0 && (-uPDG==sPDG \|\| -uPDG==nPDG) ) \|\|
3774	(mPDG<0 && (-mPDG==sPDG \|\| -mPDG==nPDG) ) ) Ord=2; // @@ Annihilation fusion
3775	#ifdef debug
3776	else G4cout<<"-Warning-G4QIonIonCollision::AnnihilatOrder:Wrong32StringFusion"<<G4endl;
3777	G4cout<<"G4QIonIonCollision::AnnihilationOrder: Ord="<<Ord<<",sPDG="<<sPDG<<",nPDG="
3778	<<nPDG<<", uPDG="<<uPDG<<",mPDG="<<mPDG<<G4endl;
3779	#endif
3780	}
3781	else if ( (LS==2 && MPS==4) \|\| (LS==4 && MPS==2) )
3782	{
3783	if (uPDG > 7) // Annihilate pLDiQ
3784	{
3785	#ifdef debug
3786	G4cout<<"G4QIonIonCollision::AnnihilOrder:pLDiQ,sPDG="<<sPDG<<",nPDG="<<nPDG<<G4endl;
3787	#endif
3788	if ( sPDG<0 && (-sPDG==uPDG/1000 \|\| -sPDG==(uPDG/100)%10) &&
3789	(nPDG==(-mPDG)/1000 \|\| nPDG==((-mPDG)/100)%10) ) Ord= 1; // LL/RR
3790	else if( nPDG<0 && (-nPDG==uPDG/1000 \|\| -nPDG==(uPDG/100)%10) &&
3791	(sPDG==(-mPDG)/1000 \|\| sPDG==((-mPDG)/100)%10) ) Ord=-1; // LR/RL
3792	else G4cerr<<"-Warning-G4QIonIonCollision::AnnihilationOrder: Wrong 24 pLDiQ, L="
3793	<<uPDG<<",R="<<mPDG<<",cL="<<sPDG<<",cR="<<nPDG<<G4endl;
3794	}
3795	else if (mPDG >7) // Annihilate pRDiQ
3796	{
3797	#ifdef debug
3798	G4cout<<"G4QIonIonCollision::AnnihilOrder:PRDiQ,sPDG="<<sPDG<<",nPDG="<<nPDG<<G4endl;
3799	#endif
3800	if ( sPDG<0 && (-sPDG==mPDG/1000 \|\| -sPDG==(mPDG/100)%10) &&
3801	(nPDG==(-uPDG)/1000 \|\| nPDG==((-uPDG)/100)%10) ) Ord=-1; // LR/RL
3802	else if( nPDG<0 && (-nPDG==mPDG/1000 \|\| -nPDG==(mPDG/100)%10) &&
3803	(sPDG==(-uPDG)/1000 \|\| sPDG==((-uPDG)/100)%10) ) Ord= 1; // LL/RR
3804	else G4cerr<<"-Warning-G4QIonIonCollision::AnnihilationOrder: Wrong 24 pLDiQ, L="
3805	<<uPDG<<",R="<<mPDG<<",cL="<<sPDG<<",cR="<<nPDG<<G4endl;
3806	}
3807	else if (sPDG > 7) // Annihilate cLDiQ
3808	{
3809	#ifdef debug
3810	G4cout<<"G4QIonIonCollision::AnnihilOrder:cLDiQ,uPDG="<<uPDG<<",mPDG="<<mPDG<<G4endl;
3811	#endif
3812	if ( uPDG<0 && (-uPDG==sPDG/1000 \|\| -uPDG==(sPDG/100)%10) &&
3813	(mPDG==(-nPDG)/1000 \|\| mPDG==((-nPDG)/100)%10) ) Ord= 1; // LL/RR
3814	else if( mPDG<0 && (-mPDG==sPDG/1000 \|\| -mPDG==(sPDG/100)%10) &&
3815	(uPDG==(-nPDG)/1000 \|\| uPDG==((-nPDG)/100)%10) ) Ord=-1; // LR/RL
3816	else G4cerr<<"-Warning-G4QIonIonCollision::AnnihilationOrder: Wrong 24 cLDiQ, L="
3817	<<uPDG<<",R="<<mPDG<<",cL="<<sPDG<<",cR="<<nPDG<<G4endl;
3818	}
3819	else if (nPDG > 7) // Annihilate cRDiQ
3820	{
3821	#ifdef debug
3822	G4cout<<"G4QIonIonCollision::AnnihilOrder:cRDiQ,uPDG="<<uPDG<<",mPDG="<<mPDG<<G4endl;
3823	#endif
3824	if ( uPDG<0 && (-uPDG==nPDG/1000 \|\| -uPDG==(nPDG/100)%10) &&
3825	(mPDG==(-sPDG)/1000 \|\| mPDG==((-sPDG)/100)%10) ) Ord=-1; // LR/RL
3826	else if( mPDG<0 && (-mPDG==nPDG/1000 \|\| -mPDG==(nPDG/100)%10) &&
3827	(uPDG==(-sPDG)/1000 \|\| uPDG==((-sPDG)/100)%10) ) Ord= 1; // LL/RR
3828	else G4cerr<<"-Warning-G4QIonIonCollision::AnnihilationOrder: Wrong 24 cRDiQ, L="
3829	<<uPDG<<",R="<<mPDG<<",cL="<<sPDG<<",cR="<<nPDG<<G4endl;
3830	}
3831	#ifdef debug
3832	else G4cout<<"-Warning-G4QIonIonCollision::AnnihilOrder: Wrong 24StringFusion"<<G4endl;
3833	G4cout<<"G4QIonIonCollision::AnnihilationOrder: Ord="<<Ord<<",sPDG="<<sPDG<<",nPDG="
3834	<<nPDG<<", uPDG="<<uPDG<<",mPDG="<<mPDG<<G4endl;
3835	#endif
3836	}
3837	else if ( LS==3 && MPS==3 )
3838	{
3839	if (uPDG > 7) // Annihilate pLDiQ
3840	{
3841	#ifdef debug
3842	G4cout<<"G4QIonIonCollision::AnnihilOrder:pLDiQ,sPDG="<<sPDG<<",nPDG="<<nPDG<<G4endl;
3843	#endif
3844	if ( sPDG<-7 && (-nPDG==uPDG/1000 \|\| -nPDG==(uPDG/100)%10) &&
3845	(mPDG==(-sPDG)/1000 \|\| mPDG==((-sPDG)/100)%10) ) Ord=-1; // LR/RL
3846	else if( nPDG<-7 && (-sPDG==uPDG/1000 \|\| -sPDG==(uPDG/100)%10) &&
3847	(mPDG==(-nPDG)/1000 \|\| mPDG==((-nPDG)/100)%10) ) Ord= 1; // LL/RR
3848	else G4cerr<<"-Warning-G4QIonIonCollision::AnnihilationOrder: Wrong 33 pLDiQ, L="
3849	<<uPDG<<",R="<<mPDG<<",cL="<<sPDG<<",cR="<<nPDG<<G4endl;
3850	}
3851	else if(mPDG > 7) // Annihilate pRDiQ
3852	{
3853	#ifdef debug
3854	G4cout<<"G4QIonIonCollision::AnnihilOrder:pRDiQ,sPDG="<<sPDG<<",nPDG="<<nPDG<<G4endl;
3855	#endif
3856	if ( sPDG<-7 && (-nPDG==mPDG/1000 \|\| -nPDG==(mPDG/100)%10) &&
3857	(uPDG==(-sPDG)/1000 \|\| uPDG==((-sPDG)/100)%10) ) Ord= 1; // LL/RR
3858	else if( nPDG<-7 && (-sPDG==mPDG/1000 \|\| -sPDG==(mPDG/100)%10) &&
3859	(uPDG==(-nPDG)/1000 \|\| uPDG==((-nPDG)/100)%10) ) Ord=-1; // LR/RL
3860	else G4cerr<<"-Warning-G4QIonIonCollision::AnnihilationOrder: Wrong33pRDiQ, L="<<uPDG
3861	<<",R="<<mPDG<<",cL="<<sPDG<<",cR="<<nPDG<<G4endl;
3862	}
3863	else if(sPDG > 7) // Annihilate cLDiQ
3864	{
3865	#ifdef debug
3866	G4cout<<"G4QIonIonCollision::AnnihilOrder:cLDiQ,uPDG="<<uPDG<<",mPDG="<<mPDG<<G4endl;
3867	#endif
3868	if ( uPDG<-7 && (-mPDG==sPDG/1000 \|\| -mPDG==(sPDG/100)%10) &&
3869	(nPDG==(-uPDG)/1000 \|\| nPDG==((-uPDG)/100)%10) ) Ord=-1; // LR/RL
3870	else if( mPDG<-7 && (-uPDG==sPDG/1000 \|\| -uPDG==(sPDG/100)%10) &&
3871	(nPDG==(-mPDG)/1000 \|\| nPDG==((-mPDG)/100)%10) ) Ord= 1; // LL/RR
3872	else G4cerr<<"-Warning-G4QIonIonCollision::AnnihilationOrder: Wrong 33 cLDiQ, L="
3873	<<uPDG<<",R="<<mPDG<<",cL="<<sPDG<<",cR="<<nPDG<<G4endl;
3874	}
3875	else if(nPDG > 7) // Annihilate cRDiQ
3876	{
3877	#ifdef debug
3878	G4cout<<"G4QIonIonCollision::AnnihilOrder:cRDiQ,uPDG="<<uPDG<<",mPDG="<<mPDG<<G4endl;
3879	#endif
3880	if ( uPDG<-7 && (-mPDG==nPDG/1000 \|\| -mPDG==(nPDG/100)%10) &&
3881	(nPDG==(-uPDG)/1000 \|\| nPDG==((-uPDG)/100)%10) ) Ord= 1; // LL/RR
3882	else if( mPDG<-7 && (-uPDG==nPDG/1000 \|\| -sPDG==(nPDG/100)%10) &&
3883	(sPDG==(-mPDG)/1000 \|\| sPDG==((-mPDG)/100)%10) ) Ord=-1; // LR/RL
3884	else G4cerr<<"-Warning-G4QIonIonCollision::AnnihilationOrder: Wrong 33 cRDiQ, L="
3885	<<uPDG<<",R="<<mPDG<<",cL="<<sPDG<<",cR="<<nPDG<<G4endl;
3886	}
3887	#ifdef debug
3888	else G4cout<<"-Warning-G4QIonIonCollision::AnnihilOrder: Wrong 33StringFusion"<<G4endl;
3889	G4cout<<"G4QIonIonCollision::AnnihilationOrder: Ord="<<Ord<<",sPDG="<<sPDG<<",nPDG="
3890	<<nPDG<<", uPDG="<<uPDG<<",mPDG="<<mPDG<<G4endl;
3891	#endif
3892	}
3893	return Ord;
3894	}
3895
3896	void G4QIonIonCollision::SwapPartons() // Swap string partons, if a string has negative M2
3897	{
3898	static const G4double baryM=800.; // Mass excess for baryons
3899	G4QStringVector::iterator ist;
3900	for(ist = strings.begin(); ist < strings.end(); ist++)
3901	{
3902	G4LorentzVector cS4M=(*ist)->Get4Momentum();
3903	G4double cSM2=cS4M.m2(); // Squared mass of the String
3904	if(cSM2<0.1) // Parton Swapping is needed
3905	{
3906	G4QParton* cLeft=(*ist)->GetLeftParton(); // Current String Left Parton
3907	G4QParton* cRight=(*ist)->GetRightParton(); // Current String Right Parton
3908	G4int cLPDG=cLeft->GetPDGCode();
3909	G4int cRPDG=cRight->GetPDGCode();
3910	G4LorentzVector cL4M=cLeft->Get4Momentum();
3911	G4LorentzVector cR4M=cRight->Get4Momentum();
3912	G4int cLT=cLeft->GetType();
3913	G4int cRT=cRight->GetType();
3914	G4QStringVector::iterator sst; // Selected partner string
3915	G4QStringVector::iterator pst; // LOOP iterator
3916	G4double maxM=-DBL_MAX; // Swapping providing the max mass
3917	G4int sDir=0; // Selected direction of swapping
3918	for(pst = strings.begin(); pst < strings.end(); pst++) if(pst != ist)
3919	{
3920	G4QParton* pLeft=(*pst)->GetLeftParton(); // Partner String Left Parton
3921	G4QParton* pRight=(*pst)->GetRightParton(); // Partner String Right Parton
3922	G4int pLPDG=pLeft->GetPDGCode();
3923	G4int pRPDG=pRight->GetPDGCode();
3924	G4LorentzVector pL4M=pLeft->Get4Momentum();
3925	G4LorentzVector pR4M=pRight->Get4Momentum();
3926	G4int pLT=pLeft->GetType();
3927	G4int pRT=pRight->GetType();
3928	G4double LM=0.;
3929	G4double RM=0.;
3930	if(((cLPDG<-7 \|\| (cLPDG>0 && cLPDG< 7) ) && (pLPDG<-7 \|\| (pLPDG>0 && pLPDG< 7) ))\|\|
3931	((cLPDG> 7 \|\| (cLPDG<0 && cLPDG>-7) ) && (pLPDG> 7 \|\| (pLPDG<0 && cLPDG>-7) )))
3932	{
3933	G4double pLM2=(cL4M+pR4M).m2(); // new partner M2
3934	G4double cLM2=(cR4M+pL4M).m2(); // new partner M2
3935	if(pLM2>0. && cLM2>0.)
3936	{
3937	G4double pLM=std::sqrt(pLM2);
3938	if(cLT+pRT==3) pLM-=baryM;
3939	G4double cLM=std::sqrt(cLM2);
3940	if(cRT+pLT==3) cLM-=baryM;
3941	LM=std::min(pLM2,cLM2);
3942	}
3943	}
3944	if(((cRPDG<-7 \|\| (cRPDG>0 && cRPDG< 7) ) && (pRPDG<-7 \|\| (pRPDG>0 && pRPDG< 7) ))\|\|
3945	((cRPDG> 7 \|\| (cRPDG<0 && cRPDG>-7) ) && (pRPDG> 7 \|\| (pRPDG<0 && cRPDG>-7) )) )
3946	{
3947	G4double pRM2=(cR4M+pL4M).m2(); // new partner M2
3948	G4double cRM2=(cL4M+pR4M).m2(); // new partner M2
3949	if(pRM2>0. && cRM2>0.)
3950	{
3951	G4double pRM=std::sqrt(pRM2);
3952	if(cRT+pLT==3) pRM-=baryM;
3953	G4double cRM=std::sqrt(cRM2);
3954	if(cLT+pRT==3) cRM-=baryM;
3955	RM=std::min(pRM,cRM);
3956	}
3957	}
3958	G4int dir=0;
3959	G4double sM=0.;
3960	if( LM && LM > RM )
3961	{
3962	dir= 1;
3963	sM=LM;
3964	}
3965	else if(RM)
3966	{
3967	dir=-1;
3968	sM=RM;
3969	}
3970	if(sM > maxM)
3971	{
3972	sst=pst;
3973	maxM=sM;
3974	sDir=dir;
3975	}
3976	}
3977	if(sDir)
3978	{
3979	G4QParton* pLeft=(*sst)->GetLeftParton(); // Partner String Left Parton
3980	G4QParton* pRight=(*sst)->GetRightParton(); // Partner String Right Parton
3981	G4QParton* swap=pLeft; // Prototype remember the partner Left
3982	if(sDir>0) // swap left partons
3983	{
3984	(*sst)->SetLeftParton(cLeft);
3985	(*ist)->SetLeftParton(swap);
3986	}
3987	else
3988	{
3989	swap=pRight;
3990	(*sst)->SetRightParton(cRight);
3991	(*ist)->SetRightParton(swap);
3992	}
3993	}
3994	#ifdef debug
3995	else G4cout<<"*G4QIonIonCollision::SwapPartons:Failed**,cLPDG="<<cLPDG<<",cRPDG="
3996	<<cRPDG<<",-->cM2="<<cSM2<<G4endl;
3997	#endif
3998
3999	}
4000	}
4001	}

Note: See TracBrowser for help on using the repository browser.

Download in other formats: