Context Navigation

source: trunk/source/processes/hadronic/models/chiral_inv_phase_space/fragmentation/src/G4QFragmentation.cc @ 1326

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

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

Download in other formats: