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