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