Geant4 9.6.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
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
56G4int G4QIonIonCollision::nCutMax=7;
57G4double G4QIonIonCollision::stringTension=1.*GeV/fermi;
58G4double G4QIonIonCollision::tubeDensity =1./fermi;
59// Parameters of diffractional fragmentation
60G4double G4QIonIonCollision::widthOfPtSquare=-0.75*GeV*GeV; // ptWidth2 forStringExcitation
61
62G4QIonIonCollision::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
1552G4QIonIonCollision::~G4QIonIonCollision()
1553{
1554 std::for_each(strings.begin(), strings.end(), DeleteQString() );
1555}
1556
1557G4QHadronVector* 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
1791void 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);//!! Try to fragment the new String !!
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);//! Try to fragment the new String !
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
3308G4bool 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
3415G4bool 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
3543void 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
3551G4double 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)
3572G4ThreeVector 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
3586G4int 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)
3651std::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)
3671G4int 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
3903void 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}
LT
@ LT
Definition: Evaluator.cc:66
FatalException
@ FatalException
Definition: G4ExceptionSeverity.hh:60
G4LorentzVector
CLHEP::HepLorentzVector G4LorentzVector
Definition: G4LorentzVector.hh:35
G4QHadronVector
std::vector< G4QHadron * > G4QHadronVector
Definition: G4QHadronVector.hh:56
G4QInteractionVector
std::vector< G4QInteraction * > G4QInteractionVector
Definition: G4QInteractionVector.hh:54
G4QPartonPairVector
std::vector< G4QPartonPair * > G4QPartonPairVector
Definition: G4QPartonPairVector.hh:46
G4ThreeVector
CLHEP::Hep3Vector G4ThreeVector
Definition: G4ThreeVector.hh:42
G4double
double G4double
Definition: G4Types.hh:64
G4int
int G4int
Definition: G4Types.hh:66
G4bool
bool G4bool
Definition: G4Types.hh:67
G4endl
#define G4endl
Definition: G4ios.hh:52
G4cerr
G4DLLIMPORT std::ostream G4cerr
G4cout
G4DLLIMPORT std::ostream G4cout
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:53
CLHEP::Hep3Vector::x
double x() const
CLHEP::Hep3Vector::mag2
double mag2() const
CLHEP::Hep3Vector::y
double y() const
CLHEP::Hep3Vector::dot
double dot(const Hep3Vector &) const
CLHEP::Hep3Vector::setZ
void setZ(double)
CLHEP::HepLorentzRotation::rotateY
HepLorentzRotation & rotateY(double delta)
Definition: LorentzRotation.cc:162
CLHEP::HepLorentzRotation::rotateZ
HepLorentzRotation & rotateZ(double delta)
Definition: LorentzRotation.cc:174
CLHEP::HepLorentzRotation::inverse
HepLorentzRotation inverse() const
CLHEP::HepLorentzVector::theta
double theta() const
CLHEP::HepLorentzVector::boostVector
Hep3Vector boostVector() const
Definition: LorentzVector.cc:177
CLHEP::HepLorentzVector::boost
HepLorentzVector & boost(double, double, double)
Definition: LorentzVector.cc:59
CLHEP::HepLorentzVector::vect
Hep3Vector vect() const
CLHEP::HepLorentzVector::minus
double minus() const
CLHEP::HepLorentzVector::transform
HepLorentzVector & transform(const HepRotation &)
G4QCHIPSWorld::GetQParticle
G4QParticle * GetQParticle(G4int PDG) const
Definition: G4QCHIPSWorld.hh:87
G4QCHIPSWorld::Get
static G4QCHIPSWorld * Get()
Definition: G4QCHIPSWorld.cc:79
G4QChipolino::GetQPDG1
G4QPDGCode GetQPDG1()
Definition: G4QChipolino.hh:91
G4QChipolino::GetQPDG2
G4QPDGCode GetQPDG2()
Definition: G4QChipolino.hh:92
G4QContent::GetSPDGCode
G4int GetSPDGCode() const
Definition: G4QContent.cc:1204
G4QContent::AddParton
G4int AddParton(G4int pPDG) const
Definition: G4QContent.cc:1655
G4QEnvironment::AddQuasmon
void AddQuasmon(G4Quasmon *Q)
Definition: G4QEnvironment.cc:10522
G4QEnvironment::Fragment
G4QHadronVector * Fragment()
Definition: G4QEnvironment.cc:4512
G4QHadron::Get4Momentum
G4LorentzVector Get4Momentum() const
Definition: G4QHadron.hh:79
G4QHadron::GetMass
G4double GetMass() const
Definition: G4QHadron.hh:176
G4QHadron::GetBaryonNumber
G4int GetBaryonNumber() const
Definition: G4QHadron.hh:181
G4QHadron::IncrementCollisionCount
void IncrementCollisionCount(G4int aCount)
Definition: G4QHadron.hh:104
G4QHadron::GetCharge
G4int GetCharge() const
Definition: G4QHadron.hh:179
G4QHadron::GetNextParton
G4QParton * GetNextParton()
Definition: G4QHadron.cc:1614
G4QHadron::GetPDGCode
G4int GetPDGCode() const
Definition: G4QHadron.hh:170
G4QHadron::DecayIn2
G4bool DecayIn2(G4LorentzVector &f4Mom, G4LorentzVector &s4Mom)
Definition: G4QHadron.cc:544
G4QHadron::Boost
void Boost(const G4LorentzVector &theBoost)
Definition: G4QHadron.cc:1293
G4QHadron::GetPosition
const G4ThreeVector & GetPosition() const
Definition: G4QHadron.hh:182
G4QHadron::GetColor
std::list< G4QParton * > GetColor()
Definition: G4QHadron.hh:93
G4QHadron::GetAntiColor
std::list< G4QParton * > GetAntiColor()
Definition: G4QHadron.hh:94
G4QHadron::GetNextAntiParton
G4QParton * GetNextAntiParton()
Definition: G4QHadron.cc:1622
G4QHadron::GetMass2
G4double GetMass2() const
Definition: G4QHadron.hh:177
G4QHadron::GetQC
G4QContent GetQC() const
Definition: G4QHadron.hh:173
G4QHadron::Set4Momentum
void Set4Momentum(const G4LorentzVector &aMom)
Definition: G4QHadron.hh:187
G4QInteraction::GetNumberOfSoftCollisions
G4int GetNumberOfSoftCollisions()
Definition: G4QInteraction.hh:59
G4QInteraction::SetNumberOfSoftCollisions
void SetNumberOfSoftCollisions(G4int nofSoft)
Definition: G4QInteraction.hh:65
G4QInteraction::GetProjectile
G4QHadron * GetProjectile() const
Definition: G4QInteraction.hh:57
G4QInteraction::GetTarget
G4QHadron * GetTarget() const
Definition: G4QInteraction.hh:58
G4QInteraction::SetTarget
void SetTarget(G4QHadron *aTarget)
Definition: G4QInteraction.hh:63
G4QInteraction::SetNumberOfDINRCollisions
void SetNumberOfDINRCollisions(G4int nofDINR)
Definition: G4QInteraction.hh:66
G4QInteraction::SetNumberOfDiffractiveCollisions
void SetNumberOfDiffractiveCollisions(G4int nofDiff)
Definition: G4QInteraction.hh:68
G4QIonIonCollision::G4QIonIonCollision
G4QIonIonCollision(G4QNucleus &pNucleus, const G4QNucleus &tNucleus)
Definition: G4QIonIonCollision.cc:62
G4QIonIonCollision::AnnihilationOrder
G4int AnnihilationOrder(G4int LS, G4int MS, G4int uP, G4int mP, G4int sP, G4int nP)
Definition: G4QIonIonCollision.cc:3671
G4QIonIonCollision::GaussianPt
G4ThreeVector GaussianPt(G4double widthSquare, G4double maxPtSquare) const
Definition: G4QIonIonCollision.cc:3572
G4QIonIonCollision::ExciteSingDiffParticipants
G4bool ExciteSingDiffParticipants(G4QHadron *aPartner, G4QHadron *bPartner) const
Definition: G4QIonIonCollision.cc:3415
G4QIonIonCollision::ReducePair
std::pair< G4int, G4int > ReducePair(G4int P1, G4int P2) const
Definition: G4QIonIonCollision.cc:3651
G4QIonIonCollision::Fragment
G4QHadronVector * Fragment()
Definition: G4QIonIonCollision.cc:1557
G4QIonIonCollision::SetParameters
static void SetParameters(G4int nC, G4double strTens, G4double tubeDens, G4double SigPt)
Definition: G4QIonIonCollision.cc:3543
G4QIonIonCollision::SumPartonPDG
G4int SumPartonPDG(G4int PDG1, G4int PFG2) const
Definition: G4QIonIonCollision.cc:3586
G4QIonIonCollision::ChooseX
G4double ChooseX(G4double Xmin, G4double Xmax) const
Definition: G4QIonIonCollision.cc:3551
G4QIonIonCollision::ExciteDiffParticipants
G4bool ExciteDiffParticipants(G4QHadron *aPartner, G4QHadron *bPartner) const
Definition: G4QIonIonCollision.cc:3308
G4QNucleus::InitByPDG
void InitByPDG(G4int newPDG)
Definition: G4QNucleus.cc:371
G4QNucleus::GetN
G4int GetN() const
Definition: G4QNucleus.hh:71
G4QNucleus::GetZ
G4int GetZ() const
Definition: G4QNucleus.hh:70
G4QNucleus::ChooseImpactXandY
std::pair< G4double, G4double > ChooseImpactXandY(G4double maxImpact)
Definition: G4QNucleus.cc:3551
G4QNucleus::GetPDG
G4int GetPDG() const
Definition: G4QNucleus.hh:69
G4QNucleus::GetOuterRadius
G4double GetOuterRadius()
Definition: G4QNucleus.cc:3947
G4QNucleus::StartLoop
G4bool StartLoop()
Definition: G4QNucleus.cc:3982
G4QNucleus::SubtractNucleon
void SubtractNucleon(G4QHadron *pNucleon)
Definition: G4QNucleus.cc:612
G4QNucleus::DoLorentzBoost
void DoLorentzBoost(const G4LorentzVector &theBoost)
Definition: G4QNucleus.hh:155
G4QNucleus::GetA
G4int GetA() const
Definition: G4QNucleus.hh:73
G4QNucleus::DoLorentzRotation
void DoLorentzRotation(const G4LorentzRotation &theLoRot)
Definition: G4QNucleus.hh:160
G4QNucleus::GetQCZNS
G4QContent GetQCZNS() const
Definition: G4QNucleus.hh:83
G4QNucleus::DeleteNucleons
void DeleteNucleons()
Definition: G4QNucleus.cc:697
G4QNucleus::GetThickness
G4double GetThickness(G4double b)
Definition: G4QNucleus.cc:4048
G4QNucleus::GetNucleons4Momentum
G4LorentzVector GetNucleons4Momentum()
Definition: G4QNucleus.hh:107
G4QNucleus::Init3D
void Init3D()
Definition: G4QNucleus.cc:3910
G4QNucleus::GetNextNucleon
G4QHadron * GetNextNucleon()
Definition: G4QNucleus.hh:100
G4QNucleus::GetGSMass
G4double GetGSMass() const
Definition: G4QNucleus.hh:82
G4QPDGCode::GetPDGCode
G4int GetPDGCode() const
Definition: G4QPDGCode.hh:326
G4QPDGCode::GetMass
G4double GetMass()
Definition: G4QPDGCode.cc:693
G4QPDGCode::MakeTwoBaryons
std::pair< G4int, G4int > MakeTwoBaryons(G4int L1, G4int L2, G4int R1, G4int R2)
Definition: G4QPDGCode.cc:2413
G4QParticle::MinMassOfFragm
G4double MinMassOfFragm()
Definition: G4QParticle.hh:113
G4QPartonPair::GetCollisionType
G4int GetCollisionType()
Definition: G4QPartonPair.hh:72
G4QParton::SetPDGCode
void SetPDGCode(G4int aPDG)
Definition: G4QParton.cc:130
G4QParton::GetQC
G4QContent GetQC() const
Definition: G4QParton.hh:82
G4QParton::GetType
const G4int & GetType() const
Definition: G4QParton.hh:88
G4QParton::GetPDGCode
G4int GetPDGCode() const
Definition: G4QParton.hh:81
G4QParton::Get4Momentum
const G4LorentzVector & Get4Momentum() const
Definition: G4QParton.hh:84
G4QParton::Set4Momentum
void Set4Momentum(const G4LorentzVector &aMomentum)
Definition: G4QParton.hh:75
G4QParton::ReduceDiQADiQ
G4bool ReduceDiQADiQ(G4QParton *d1, G4QParton *d2)
Definition: G4QParton.cc:246
G4QProbability::GetPomDiffProbability
G4double GetPomDiffProbability(const G4double s_value, const G4double imp2)
Definition: G4QProbability.hh:78
G4QProbability::GetCutPomProbability
G4double GetCutPomProbability(const G4double s, const G4double ip2, const G4int nPom)
Definition: G4QProbability.cc:74
G4QProbability::GetPomInelProbability
G4double GetPomInelProbability(const G4double s_value, const G4double imp2)
Definition: G4QProbability.hh:88
G4QString::FragmentString
G4QHadronVector * FragmentString(G4bool QL)
Definition: G4QString.cc:348
G4QString::GetLeftParton
G4QParton * GetLeftParton() const
Definition: G4QString.hh:87
G4QString::GetCharge
G4int GetCharge() const
Definition: G4QString.hh:92
G4QString::GetDirection
G4int GetDirection() const
Definition: G4QString.hh:89
G4QString::Get4Momentum
G4LorentzVector Get4Momentum() const
Definition: G4QString.cc:124
G4QString::GetBaryonNumber
G4int GetBaryonNumber() const
Definition: G4QString.hh:93
G4QString::GetRightParton
G4QParton * GetRightParton() const
Definition: G4QString.hh:88
G4QString::GetQC
G4QContent GetQC() const
Definition: G4QString.hh:91
G4QString::Boost
void Boost(G4ThreeVector &Velocity)
Definition: G4QString.cc:152
G4Quasmon::GetQC
G4QContent GetQC() const
Definition: G4Quasmon.hh:162
G4Quasmon::Set4Momentum
void Set4Momentum(G4LorentzVector Q4M)
Definition: G4Quasmon.hh:96
G4Quasmon::DecayQHadron
G4QHadronVector * DecayQHadron(G4QHadron *hadron)
Definition: G4Quasmon.cc:5779
G4Quasmon::Boost
void Boost(const G4LorentzVector &theBoost)
Definition: G4Quasmon.cc:6202
G4Quasmon::Get4Momentum
G4LorentzVector Get4Momentum() const
Definition: G4Quasmon.hh:161
G4Quasmon::SetQC
void SetQC(G4QContent QQC)
Definition: G4Quasmon.hh:97
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
sqr
T sqr(const T &x)
Definition: templates.hh:145
DBL_MAX
#define DBL_MAX
Definition: templates.hh:83