Context Navigation

G4QFragmentation.cc @ 1358

Last change on this file since 1358 was 1340, checked in by garnier, 14 years ago
update ti head
File size: 203.5 KB

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

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

Context Navigation

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

Download in other formats: