Context Navigation

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

Last change on this file since 1199 was 1199, checked in by garnier, 15 years ago
nvx fichiers dans CVS
File size: 197.6 KB

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

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

Download in other formats: