G4QIonIonCollision Class Reference

#include <G4QIonIonCollision.hh>


Public Member Functions

 G4QIonIonCollision (G4QNucleus &pNucleus, const G4QNucleus &tNucleus)
 ~G4QIonIonCollision ()
G4QHadronVectorFragment ()

Static Public Member Functions

static void SetParameters (G4int nC, G4double strTens, G4double tubeDens, G4double SigPt)

Protected Member Functions

G4bool ExciteDiffParticipants (G4QHadron *aPartner, G4QHadron *bPartner) const
G4bool ExciteSingDiffParticipants (G4QHadron *aPartner, G4QHadron *bPartner) const
std::pair< G4int, G4intReducePair (G4int P1, G4int P2) const
void Breeder ()
G4bool IsSingleDiffractive ()
G4int SumPartonPDG (G4int PDG1, G4int PFG2) const
G4double ChooseX (G4double Xmin, G4double Xmax) const
G4ThreeVector GaussianPt (G4double widthSquare, G4double maxPtSquare) const
G4int AnnihilationOrder (G4int LS, G4int MS, G4int uP, G4int mP, G4int sP, G4int nP)
void SwapPartons ()


Detailed Description

Definition at line 59 of file G4QIonIonCollision.hh.


Constructor & Destructor Documentation

G4QIonIonCollision::G4QIonIonCollision ( G4QNucleus pNucleus,
const G4QNucleus tNucleus 
)

Definition at line 62 of file G4QIonIonCollision.cc.

References G4QString::Boost(), G4QNucleus::ChooseImpactXandY(), DBL_MAX, G4QNucleus::DeleteNucleons(), G4QPartonPair::DIFFRACTIVE, G4QNucleus::DoLorentzBoost(), G4QNucleus::DoLorentzRotation(), ExciteDiffParticipants(), ExciteSingDiffParticipants(), FatalException, G4cerr, G4cout, G4endl, G4Exception(), G4UniformRand, G4QCHIPSWorld::Get(), G4QParton::Get4Momentum(), G4QString::Get4Momentum(), G4QHadron::Get4Momentum(), G4QNucleus::GetA(), G4QHadron::GetAntiColor(), G4QHadron::GetBaryonNumber(), G4QPartonPair::GetCollisionType(), G4QHadron::GetColor(), G4QProbability::GetCutPomProbability(), G4QNucleus::GetGSMass(), G4QHadron::GetMass(), G4QNucleus::GetN(), G4QHadron::GetNextAntiParton(), G4QNucleus::GetNextNucleon(), G4QHadron::GetNextParton(), G4QNucleus::GetNucleons4Momentum(), G4QNucleus::GetOuterRadius(), G4QParton::GetPDGCode(), G4QHadron::GetPDGCode(), G4QProbability::GetPomDiffProbability(), G4QProbability::GetPomInelProbability(), G4QHadron::GetPosition(), G4QHadron::GetQC(), G4QContent::GetSPDGCode(), G4QNucleus::GetThickness(), G4QParton::GetType(), G4QNucleus::GetZ(), G4QHadron::IncrementCollisionCount(), G4QNucleus::Init3D(), G4QNucleus::InitByPDG(), IsSingleDiffractive(), LT, G4QPDGCode::MakeTwoBaryons(), G4QPartonPair::PROJECTILE, G4QParton::ReduceDiQADiQ(), G4QParton::Set4Momentum(), G4QHadron::Set4Momentum(), G4QInteraction::SetNumberOfDiffractiveCollisions(), G4QInteraction::SetNumberOfDINRCollisions(), G4QInteraction::SetNumberOfSoftCollisions(), G4QParton::SetPDGCode(), G4QInteraction::SetTarget(), G4InuclParticleNames::sm, G4QPartonPair::SOFT, sqr(), G4QNucleus::StartLoop(), G4QNucleus::SubtractNucleon(), SwapPartons(), and G4QPartonPair::TARGET.

00063 {
00064   static const G4double mProt= G4QPDGCode(2212).GetMass(); // Mass of proton
00065   static const G4double mPi0= G4QPDGCode(111).GetMass();   // Mass of Pi0
00066   theWorld= G4QCHIPSWorld::Get();         // Get a pointer to the CHIPS World
00067   theResult = new G4QHadronVector;        // Define theResultingHadronVector
00068   G4bool stringsInitted=false;            // Strings are initiated
00069   G4LorentzRotation toZ;                  // Lorentz Transformation to the projectileSystem
00070   G4LorentzVector proj4M=pNucleus.Get4Momentum(); // Projectile nucleus 4-momentum in LS
00071   //G4double projM=pNucleus.GetGSMass();    // Ground state mass of the projectile nucleus
00072   G4double targM=tNucleus.GetGSMass();    // Ground state mass of the target nucleus
00073 #ifdef edebug
00074   G4cout<<"G4QIn::Constr:*Called*,pA="<<pNucleus<<proj4M<<",tA="<<tNucleus<<targM<<G4endl;
00075 #endif
00076   G4int pZ=pNucleus.GetZ();
00077   G4int pN=pNucleus.GetN();
00078   G4int pA=pZ+pN;
00079   G4int pPDG=90000000+pZ*1000+pN;
00080   G4int tZ=tNucleus.GetZ();
00081   G4int tN=tNucleus.GetN();
00082   G4int tA=tZ+tN;
00083   G4int tPDG=90000000+tZ*1000+tN;
00084   toZ.rotateZ(-proj4M.phi());
00085   toZ.rotateY(-proj4M.theta());
00086   G4LorentzVector zProj4M=toZ*proj4M;     // Proj 4-momentum in LS rotated to Z axis
00087   pNucleus.Set4Momentum(zProj4M);         // Now the projectile nucleus moves along Z axis
00088 #ifdef edebug
00089   G4int totChg=pZ+tZ;                     // Charge of the projectile+target for the CHECK
00090   G4int totBaN=pA+tA;                     // Baryon Number of Proj+Targ for CHECK
00091   G4LorentzVector tgLS4M(0.,0.,0.,targM); // Target 4-momentum in LS
00092   G4LorentzVector totLS4M=proj4M+tgLS4M;  // Total 4-momentum in LS
00093   G4LorentzVector totZLS4M=zProj4M+tgLS4M;// Total 4-momentum in LS with momentum along Z
00094   G4cout<<"-EMC-G4QIonIonCollision::Constr: tLS4M="<<totLS4M<<",tZLS4M="<<totZLS4M<<G4endl;
00095   // === From nere all consideration is made in the rotated LS frame (proj is along Z) ===
00096 #endif
00097   G4LorentzRotation toLab(toZ.inverse()); // Lorentz Transfornation "ZLS"->LS (at the end)
00098   G4QInteractionVector theInteractions;   // A vector of interactions (taken from the Body)
00099   G4QPartonPairVector  thePartonPairs;    // The parton pairs (taken from the Body)
00100   G4int                ModelMode=SOFT;    // The current model type (taken from the Body)
00101   theTargNucleus=G4QNucleus(tZ,tN);       // Create theTargNucleus to Move From LS to CM
00102   theTargNucleus.InitByPDG(tPDG);         // Reinit the Nucleus for the new Attempt?
00103   theTargNucleus.Init3D();                // 3D-initialisation(nucleons)ofTheTGNucleusClone
00104 #ifdef debug
00105   G4cout<<"G4QIonIonCollision::Constr:TargNuclMom="<<theTargNucleus.Get4Momentum()<<G4endl;
00106 #endif
00107   theProjNucleus=G4QNucleus(pZ,pN);       // Create theProjNucleus to Move From LS to CM
00108   theProjNucleus.InitByPDG(pPDG);         // Reinit the Nucleus for the new Attempt?
00109   theProjNucleus.Init3D();                // 3D-initialisation(nucleons)ofThePrNucleusClone
00110 #ifdef debug
00111   G4cout<<"G4QIonIonCollision::Constr:ProjNuclMom="<<theProjNucleus.Get4Momentum()<<G4endl;
00112 #endif
00113 #ifdef edebug
00114   G4LorentzVector sumP1=theProjNucleus.GetNucleons4Momentum();// Sum ofPrNucleons4M inRotLS
00115   G4LorentzVector sumT1=theTargNucleus.GetNucleons4Momentum();// Sum ofTgNucleons4M inRotLS
00116   G4cout<<"-EMC-G4QIonIonCollision::Construct: pNuc4M="<<sumP1<<", tNuc4M="<<sumT1<<G4endl;
00117 #endif
00118   G4ThreeVector theCMVelocity(0.,0.,0.);             // Prototype of the "CM" velocity
00119   G4ThreeVector theALVelocity(0.,0.,0.);             // Prototype of "CMAntiLab" velocity
00120   // Very important! Here the projectile 4M is recalculated in the ZLS (previously in LS)
00121   proj4M = pNucleus.Get4Momentum();                  // 4-mom of theProjectileHadron inZLS
00122   G4double pz_per_projectile = proj4M.pz()/pA;       // Momentum per nucleon
00123   G4double e_per_projectile = proj4M.e()/pA+targM/tA;// Add MassOfTargetNucleon
00124   G4double vz = pz_per_projectile/e_per_projectile;  // Speed of the "oneNuclenCM"
00125 #ifdef debug
00126   G4cout<<"G4QIonIonCollision::Construct: (КщеЯ)Projectile4M="<<proj4M<<", Vz="<<vz
00127         <<", pA="<<pA<<", pE="<<e_per_projectile<<G4endl;
00128 #endif
00129   theCMVelocity.setZ(vz);                            // CM (w/r to one nucleon) velosity
00130   theProjNucleus.DoLorentzBoost(-theCMVelocity);     // BoostProjNucleus to "CM"
00131   theTargNucleus.DoLorentzBoost(-theCMVelocity);     // BoostProjNucleus to "CM"
00132   G4LorentzVector prN4M=theProjNucleus.Get4Momentum();
00133   G4double avz=prN4M.pz()/prN4M.e();                 // Speed of AntiLabSys in CMS
00134   theALVelocity.setZ(avz);                           // CM (w/r to one nucleon) velosity
00135 #ifdef edebug
00136   G4LorentzVector sumP2=theProjNucleus.GetNucleons4Momentum();// Sum of Nucleons4M in RotCM
00137   G4LorentzVector sumT2=theTargNucleus.GetNucleons4Momentum();// Sum of Nucleons4M in RotCM
00138   G4cout<<"-EMC-G4QIonIonCollision::Construct: Boosted, vCM="<<vz<<", vAL="<<avz<<", tN4M="
00139         <<sumT2<<", pN4M="<<sumP2<<G4endl;
00140 #endif
00141   G4int maxCuts = 7;                                 // @@ Can be reduced for low energies
00142   //
00143   // >>----------->> Find collisions meeting collision conditions and the NN interaction XS
00144   //
00145   G4double outerRadius = theProjNucleus.GetOuterRadius()+theTargNucleus.GetOuterRadius();
00146   G4QProbability theProbability(2212);               // *** proj/targ nucleons ***
00147   // Clean up all previous interactions and reset the counters
00148 #ifdef debug
00149   G4cout<<"G4QIonIonCollision::Construct: theIntSize="<<theInteractions.size()<<G4endl;
00150 #endif
00151   G4int attCnt=-1;
00152   G4int maxAtt=27;
00153   // From here a loop over nucleons of the projectile Nucleus (@@ Z-sorting isn't implem!!)
00154   //
00155   G4QHadron* pNucleon; // Proto of the Projectile Nucleon pointer
00156   G4QHadron* tNucleon; // Proto of the Target Nucleon pointer
00157   G4QNucleus curProjNucleus;
00158   G4QNucleus curTargNucleus;
00159   while(!theInteractions.size() && ++attCnt < maxAtt)// TillAtLeastOneInteractionIsCreated
00160   {
00161     // std::for_each(theInteractions.begin(),theInteractions.end(), DeleteQInteraction());
00162     // Here we need to clean up the ProjNucleon and the TargNucleon in the Interactions !
00163     G4int nint=theInteractions.size();               // For the 1st attempt should be zero
00164     for(G4int ni=0; ni < nint; ++ni)
00165     {
00166       G4QInteraction* curInt = theInteractions[ni];
00167       delete curInt->GetProjectile();
00168       delete curInt->GetTarget();
00169       delete curInt;
00170     }
00171     theInteractions.clear();
00172     // Now we need to make a copy of the projectile and the target nuclei with 3D info
00173     if(attCnt)                                       // This is not theFirstAttempt->Clean
00174     {
00175       curProjNucleus.DeleteNucleons();
00176       curTargNucleus.DeleteNucleons();
00177     }
00178     curProjNucleus = theProjNucleus;
00179     curTargNucleus = theTargNucleus;
00180     // choose random impact parameter
00181     std::pair<G4double, G4double> theImpactParameter;
00182     theImpactParameter = curProjNucleus.ChooseImpactXandY(outerRadius);
00183     G4double impactX = theImpactParameter.first;
00184     G4double impactY = theImpactParameter.second;
00185 #ifdef debug
00186     G4cout<<"G4QIonIonCollision::Con:At#="<<attCnt<<",X="<<impactX<<",Y="<<impactY<<G4endl;
00187 #endif
00188     curProjNucleus.StartLoop();                      // Prepare Loop ovder nucleons
00189     G4int pnCnt=0;
00190     G4int pnA=curProjNucleus.GetA();
00191 #ifdef debu
00192     G4cout<<"G4QIonIonCollision::Constr: Before the WHILE over pNucleons,pA="<<pnA<<G4endl;
00193 #endif
00194     while((pNucleon=curProjNucleus.GetNextNucleon()) && pnCnt < pnA) // @@ can be for LOOP?
00195     {
00196       ++pnCnt;
00197       G4QInteractionVector curInteractions;          // A temporary vector of interactions
00198       G4LorentzVector pNuc4M=pNucleon->Get4Momentum();
00199       G4ThreeVector pNucR=pNucleon->GetPosition();   // Position of the pNucleon WRTo pNucl
00200       G4double pImpX = impactX+pNucR.x();
00201       G4double pImpY = impactY+pNucR.y();
00202       G4double prEn=proj4M.e();                      // For mesons
00203       G4int proB=pNucleus.GetBaryonNumber();
00204       if     (proB>0) prEn-=pNucleus.GetMass();      // For baryons
00205       else if(proB<0) prEn+=mProt;                   // For anti-baryons
00206 #ifdef debug
00207       G4double impactUsed = 0.;
00208       G4cout<<"G4QIonIonCollision::Construct: estimated energy, prEn="<<prEn<<G4endl;
00209 #endif
00210       G4int totalCuts = 0;
00211       // @@ This is a fake (random) s calculation @@ can be done inside the TARG-while
00212       G4int tnA=curTargNucleus.GetA();
00213       G4double pImp=std::sqrt(pImpX*pImpX+pImpY*pImpY);   
00214       G4double eflen=curTargNucleus.GetThickness(pImp); // EffectiveLength
00215       maxEn=eflen*stringTension;                     // maxAbsorbedEnergy in IndUnits=MeV
00216       maxNuc=static_cast<int>(eflen*tubeDensity+0.5);// max #0f involved nuclear nucleons
00217 #ifdef debug
00218       G4cout<<"G4QIonIonCollision::Con:pE="<<prEn<<" < mE="<<maxEn<<",mN="<<maxNuc<<G4endl;
00219 #endif
00220       if(prEn < maxEn)                               // Create DIN interaction and go out
00221       {
00222         G4QHadron* aProjectile = new G4QHadron(*pNucleon);// Copy selected PrNuc for String
00223         curTargNucleus.StartLoop();                  // Initialize newSelection forNucleons
00224         tNucleon=curTargNucleus.GetNextNucleon();    // Select a nucleon
00225         G4QHadron* aTarget = new G4QHadron(*tNucleon);// Copy selected nucleon for String
00226         G4QInteraction* anInteraction = new G4QInteraction(aProjectile);
00227         anInteraction->SetTarget(aTarget); 
00228         anInteraction->SetNumberOfDINRCollisions(1); // Consider this interaction as DINR
00229         curInteractions.push_back(anInteraction);    //--> now the Interaction is not empty
00230         curProjNucleus.DoLorentzBoost(-theALVelocity);// Boost theResPrNucleus toRotAntiLab
00231         curProjNucleus.SubtractNucleon(pNucleon);    // Pointer to the used ProjNucleon
00232         curProjNucleus.DoLorentzBoost(theALVelocity);// Boost theResProjNucleus back to CM
00233         curTargNucleus.DoLorentzBoost(theCMVelocity);// Boost theResTgNucleus to RotatedLS
00234         curTargNucleus.SubtractNucleon(tNucleon);    // Pointer to the used TargNucleon
00235         curTargNucleus.DoLorentzBoost(-theCMVelocity);// Boost theResTargNucleus back to CM
00236 #ifdef debug
00237         G4cout<<"G4QIonIonCollision::Construct: DINR interaction is created"<<G4endl;
00238 #endif
00239         break;                                       // Break the WHILE of the pNucleons
00240       }
00241       // LOOP over nuclei of the target nucleus to select collisions
00242       curTargNucleus.StartLoop();                    // To get the same nucleon
00243       G4int tnCnt = 0;
00244 #ifdef debu
00245       G4cout<<"G4QIonIonCollision::Constr: Before WHILE over tNucleons, tA="<<tnA<<G4endl;
00246 #endif
00247       while((tNucleon=curTargNucleus.GetNextNucleon()) && tnCnt<tnA && totalCuts<maxCuts)
00248       {
00249         ++tnCnt;
00250         G4LorentzVector tNuc4M=tNucleon->Get4Momentum();
00251 #ifdef debug
00252         G4cout<<"G4QIonIonCollision::Constr: OuterR="<<outerRadius<<", mC="<<maxCuts
00253               <<", pA="<<curProjNucleus<<", tA="<<curTargNucleus<<G4endl;
00254 #endif
00255         // Check the reaction threshold 
00256         G4double s_value = (tNuc4M + pNuc4M).mag2();   // Squared CM Energy of compound
00257         G4double ThresholdMass = pNucleon->GetMass() + tNucleon->GetMass();
00258 #ifdef debug
00259         G4cout<<"G4QInel::Constr: s="<<s_value<<", ThreshM="<<sqr(ThresholdMass)<<G4endl;
00260 #endif
00261         ModelMode = SOFT;                              // NOT-Diffractive hadronization
00262         if (s_value < 0.)                              // At ThP=0 is impossible(virtNucl)
00263         {
00264           G4cerr<<"*G4QInelast::Constr:s="<<s_value<<",pN4M="<<pNuc4M<<",tN4M="<<tNuc4M<<G4endl;
00265           G4Exception("G4QIonIonCollision::Construct:","72",FatalException,"LowEn(NegS)");
00266         }
00267         if (s_value < sqr(ThresholdMass))              // --> Only diffractive interaction
00268         {
00269 #ifdef debug
00270           G4cout<<"G4QIonIonCollision::Constr: ***OnlyDiffraction***, ThM="<<ThresholdMass
00271                 <<">sqrt(s)="<<std::sqrt(s_value)<<" -> only Diffraction is possible"<<G4endl;
00272 #endif
00273           ModelMode = DIFFRACTIVE;
00274         }
00275         G4ThreeVector tNucR=tNucleon->GetPosition(); // Position of the tNucleon WRTo tNucl
00276         G4double dImpX = pImpX-tNucR.x();
00277         G4double dImpY = pImpY-tNucR.y();
00278         G4double Distance2=dImpX*dImpX+dImpY*dImpY;
00279 #ifdef sdebug
00280         G4cout<<"G4QIonIonCollision::Construct: s="<<s_value<<", D2="<<Distance2<<G4endl;
00281 #endif
00282         // Needs to be moved to Probability class @@
00283         if(s_value<=10000.)
00284         {
00285           G4cout<<"-Warning-G4QIonIonCollision::Construct: s < .01 GeV^2, p4M="
00286                 <<pNucleon->Get4Momentum()<<",t4M="<<tNucleon->Get4Momentum()<<G4endl;
00287           continue;                                  // skip the rest of the targetNucleons
00288         }
00289         G4double Probability = theProbability.GetPomInelProbability(s_value, Distance2);// P_INEL
00290         // test for inelastic collision
00291 #ifdef sdebug
00292         G4cout<<"G4QIonIonCollision::Construct: Probubility="<<Probability<<G4endl;
00293 #endif
00294         G4double rndNumber = G4UniformRand();        // For the printing purpose
00295         // ModelMode = DIFFRACTIVE;
00296 #ifdef sdebug
00297         G4cout<<"G4QIonIonCollision::Construct: NLOOP prob="<<Probability<<", rndm="
00298               <<rndNumber<<", d="<<std::sqrt(Distance2)<<G4endl;
00299 #endif
00300         if (Probability > rndNumber) // Inelastic (diffractive or soft) interaction (JOB)
00301         {
00302           G4QHadron* aProjectile = new G4QHadron(*pNucleon);// Copy selected pNuc to String
00303           G4QHadron* aTarget = new G4QHadron(*tNucleon);// Copy selected tNucleon to String
00304 #ifdef edebug
00305           G4cout<<"--->EMC-->G4QIonIonCollision::Construct: TargNucleon is filled, 4M/PDG="
00306                 <<aTarget->Get4Momentum()<<aTarget->GetPDGCode()<<G4endl;
00307 #endif
00308           // Now the energy of the nucleons must be updated in CMS
00309           curProjNucleus.DoLorentzBoost(-theALVelocity);// Boost theResNucleus toRotatedLS
00310           curProjNucleus.SubtractNucleon(pNucleon);     // Pointer to the used nucleon
00311           curProjNucleus.DoLorentzBoost(theALVelocity); // Boost theResNucleus back to CM
00312           curTargNucleus.DoLorentzBoost(theCMVelocity); // Boost theResNucleus toRotatedLS
00313           curTargNucleus.SubtractNucleon(tNucleon);     // Pointer to the used nucleon
00314           curTargNucleus.DoLorentzBoost(-theCMVelocity);// Boost theResNucleus back to CM
00315           if((theProbability.GetPomDiffProbability(s_value,Distance2)/Probability >
00316               G4UniformRand() && ModelMode==SOFT ) || ModelMode==DIFFRACTIVE)
00317           { 
00318             // ------------->> diffractive interaction @@ IsSingleDiffractive called once
00319             if(IsSingleDiffractive()) ExciteSingDiffParticipants(aProjectile, aTarget);
00320             else                          ExciteDiffParticipants(aProjectile, aTarget);
00321             G4QInteraction* anInteraction = new G4QInteraction(aProjectile);
00322             anInteraction->SetTarget(aTarget); 
00323             anInteraction->SetNumberOfDiffractiveCollisions(1); // Why not increment? M.K.?
00324             curInteractions.push_back(anInteraction);//--> now theInteractions not empty
00325             // @@ Why not breake the NLOOP, if only one diffractive can happend?
00326             totalCuts++;                             // UpdateOfNucleons in't necessary
00327 #ifdef debug
00328             G4cout<<"G4QIonIonCollision::Constr:NLOOP DiffInteract,tC="<<totalCuts<<G4endl;
00329 #endif
00330           }
00331           else
00332           {
00333             // -------------->> nondiffractive = soft interaction
00334             // sample nCut+1 (cut Pomerons) pairs of strings can be produced
00335             G4int nCut;                              // Result in a chosen number of cuts
00336             G4double* running = new G4double[nCutMax];// @@ This limits the max cuts
00337             for(nCut = 0; nCut < nCutMax; nCut++)    // Calculates multiCut probabilities
00338             {
00339               running[nCut]= theProbability.GetCutPomProbability(s_value, Distance2, nCut+1);
00340               if(nCut) running[nCut] += running[nCut-1];// Sum up with the previous one
00341             }
00342             G4double random = running[nCutMax-1]*G4UniformRand();
00343             for(nCut = 0; nCut < nCutMax; nCut++) if(running[nCut] > random) break; // tNuc
00344             delete [] running;
00345 #ifdef debug
00346             G4cout<<"G4QIonIonCollision::Construct: NLOOP-Soft Chosen nCut="<<nCut<<G4endl;
00347 #endif
00348             // @@ If nCut>0 interaction with a few nucleons is possible
00349             // @@ nCut is found with big efforts and now nCut=0 ?? M.K. ?? !!
00350             //nCut = 0; // @@ in original code ?? @@
00351             aTarget->IncrementCollisionCount(nCut+1); // @@ What about multyNucleon target?
00352             aProjectile->IncrementCollisionCount(nCut+1);
00353             G4QInteraction* anInteraction = new G4QInteraction(aProjectile);
00354             anInteraction->SetTarget(aTarget);
00355             anInteraction->SetNumberOfSoftCollisions(nCut+1);
00356             curInteractions.push_back(anInteraction); // Now curInteractions are not empty
00357             totalCuts += nCut+1;
00358 #ifdef debug
00359             G4cout<<"G4QIonIonCollision::Constr:NLOOP SoftInteract,tC="<<totalCuts<<G4endl;
00360             impactUsed=Distance2;
00361 #endif
00362           }
00363         }// Probability selection
00364       } // End of While over target nucleons
00365       G4int nCurInt=curInteractions.size();
00366       for(G4int ni=0; ni < nCurInt; ++ni) theInteractions.push_back(curInteractions[ni]);
00367       curInteractions.clear();                      // Probably, not necessary...
00368     } // End of WHILE over projectile nucleons
00369   } // End of WHILE over attempts to find at least one nucleus-nucleus interaction
00370   theProjNucleus.DeleteNucleons();
00371   theTargNucleus.DeleteNucleons();
00372   theProjNucleus = curProjNucleus;
00373   theTargNucleus = curTargNucleus;
00374   curProjNucleus.DeleteNucleons();
00375   curTargNucleus.DeleteNucleons();
00376   G4int nInt=theInteractions.size();
00377 #ifdef debug
00378   G4cout<<"G4QIonIonCollision::Constr: #ofInteractions = "<<nInt<<", #OfDINR = "
00379         <<theInteractions[0]->GetNumberOfDINRCollisions()<<G4endl;
00380 #endif
00381   if(!nInt || (nInt==1 && theInteractions[0]->GetNumberOfDINRCollisions()==1)) // QFreeInel
00382   {
00383     G4QHadron* aTarget;
00384     G4QHadron* aProjectile;
00385     if(nInt)                                         // Take Targ/Proj from the Interaction
00386     {
00387         aTarget=theInteractions[0]->GetTarget();
00388              aProjectile=theInteractions[0]->GetProjectile();
00389       delete theInteractions[0];
00390       theInteractions.clear();
00391     }
00392     else                                             // Create a new target nucleon (?)
00393     {
00394       theProjNucleus.StartLoop();                    // To get the same nucleon from projec
00395       pNucleon=theProjNucleus.GetNextNucleon();      // Get theNucleon to create projectNuc
00396       aProjectile = new G4QHadron(*pNucleon);        // Copy selected pNucleon for String
00397       theProjNucleus.DoLorentzBoost(-theALVelocity); // Boost theResProjNucleus toRotatedLS
00398       theProjNucleus.SubtractNucleon(pNucleon);      // Pointer to SelProjNucleon to delete
00399       theProjNucleus.DoLorentzBoost(theALVelocity);  // Boost theResProNucleus back to CMS
00400       theTargNucleus.StartLoop();                    // To get the same nucleon from target
00401       tNucleon=theTargNucleus.GetNextNucleon();      // Get theNucleon to create targetNucl
00402       aTarget = new G4QHadron(*tNucleon);            // Copy selected tNucleon for String
00403       theTargNucleus.DoLorentzBoost(theCMVelocity);  // Boost theResTargNucleus toRotatedLS
00404       theTargNucleus.SubtractNucleon(tNucleon);      // Pointer to SelTargNucleon to delete
00405       theTargNucleus.DoLorentzBoost(-theCMVelocity); // Boost theResTargNucleus back to CMS
00406     }
00407     G4QContent QQC=aTarget->GetQC()+aProjectile->GetQC(); // QContent of the compound
00408     G4LorentzVector Q4M=aTarget->Get4Momentum()+aProjectile->Get4Momentum(); // 4-mom of Q
00409     delete aTarget;
00410     delete aProjectile;
00411     // @@ 4-Mom should be converted to LS for the TargQuasmon and to AL for the ProjQuasmon
00412     Q4M.boost(theCMVelocity);
00413     Q4M=toLab*Q4M;
00414     G4Quasmon* stringQuasmon = new G4Quasmon(QQC, Q4M);
00415     theQuasmons.push_back(stringQuasmon);
00416     theTargNucleus.DoLorentzBoost(theCMVelocity);   // BoostTheResidualNucleus toRotatedLS
00417     theTargNucleus.DoLorentzRotation(toLab);        // Recove Z-direction in LS ("LS"->LS)
00418     return;
00419   }
00420   //
00421   // ------------------ now build the parton pairs for the strings ------------------
00422   //
00423   for(G4int i=0; i<nInt; i++)
00424   {
00425     theInteractions[i]->SplitHadrons();
00426 #ifdef edebug
00427     G4QHadron* projH=theInteractions[i]->GetProjectile(); // Projectile of theInteraction
00428     G4QHadron* targH=theInteractions[i]->GetTarget();     // Target of the Interaction
00429     G4LorentzVector pSP(0.,0.,0.,0.);               // Sum of parton's 4mom's for proj
00430     G4LorentzVector tSP(0.,0.,0.,0.);               // Sum of parton's 4mom's for proj
00431     std::list<G4QParton*> projCP=projH->GetColor(); // Pointers to proj Color-partons
00432     std::list<G4QParton*> projAC=projH->GetAntiColor();// PointersTo projAntiColorPartons
00433     std::list<G4QParton*> targCP=targH->GetColor(); // Pointers to targ Color-partons
00434     std::list<G4QParton*> targAC=targH->GetAntiColor();// PointersTo targAntiColorPartons
00435     std::list<G4QParton*>::iterator picp = projCP.begin();
00436     std::list<G4QParton*>::iterator pecp = projCP.end();
00437     std::list<G4QParton*>::iterator piac = projAC.begin();
00438     std::list<G4QParton*>::iterator peac = projAC.end();
00439     std::list<G4QParton*>::iterator ticp = targCP.begin();
00440     std::list<G4QParton*>::iterator tecp = targCP.end();
00441     std::list<G4QParton*>::iterator tiac = targAC.begin();
00442     std::list<G4QParton*>::iterator teac = targAC.end();
00443     for(; picp!=pecp&& piac!=peac&& ticp!=tecp&& tiac!=teac; ++picp,++piac,++ticp,++tiac)
00444     {
00445       pSP+=(*picp)->Get4Momentum();
00446       pSP+=(*piac)->Get4Momentum();
00447       tSP+=(*ticp)->Get4Momentum();
00448       tSP+=(*tiac)->Get4Momentum();
00449     }
00450     G4cout<<"-EMC-G4QIonIonCollision::Construct: Interaction#"<<i<<",dP4M="
00451           <<projH->Get4Momentum()-pSP<<",dT4M="<<targH->Get4Momentum()-tSP<<G4endl;
00452 #endif
00453   }  
00454   // 
00455   // >>........>> make soft collisions (ordering is vital)
00456   //
00457   G4QInteractionVector::iterator it;
00458 #ifdef debug
00459   G4cout<<"G4QIonIonCollision::Constr: Creation ofSoftCollisionPartonPair STARTS"<<G4endl;
00460 #endif
00461   G4bool rep=true;
00462   while(rep && theInteractions.size())
00463   {
00464    for(it = theInteractions.begin(); it != theInteractions.end(); ++it)   
00465    {
00466     G4QInteraction* anIniteraction = *it;
00467     G4QPartonPair*  aPair=0;
00468     G4int nSoftCollisions = anIniteraction->GetNumberOfSoftCollisions();
00469 #ifdef debug
00470     G4cout<<"G4QIonIonCollision::Construct: #0f SoftCollisions ="<<nSoftCollisions<<G4endl;
00471 #endif
00472     if (nSoftCollisions)
00473     { 
00474       G4QHadron* pProjectile = anIniteraction->GetProjectile();
00475       G4QHadron* pTarget     = anIniteraction->GetTarget();
00476       for (G4int j = 0; j < nSoftCollisions; j++)
00477       {
00478         aPair = new G4QPartonPair(pTarget->GetNextParton(),
00479                                   pProjectile->GetNextAntiParton(),
00480                                   G4QPartonPair::SOFT, G4QPartonPair::TARGET);
00481         thePartonPairs.push_back(aPair); // A target pair (Why TAGRET?)
00482         aPair = new G4QPartonPair(pProjectile->GetNextParton(),
00483                                   pTarget->GetNextAntiParton(),
00484                                   G4QPartonPair::SOFT, G4QPartonPair::PROJECTILE);
00485         thePartonPairs.push_back(aPair); // A projectile pair (Why Projectile?)
00486 #ifdef debug
00487         G4cout<<"--->G4QIonIonCollision::Constr: SOFT, 2 parton pairs are filled"<<G4endl;
00488 #endif
00489       }  
00490       delete *it;
00491       it=theInteractions.erase(it);      // SoftInteractions're converted&erased
00492       if( it != theInteractions.begin() )// To avoid going below begin() (for Windows)
00493       {
00494         it--;
00495         rep=false;
00496       }
00497       else
00498       {
00499         rep=true;
00500         break;
00501       }
00502     }
00503     else rep=false;
00504    }
00505   }
00506 #ifdef debug
00507   G4cout<<"G4QIonIonCollision::Constr: -> Parton pairs for SOFT strings are made"<<G4endl;
00508 #endif  
00509   //
00510   // >>.............>> make the rest as the diffractive interactions
00511   //
00512   for(unsigned i = 0; i < theInteractions.size(); i++) // Interactions are reduced bySoft
00513   {
00514     // The double or single diffraction is defined by the presonce of proj/targ partons
00515     G4QInteraction* anIniteraction = theInteractions[i];
00516     G4QPartonPair* aPartonPair;
00517 #ifdef debug
00518     G4cout<<"G4QIonIonCollision::Constr: CreationOfDiffractivePartonPairs, i="<<i<<G4endl;
00519 #endif
00520     // the projectile diffraction parton pair is created first
00521     G4QHadron* aProjectile = anIniteraction->GetProjectile();
00522     G4QParton* aParton = aProjectile->GetNextParton();
00523     if (aParton)
00524     {
00525       aPartonPair = new G4QPartonPair(aParton, aProjectile->GetNextAntiParton(), 
00526                                       G4QPartonPair::DIFFRACTIVE,
00527                                       G4QPartonPair::PROJECTILE);
00528       thePartonPairs.push_back(aPartonPair);
00529 #ifdef debug
00530       G4cout<<"G4QIonIonCollision::Constr: proj Diffractive PartonPair is filled"<<G4endl;
00531 #endif
00532     }
00533     // then the target diffraction parton pair is created
00534     G4QHadron* aTarget = anIniteraction->GetTarget();
00535     aParton = aTarget->GetNextParton();
00536     if (aParton)
00537     {
00538       aPartonPair = new G4QPartonPair(aParton, aTarget->GetNextAntiParton(), 
00539                                       G4QPartonPair::DIFFRACTIVE, G4QPartonPair::TARGET);
00540       thePartonPairs.push_back(aPartonPair);
00541 #ifdef debug
00542       G4cout<<"G4QIonIonCollision::Constr: target Diffractive PartonPair's filled"<<G4endl;
00543 #endif
00544     }
00545   }
00546 #ifdef debug
00547   G4cout<<"G4QIonIonCollision::Construct: DiffractivePartonPairs are created"<<G4endl;
00548 #endif  
00549   //
00550   // >>.......>> clean-up  Interactions, if necessary
00551   //
00552   std::for_each(theInteractions.begin(),theInteractions.end(), DeleteQInteraction());
00553   theInteractions.clear();
00554 #ifdef debug
00555   G4cout<<"G4QIonIonCollision::Construct: Temporary objects are cleaned up"<<G4endl;
00556 #endif  
00557   // This function prepares theBoost for transformation of secondaries to LS (-ProjRot!)
00558   theProjNucleus.DoLorentzBoost(theCMVelocity);// Boost theResidualProjNucleus to RotatedLS
00559   theTargNucleus.DoLorentzBoost(theCMVelocity);// Boost theResidualTargNucleus to RotatedLS
00560   // @@ Nucleus isn't completely in LS, it needs the toZ (-ProjRot) rotation to consE/M
00561 #ifdef debug
00562   G4cout<<">>........>>G4QIonIonCollision::Construct: >>..>> Strings are created "<<G4endl;
00563 #endif
00564   G4QPartonPair* aPair;
00565   G4QString* aString=0;
00566   while(thePartonPairs.size()) // @@ At present noDifference in stringBuild (? M.K.)
00567   {
00568     aPair = thePartonPairs.back();           // Get the parton pair
00569     thePartonPairs.pop_back();               // Clean up thePartonPairPointer in the Vector
00570 #ifdef debug
00571     G4cout<<"G4QIonIonCollision::Construct:StringType="<<aPair->GetCollisionType()<<G4endl;
00572 #endif
00573     aString= new G4QString(aPair);
00574 #ifdef debug
00575     G4cout<<"G4QIonIonCollision::Construct:NewString4M="<<aString->Get4Momentum()<<G4endl;
00576 #endif
00577     aString->Boost(theCMVelocity);       // ! Strings are moved to ZLS when pushed !
00578     strings.push_back(aString);
00579     stringsInitted=true;
00580     delete aPair;
00581   } // End of the String Creation LOOP
00582 #ifdef edebug
00583   G4LorentzVector sum=theProjNucleus.Get4Momentum()+theTargNucleus.Get4Momentum(); // in LS
00584   G4int rChg=totChg-theProjNucleus.GetZ()-theTargNucleus.GetZ();
00585   G4int rBaN=totBaN-theProjNucleus.GetA()-theTargNucleus.GetA();
00586   G4int nStrings=strings.size();
00587   G4cout<<"-EMC-G4QIonIonCollision::Constr:#ofString="<<nStrings<<",tNuc4M="<<sum<<G4endl;
00588   for(G4int i=0; i<nStrings; i++)
00589   {
00590     G4QString* prString=strings[i];
00591     G4LorentzVector strI4M=prString->Get4Momentum();
00592     sum+=strI4M;
00593     G4int      sChg=prString->GetCharge();
00594     G4int      sBaN=prString->GetBaryonNumber();
00595     G4int      LPDG=prString->GetLeftParton()->GetPDGCode();
00596     G4int      RPDG=prString->GetRightParton()->GetPDGCode();
00597     G4QContent LQC =prString->GetLeftParton()->GetQC();
00598     G4QContent RQC =prString->GetRightParton()->GetQC();
00599     rChg-=sChg;
00600     rBaN-=sBaN;
00601     G4cout<<"-EMC-G4QIonIonCollision::Construct: String#"<<i<<", 4M="<<strI4M<<", LPDG="
00602           <<LPDG<<LQC<<",RPDG="<<RPDG<<RQC<<", Ch="<<sChg<<", BN="<<sBaN<<G4endl;
00603   }
00604   G4cout<<"-EMC-G4QInel::Constr: r4M="<<sum-totZLS4M<<", rC="<<rChg<<", rB="<<rBaN<<G4endl;
00605 #endif
00606   if(!stringsInitted)
00607   {
00608     G4cerr<<"*****G4QIonIonCollision::Construct:**** No strings are created *****"<<G4endl;
00609     G4Exception("G4QIonIonCollision::Construct:","72",FatalException,"No Strings created");
00610   }
00611 #ifdef debug
00612   G4cout<<"G4QIonIonCollision::Constr:BeforeRotation, #0fStrings="<<strings.size()<<G4endl;
00613 #endif
00614   //
00615   // ---------------- At this point the strings must be already created in "LS" -----------
00616   //
00617   for(unsigned astring=0; astring < strings.size(); astring++)
00618             strings[astring]->LorentzRotate(toLab); // Recove Z-direction in LS ("LS"->LS)
00619   theProjNucleus.DoLorentzRotation(toLab); // Recove Z-dir in LS ("LS"->LS) for ProjNucleus
00620   theTargNucleus.DoLorentzRotation(toLab); // Recove Z-dir in LS ("LS"->LS) for TargNucleus
00621   // Now everything is in LS system
00622 #ifdef edebug
00623   G4LorentzVector sm=theProjNucleus.Get4Momentum()+theTargNucleus.Get4Momentum();// NucInLS
00624   G4int rCg=totChg-theProjNucleus.GetZ()-theTargNucleus.GetZ();
00625   G4int rBC=totBaN-theProjNucleus.GetA()-theTargNucleus.GetA();
00626   G4int nStrs=strings.size();
00627   G4cout<<"-EMCLS-G4QIonIonCollision::Constr:#ofS="<<nStrings<<",Nuc4M(E=M)="<<sum<<G4endl;
00628   for(G4int i=0; i<nStrs; i++)
00629   {
00630     G4LorentzVector strI4M=strings[i]->Get4Momentum();
00631     sm+=strI4M;
00632     G4int sChg=strings[i]->GetCharge();
00633     rCg-=sChg;
00634     G4int sBaN=strings[i]->GetBaryonNumber();
00635     rBC-=sBaN;
00636     G4cout<<"-EMCLS-G4QInel::Construct:String#"<<i<<",4M="<<strI4M<<strI4M.m()<<",Charge="
00637           <<sChg<<",BaryN="<<sBaN<<G4endl;
00638   }
00639   G4cout<<"-EMCLS-...G4QInel::Constr: r4M="<<sm-totLS4M<<", rC="<<rCg<<",rB="<<rBC<<G4endl;
00640 #endif
00641   //
00642   // --- Strings are created, but we should try to get rid of negative mass strings -----
00643   //
00644   SwapPartons();
00645   //
00646   // --- Strings are created, but we should get rid of too light strings (Mmin+MPi0) -----
00647   //
00648   G4int problem=0;                                   // 0="no problem", incremented by ASIS
00649   G4QStringVector::iterator ist;
00650   G4bool con=true;
00651   while(con && strings.size())
00652   {
00653    for(ist = strings.begin(); ist < strings.end(); ++ist)
00654    {
00655     G4bool bad=true;
00656     G4LorentzVector cS4M=(*ist)->Get4Momentum();
00657     G4double cSM2=cS4M.m2();                         // Squared mass of the String
00658     G4QParton* cLeft=(*ist)->GetLeftParton();
00659     G4QParton* cRight=(*ist)->GetRightParton();
00660     G4int cLT=cLeft->GetType();
00661     G4int cRT=cRight->GetType();
00662     G4int cLPDG=cLeft->GetPDGCode();
00663     G4int cRPDG=cRight->GetPDGCode();
00664     G4int aLPDG=0;
00665     G4int aRPDG=0;
00666     if     (cLPDG > 7) aLPDG=  cLPDG/100;
00667     else if(cLPDG <-7) aLPDG=(-cLPDG)/100;
00668     if     (cRPDG > 7) aRPDG=  cRPDG/100;
00669     else if(cRPDG <-7) aRPDG=(-cRPDG)/100;
00670     G4int L1=0;
00671     G4int L2=0;
00672     if(aLPDG)
00673     {
00674       L1=aLPDG/10;
00675       L2=aLPDG%10;
00676     }
00677     G4int R1=0;
00678     G4int R2=0;
00679     if(aRPDG)
00680     {
00681       R1=aRPDG/10;
00682       R2=aRPDG%10;
00683     }
00684     G4double cSM=cSM2;
00685     if(cSM2>0.) cSM=std::sqrt(cSM2);
00686 #ifdef debug
00687     G4cout<<"G4QIonIonCollision::Construct: Needs Fusion? cLPDG="<<cLPDG<<",cRPDG="<<cRPDG
00688           <<",cM(cM2 if neg)="<<cSM<<G4endl;
00689 #endif
00690     if(cSM>0.)                                       // Mass can be calculated
00691     {
00692       G4bool single=true;
00693       G4double miM=0.;                               // Proto of the Min HadronString Mass
00694       if(cLT==2 && cRT==2)
00695       {
00696         if(L1!=R1 && L1!=R2 && L2!=R1 && L2!=R2)     // Unreducable DiQ-aDiQ
00697         {
00698           single=false;
00699           G4QPDGCode tmp;
00700           std::pair<G4int,G4int> paB=tmp.MakeTwoBaryons(L1, L2, R1, R2);
00701           miM=G4QPDGCode(paB.first).GetMass()+G4QPDGCode(paB.second).GetMass();
00702         }
00703       }
00704         if(single) miM=G4QPDGCode((*ist)->GetQC().GetSPDGCode()).GetMass() + mPi0;//MinHSMass
00705         //if(single) miM=G4QPDGCode((*ist)->GetQC().GetSPDGCode()).GetMass();//MinHSMass
00706 #ifdef debug
00707       G4cout<<"G4QInel::Const:*IsItGood? realM="<<std::sqrt(cSM2)<<" > GSM="<<miM<<G4endl;
00708 #endif
00709       if(std::sqrt(cSM2) > miM) bad=false;           // String is OK
00710     }
00711     if(bad)                                          // String should be merged with others
00712     {
00713 #ifdef debug
00714       G4cout<<"G4QInel::Const:TryFuse,L1="<<L1<<",L2="<<L2<<",R1="<<R1<<",R2="<<R2<<G4endl;
00715 #endif
00716       G4int cST=cLT+cRT;
00717       G4double excess=-DBL_MAX;                      // The value to be maximized excess M
00718       G4double maxiM2=-DBL_MAX;                      // The value to be maximized M2
00719       G4QStringVector::iterator sst;                 // Selected partner string
00720       G4QStringVector::iterator pst;
00721       G4int sLPDG=0;                                 // selectedLeft (like inStringPartner)
00722       G4int sRPDG=0;                                 // selectedRight(like inStringPartner)
00723       G4int sOrd=0;                                  // selected Order of PartonFusion
00724       G4bool minC=true;                              // for the case when M2<0
00725       if(cSM2>0.) minC=false;                        // If M2>0 already don'tSearchFor M2>0
00726       for(pst = strings.begin(); pst < strings.end(); pst++) if(pst != ist)
00727       {
00728         G4LorentzVector pS4M=(*pst)->Get4Momentum()+cS4M; // Summed 4-momentum
00729         G4int nLPDG=0;                               // new Left (like in theStringPartner)
00730         G4int nRPDG=0;                               // new Right(like in theStringPartner)
00731         G4double pExcess=-DBL_MAX;                   // Prototype of the excess
00732         G4double pSM2=pS4M.m2();                     // Squared mass of the Fused Strings
00733 #ifdef debug
00734         G4cout<<"->G4QIonIonCollision::Construct: sum4M="<<pS4M<<", M2="<<pSM2<<G4endl;
00735 #endif
00736         //if(pSM2>0.)                                  // The partner can be a candidate
00737         //{
00738         G4QParton* pLeft=(*pst)->GetLeftParton();
00739         G4QParton* pRight=(*pst)->GetRightParton();
00740         G4int pLT=pLeft->GetType();
00741         G4int pRT=pRight->GetType();
00742         G4int pLPDG=pLeft->GetPDGCode();
00743         G4int pRPDG=pRight->GetPDGCode();
00744         G4int LPDG=0;
00745         G4int RPDG=0;
00746         if     (pLPDG > 7) LPDG=  pLPDG/100;
00747         else if(pLPDG <-7) LPDG=(-pLPDG)/100;
00748         if     (pRPDG > 7) RPDG=  pRPDG/100;
00749         else if(pRPDG <-7) RPDG=(-pRPDG)/100;
00750         G4int pL1=0;
00751         G4int pL2=0;
00752         if(LPDG)
00753         {
00754           pL1=LPDG/10;
00755           pL2=LPDG%10;
00756         }
00757         G4int pR1=0;
00758         G4int pR2=0;
00759         if(RPDG)
00760         {
00761           pR1=RPDG/10;
00762           pR2=RPDG%10;
00763         }
00764         G4int pST=pLT+pRT;
00765 #ifdef debug
00766         G4cout<<"G4QInel::Construct: Partner/w pLPDG="<<pLPDG<<", pRPDG="<<pRPDG<<", pM2="
00767               <<pSM2<<G4endl;
00768 #endif
00769         // Different fromCompactAlrorithm ofStringFusionAfterDecay (no DiQaDiQ reduction)
00770         G4int tf=0;                                // Type combination flag
00771         G4int af=0;                                // Annihilatio combination flag
00772         if     (cST==2 && pST==2)                  // QaQ + QaQ = DiQaDiQ (always)
00773         {
00774           tf=1;
00775           af=1;
00776         }
00777         else if(cST==2 && pST==3)                  // QaQ + QDiQ/aQaDiQ = QDiQ/aQaDiQ (s)
00778         {
00779           tf=2;
00780           if     (pLPDG > 7 &&
00781                   ( (cLPDG<0 && (-cLPDG==pL1 || -cLPDG==pL2 || -cLPDG==pRPDG) ) ||
00782                     (cRPDG<0 && (-cRPDG==pL1 || -cRPDG==pL2 || -cRPDG==pRPDG) )
00783                   )
00784                  ) af=1;
00785           else if(pRPDG > 7 &&
00786                   ( (cLPDG<0 && (-cLPDG==pR1 || -cLPDG==pR2 || -cLPDG==pLPDG) ) ||
00787                     (cRPDG<0 && (-cRPDG==pR1 || -cRPDG==pR2 || -cRPDG==pLPDG) )
00788                   )
00789                  ) af=2;
00790           else if(pLPDG <-7 &&
00791                   ( (cLPDG>0 && ( cLPDG==pL1 || cLPDG==pL2 || cLPDG==-pRPDG) ) ||
00792                     (cRPDG>0 && ( cRPDG==pL1 || cRPDG==pL2 || cRPDG==-pRPDG) )
00793                   )
00794                  ) af=3;
00795           else if(pRPDG <-7 &&
00796                   ( (cLPDG>0 && ( cLPDG==pR1 || cLPDG==pR2 || cLPDG==-pLPDG) ) ||
00797                     (cRPDG>0 && ( cRPDG==pR1 || cRPDG==pR2 || cRPDG==-pLPDG) )
00798                   )
00799                  ) af=4;
00800 #ifdef debug
00801           else G4cout<<"G4QIonIonCollision::Constr:2(QaQ+QDiQ/aQaDiQ) Can't fuse"<<G4endl;
00802 #endif
00803         }
00804         else if(cST==3 && pST==2)                  // QDiQ/aQaDiQ + QaQ = QDiQ/aQaDiQ (s)
00805         {
00806           tf=3;
00807           if     (cLPDG > 7 &&
00808                   ( (pLPDG<0 && (-pLPDG==L1 || -pLPDG==L2) ) ||
00809                     (pRPDG<0 && (-pRPDG==L1 || -pRPDG==L2) )
00810                   )
00811                  ) af=1;
00812           else if(cRPDG > 7 &&
00813                   ( (pLPDG<0 && (-pLPDG==R1 || -pLPDG==R2) ) ||
00814                     (pRPDG<0 && (-pRPDG==R1 || -pRPDG==R2) )
00815                   )
00816                  ) af=2;
00817           else if(cLPDG <-7 &&
00818                   ( (pLPDG>0 && ( pLPDG==L1 || pLPDG==L2) ) ||
00819                     (pRPDG>0 && ( pRPDG==L1 || pRPDG==L2) )
00820                   )
00821                  ) af=3;
00822           else if(cRPDG <-7 &&
00823                   ( (pLPDG>0 && ( pLPDG==R1 || pLPDG==R2) ) ||
00824                     (pRPDG>0 && ( pRPDG==R1 || pRPDG==R2) )
00825                   )
00826                  ) af=4;
00827 #ifdef debug
00828           else G4cout<<"G4QIonIonCollision::Constr:3(QDiQ/aQaDiQ+QaQ) Can't fuse"<<G4endl;
00829 #endif
00830         }
00831         else if(cST==2 && pST==4)                  // QaQ + aDiQDiQ = QaQ (double)
00832         {
00833           tf=4;
00834           if     (pLPDG > 7) // pRPDG <-7
00835           {
00836             if     ( (-cLPDG==pL1 || -cLPDG==pL2) && (cRPDG==pR1 || cRPDG==pR2) ) af=1;
00837             else if( (-cRPDG==pL1 || -cRPDG==pL2) && (cLPDG==pR1 || cLPDG==pR2) ) af=2;
00838           }
00839           else if(pRPDG > 7) // pLPDG <-7
00840           {
00841             if     ( (-cRPDG==pR1 || -cRPDG==pR2) && (cLPDG==pL1 || cLPDG==pL2) ) af=3;
00842             else if( (-cLPDG==pR1 || -cLPDG==pR2) && (cRPDG==pL1 || cRPDG==pL2) ) af=4;
00843           }
00844 #ifdef debug
00845           else G4cout<<"-G4QIonIonCollision::Constr: 4 (QaQ+aQDiQDiQ) Can't fuse"<<G4endl;
00846 #endif
00847         }
00848         else if(cST==4 && pST==2)                  // aDiQDiQ + QaQ = QaQ (double)
00849         {
00850           tf=5;
00851           if     (cLPDG > 7) // cRPDG<-7
00852           {
00853             if     ( (-pLPDG==L1 || -pLPDG==L2) && (pRPDG==R1 || pRPDG==R2) ) af=1;
00854             else if( (-pRPDG==L1 || -pRPDG==L2) && (pLPDG==R1 || pLPDG==R2) ) af=2;
00855           }
00856           else if(cRPDG > 7) // cLPDG<-7
00857           {
00858             if     ( (-pRPDG==R1 || -pRPDG==R2) && (pLPDG==L1 || pLPDG==L2) ) af=3;
00859             else if( (-pLPDG==R1 || -pLPDG==R2) && (pRPDG==L1 || pRPDG==L2) ) af=4;
00860           }
00861 #ifdef debug
00862           else G4cout<<"-G4QIonIonCollision::Constr: 5 (aQDiQDiQ+QaQ) Can't fuse"<<G4endl;
00863 #endif
00864         }
00865         else if(cST==3 && pST==3)                  // QDiQ + aQaDiQ = QaQ (double)
00866         {
00867           tf=6;
00868           if(pLPDG > 7)
00869           {
00870             if     (cLPDG<-7 && (-cRPDG==pL1 || -cRPDG==pL2) && (pRPDG==L1 || pRPDG==L2))
00871               af=1;
00872             else if(cRPDG<-7 && (-cLPDG==pL1 || -cLPDG==pL2) && (pRPDG==R1 || pRPDG==R2))
00873               af=2;
00874           }
00875           else if(pRPDG > 7)
00876           {
00877             if     (cLPDG<-7 && (-cRPDG==pR1 || -cRPDG==pR2) && (pLPDG==L1 || pLPDG==L2))
00878               af=3;
00879             else if(cRPDG<-7 && (-cLPDG==pR1 || -cLPDG==pR2) && (pLPDG==R1 || pLPDG==R2))
00880               af=4;
00881           }
00882           else if(cLPDG > 7)
00883           {
00884             if     (pLPDG<-7 && (-pRPDG==L1 || -pRPDG==L2) && (cRPDG==pL1 || cRPDG==pL2))
00885               af=5;
00886             else if(pRPDG<-7 && (-pLPDG==L1 || -pLPDG==L2) && (cRPDG==pR1 || cRPDG==pR2))
00887               af=6;
00888           }
00889           else if(cRPDG > 7)
00890           {
00891             if     (pLPDG<-7 && (-pRPDG==R1 || -pRPDG==R2) && (cLPDG==pL1 || cLPDG==pL2))
00892               af=7;
00893             else if(pRPDG<-7 && (-pLPDG==R1 || -pLPDG==R2) && (cLPDG==pR1 || cLPDG==pR2))
00894               af=8;
00895           }
00896 #ifdef debug
00897           else G4cout<<"-G4QIonIonCollision::Constr: 6 (QDiQ+aQaDiQ) Can't fuse"<<G4endl;
00898 #endif
00899         }
00900 #ifdef debug
00901         G4cout<<"G4QIonIonCollision::Constr: **Possibility**, tf="<<tf<<",af="<<af<<G4endl;
00902 #endif
00903         if(tf && af)
00904         {
00905           // Strings can be fused, update the max excess and remember usefull switches
00906           G4int order=0;                           // LL/RR (1) or LR/RL (2) PartonFusion
00907           switch (tf)
00908           {
00909             case 1: // ------------------------------------> QaQ + QaQ = DiQaDiQ (always)
00910               if     (cLPDG > 0 && pLPDG > 0)
00911               {
00912                 order= 1;
00913                 if     (cLPDG > pLPDG) nLPDG=cLPDG*1000+pLPDG*100+1;
00914                 else if(cLPDG < pLPDG) nLPDG=pLPDG*1000+cLPDG*100+1;
00915                 else                   nLPDG=pLPDG*1000+cLPDG*100+3;
00916                 if     (cRPDG < pRPDG) nRPDG=cRPDG*1000+pRPDG*100-1;
00917                 else if(cRPDG > pRPDG) nRPDG=pRPDG*1000+cRPDG*100-1;
00918                 else                   nRPDG=pRPDG*1000+cRPDG*100-3;
00919               }
00920               else if(cLPDG < 0 && pLPDG < 0)
00921               {
00922                 order= 1;
00923                 if     (cRPDG > pRPDG) nRPDG=cRPDG*1000+pRPDG*100+1;
00924                 else if(cRPDG < pRPDG) nRPDG=pRPDG*1000+cRPDG*100+1;
00925                 else                   nRPDG=pRPDG*1000+cRPDG*100+3;
00926                 if     (cLPDG < pLPDG) nLPDG=cLPDG*1000+pLPDG*100-1;
00927                 else if(cLPDG > pLPDG) nLPDG=pLPDG*1000+cLPDG*100-1;
00928                 else                   nLPDG=pLPDG*1000+cLPDG*100-3;
00929               }
00930               else if(cRPDG > 0 && pLPDG > 0)
00931               {
00932                 order=-1;
00933                 if     (cRPDG > pLPDG) nLPDG=cRPDG*1000+pLPDG*100+1;
00934                 else if(cRPDG < pLPDG) nLPDG=pLPDG*1000+cRPDG*100+1;
00935                 else                   nLPDG=pLPDG*1000+cRPDG*100+3;
00936                 if     (cLPDG < pRPDG) nRPDG=cLPDG*1000+pRPDG*100-1;
00937                 else if(cLPDG > pRPDG) nRPDG=pRPDG*1000+cLPDG*100-1;
00938                 else                   nRPDG=pRPDG*1000+cLPDG*100-3;
00939               }
00940               else if(cRPDG < 0 && pLPDG < 0)
00941               {
00942                 order=-1;
00943                 if     (cLPDG > pRPDG) nRPDG=cLPDG*1000+pRPDG*100+1;
00944                 else if(cLPDG < pRPDG) nRPDG=pRPDG*1000+cLPDG*100+1;
00945                 else                   nRPDG=pRPDG*1000+cLPDG*100+3;
00946                 if     (cRPDG < pLPDG) nLPDG=cRPDG*1000+pLPDG*100-1;
00947                 else if(cRPDG > pLPDG) nLPDG=pLPDG*1000+cRPDG*100-1;
00948                 else                   nLPDG=pLPDG*1000+cRPDG*100-3;
00949               }
00950               break;
00951             case 2: // ------------------------> QaQ + QDiQ/aQaDiQ = QDiQ/aQaDiQ (single)
00952              switch (af)
00953              {
00954                case 1: // ....................... pLPDG > 7
00955                  if(cLPDG < 0)
00956                  {
00957                    order= 1;
00958                    if(-cLPDG==pRPDG)
00959                    {
00960                      nLPDG=pLPDG;
00961                      nRPDG=cRPDG;
00962                    }
00963                    else
00964                    {
00965                      if     (cRPDG > pRPDG) nRPDG=cRPDG*1000+pRPDG*100+1;
00966                      else if(cRPDG < pRPDG) nRPDG=pRPDG*1000+cRPDG*100+1;
00967                      else                   nRPDG=pRPDG*1000+cRPDG*100+3;
00968                      if  (-cLPDG == pL1)    nLPDG=pL2;
00969                      else                   nLPDG=pL1; // -cLPDG == pL2
00970                    }
00971                  }
00972                  else // cRPDG < 0
00973                  {
00974                    order=-1;
00975                    if(-cRPDG==pRPDG)
00976                    {
00977                      nLPDG=pLPDG;
00978                      nRPDG=cLPDG;
00979                    }
00980                    else
00981                    {
00982                      if     (cLPDG > pRPDG) nRPDG=cLPDG*1000+pRPDG*100+1;
00983                      else if(cLPDG < pRPDG) nRPDG=pRPDG*1000+cLPDG*100+1;
00984                      else                   nRPDG=pRPDG*1000+cLPDG*100+3;
00985                      if  (-cRPDG == pL1)    nLPDG=pL2;
00986                      else                   nLPDG=pL1; // -cRPDG == pL2
00987                    }
00988                  }
00989                  break;
00990                case 2: // ....................... pRPDG > 7
00991                  if(cLPDG < 0)
00992                  {
00993                    order=-1;
00994                    if(-cLPDG==pLPDG)
00995                    {
00996                      nLPDG=cRPDG;
00997                      nRPDG=pRPDG;
00998                    }
00999                    else
01000                    {
01001                      if     (cRPDG > pLPDG) nLPDG=cRPDG*1000+pLPDG*100+1;
01002                      else if(cRPDG < pLPDG) nLPDG=pLPDG*1000+cRPDG*100+1;
01003                      else                   nLPDG=pLPDG*1000+cRPDG*100+3;
01004                      if  (-cLPDG == pR1)    nRPDG=pR2;
01005                      else                   nRPDG=pR1; // -cLPDG == pR2
01006                    }
01007                  }
01008                  else // cRPDG < 0
01009                  {
01010                    order= 1;
01011                    if(-cRPDG==pLPDG)
01012                    {
01013                      nLPDG=cLPDG;
01014                      nRPDG=pRPDG;
01015                    }
01016                    else
01017                    {
01018                      if     (cLPDG > pLPDG) nLPDG=cLPDG*1000+pLPDG*100+1;
01019                      else if(cLPDG < pLPDG) nLPDG=pLPDG*1000+cLPDG*100+1;
01020                      else                   nLPDG=pLPDG*1000+cLPDG*100+3;
01021                      if  (-cRPDG == pR1)    nRPDG=pR2;
01022                      else                   nRPDG=pR1; // -cRPDG == pR2
01023                    }
01024                  }
01025                  break;
01026                case 3: // ....................... pLPDG <-7
01027                  if(cLPDG > 0)
01028                  {
01029                    order= 1;
01030                    if(cLPDG==-pRPDG)
01031                    {
01032                      nLPDG=pLPDG;
01033                      nRPDG=cRPDG;
01034                    }
01035                    else
01036                    {
01037                      if     (cRPDG < pRPDG) nRPDG=cRPDG*1000+pRPDG*100-1;
01038                      else if(cRPDG > pRPDG) nRPDG=pRPDG*1000+cRPDG*100-1;
01039                      else                   nRPDG=pRPDG*1000+cRPDG*100-3;
01040                      if  ( cLPDG == pL1)    nLPDG=-pL2;
01041                      else                   nLPDG=-pL1; // cLPDG == pL2
01042                    }
01043                  }
01044                  else // cRPDG > 0
01045                  {
01046                    order=-1;
01047                    if(cRPDG==-pRPDG)
01048                    {
01049                      nLPDG=pLPDG;
01050                      nRPDG=cLPDG;
01051                    }
01052                    else
01053                    {
01054                      if     (cLPDG < pRPDG) nRPDG=cLPDG*1000+pRPDG*100-1;
01055                      else if(cLPDG > pRPDG) nRPDG=pRPDG*1000+cLPDG*100-1;
01056                      else                   nRPDG=pRPDG*1000+cLPDG*100-3;
01057                      if  ( cRPDG == pL1)    nLPDG=-pL2;
01058                      else                   nLPDG=-pL1; // cRPDG == pL2
01059                    }
01060                  }
01061                  break;
01062                case 4: // ....................... pRPDG <-7
01063                  if(cLPDG > 0)
01064                  {
01065                    order=-1;
01066                    if(cLPDG==-pLPDG)
01067                    {
01068                      nLPDG=cRPDG;
01069                      nRPDG=pRPDG;
01070                    }
01071                    else
01072                    {
01073                      if     (cRPDG < pLPDG) nLPDG=cRPDG*1000+pLPDG*100-1;
01074                      else if(cRPDG > pLPDG) nLPDG=pLPDG*1000+cRPDG*100-1;
01075                      else                   nLPDG=pLPDG*1000+cRPDG*100-3;
01076                      if  ( cLPDG == pR1)    nRPDG=-pR2;
01077                      else                   nRPDG=-pR1; // cLPDG == pR2
01078                    }
01079                  }
01080                  else // cRPDG > 0
01081                  {
01082                    order= 1;
01083                    if(cRPDG==-pLPDG)
01084                    {
01085                      nLPDG=cLPDG;
01086                      nRPDG=pRPDG;
01087                    }
01088                    else
01089                    {
01090                      if     (cLPDG < pLPDG) nLPDG=cLPDG*1000+pLPDG*100-1;
01091                      else if(cLPDG > pLPDG) nLPDG=pLPDG*1000+cLPDG*100-1;
01092                      else                   nLPDG=pLPDG*1000+cLPDG*100-3;
01093                      if  ( cRPDG == pR1)    nRPDG=-pR2;
01094                      else                   nRPDG=-pR1; // cRPDG == pR2
01095                    }
01096                  }
01097                  break;
01098              }
01099              break;
01100             case 3: // ------------------------> QDiQ/aQaDiQ + QaQ = QDiQ/aQaDiQ (single)
01101              switch (af)
01102              {
01103                case 1: // ....................... cLPDG > 7
01104                  if(pLPDG < 0)
01105                  {
01106                    order= 1;
01107                    if     (pRPDG > cRPDG) nRPDG=pRPDG*1000+cRPDG*100+1;
01108                    else if(pRPDG < cRPDG) nRPDG=cRPDG*1000+pRPDG*100+1;
01109                    else                   nRPDG=cRPDG*1000+pRPDG*100+3;
01110                    if  (-pLPDG == L1)     nLPDG=L2;
01111                    else                   nLPDG=L1; // -pLPDG == L2
01112                  }
01113                  else // pRPDG < 0
01114                  {
01115                    order=-1;
01116                    if     (pLPDG > cRPDG) nLPDG=pLPDG*1000+cRPDG*100+1;
01117                    else if(pLPDG < cRPDG) nLPDG=cRPDG*1000+pLPDG*100+1;
01118                    else                   nLPDG=cRPDG*1000+pLPDG*100+3;
01119                    if  (-pRPDG == L1)     nRPDG=L2;
01120                    else                   nRPDG=L1; // -pRPDG == L2
01121                  }
01122                  break;
01123                case 2: // ....................... cRPDG > 7
01124                  if(pLPDG < 0)
01125                  {
01126                    order=-1;
01127                    if     (pRPDG > cLPDG) nRPDG=pRPDG*1000+cLPDG*100+1;
01128                    else if(pRPDG < cLPDG) nRPDG=cLPDG*1000+pRPDG*100+1;
01129                    else                   nRPDG=cLPDG*1000+pRPDG*100+3;
01130                    if  (-pLPDG == R1)     nLPDG=R2;
01131                    else                   nLPDG=R1; // -pLPDG == R2
01132                  }
01133                  else // pRPDG < 0
01134                  {
01135                    order= 1;
01136                    if     (pLPDG > cLPDG) nLPDG=pLPDG*1000+cLPDG*100+1;
01137                    else if(pLPDG < cLPDG) nLPDG=cLPDG*1000+pLPDG*100+1;
01138                    else                   nLPDG=cLPDG*1000+pLPDG*100+3;
01139                    if  (-pRPDG == R1)     nRPDG=R2;
01140                    else                   nRPDG=R1; // -pRPDG == R2
01141                  }
01142                  break;
01143                case 3: // ....................... cLPDG <-7 (cRPDG <0)
01144                  if(pLPDG > 0)
01145                  {
01146                    order= 1;
01147                    if     (pRPDG < cRPDG) nRPDG=pRPDG*1000+cRPDG*100-1;
01148                    else if(pRPDG > cRPDG) nRPDG=cRPDG*1000+pRPDG*100-1;
01149                    else                   nRPDG=cRPDG*1000+pRPDG*100-3;
01150                    if  ( pLPDG == L1)     nLPDG=-L2;
01151                    else                   nLPDG=-L1; // pLPDG == L2
01152                  }
01153                  else // pRPDG > 0
01154                  {
01155                    order=-1;
01156                    if     (pLPDG < cRPDG) nLPDG=pLPDG*1000+cRPDG*100-1;
01157                    else if(pLPDG > cRPDG) nLPDG=cRPDG*1000+pLPDG*100-1;
01158                    else                   nLPDG=cRPDG*1000+pLPDG*100-3;
01159                    if  ( pRPDG == L1)     nRPDG=-L2;
01160                    else                   nRPDG=-L1; // pRPDG == L2
01161                  }
01162                  break;
01163                case 4: // ....................... cRPDG <-7 (cLPDG <0)
01164                  if(pLPDG > 0)                       // pRPDG & cLPDG are anti-quarks
01165                  {
01166                    order=-1;
01167                    if     (pRPDG < cLPDG) nRPDG=pRPDG*1000+cLPDG*100-1;
01168                    else if(pRPDG > cLPDG) nRPDG=cLPDG*1000+pRPDG*100-1;
01169                    else                   nRPDG=cLPDG*1000+pRPDG*100-3;
01170                    if  ( pLPDG == R1)     nLPDG=-R2;
01171                    else                   nLPDG=-R1; // pLPDG == R2
01172                  }
01173                  else // pRPDG > 0
01174                  {
01175                    order= 1;
01176                    if     (pLPDG < cLPDG) nLPDG=pLPDG*1000+cLPDG*100-1;
01177                    else if(pLPDG > cLPDG) nLPDG=cLPDG*1000+pLPDG*100-1;
01178                    else                   nLPDG=cLPDG*1000+pLPDG*100-3;
01179                    if  ( pRPDG == R1)     nRPDG=-R2;
01180                    else                   nRPDG=-R1; // pRPDG == R2
01181                  }
01182                  break;
01183              }
01184              break;
01185             case 4: // ------------------------------------> QaQ + aDiQDiQ = QaQ (double)
01186              switch (af)
01187              {
01188                case 1: // ....................... pLPDG > 7 && pRPDG <-7 === LL/RR
01189                  order= 1;
01190                  if(-cLPDG == pL1) nLPDG= pL2;
01191                  else              nLPDG= pL1;
01192                  if( cRPDG == pR1) nRPDG=-pR2;
01193                  else              nRPDG=-pR1;
01194                  break;
01195                case 2: // ...................... pLPDG > 7 && pRPDG <-7 === LR/RL
01196                  order=-1;
01197                  if(-cRPDG == pL1) nLPDG= pL2;
01198                  else              nLPDG= pL1;
01199                  if( cLPDG == pR1) nRPDG=-pR2;
01200                  else              nRPDG=-pR1;
01201                  break;
01202                case 3: // ...................... pRPDG > 7 && pLPDG <-7 === LL/RR
01203                  order= 1;
01204                  if( cLPDG == pL1) nLPDG=-pL2;
01205                  else              nLPDG=-pL1;
01206                  if(-cRPDG == pR1) nRPDG= pR2;
01207                  else              nRPDG= pR1;
01208                  break;
01209                case 4: // ...................... pRPDG > 7 && pLPDG <-7 === LR/RL
01210                  order=-1;
01211                  if( cRPDG == pL1) nLPDG=-pL2;
01212                  else              nLPDG=-pL1;
01213                  if(-cLPDG == pR1) nRPDG= pR2;
01214                  else              nRPDG= pR1;
01215                  break;
01216              }
01217              break;
01218             case 5: // ------------------------------------> aDiQDiQ + QaQ = QaQ (double)
01219              switch (af)
01220              {
01221                case 1: // ...................... cLPDG > 7 && cRPDG <-7 === LL/RR
01222                  order= 1;
01223                  if(-pLPDG == L1) nLPDG= L2;
01224                  else             nLPDG= L1;
01225                  if( pRPDG == R1) nRPDG=-R2;
01226                  else             nRPDG=-R1;
01227                  break;
01228                case 2: // ...................... cLPDG > 7 && cRPDG <-7 === LR/RL
01229                  order=-1;
01230                  if(-pRPDG == L1) nRPDG= L2;
01231                  else             nRPDG= L1;
01232                  if( pLPDG == R1) nLPDG=-R2;
01233                  else             nLPDG=-R1;
01234                  break;
01235                case 3: // ...................... cRPDG > 7 && cLPDG <-7 === LL/RR
01236                  order= 1;
01237                  if( pLPDG == L1) nLPDG=-L2;
01238                  else             nLPDG=-L1;
01239                  if(-pRPDG == R1) nRPDG= R2;
01240                  else             nRPDG= R1;
01241                  break;
01242                case 4: // ...................... cRPDG > 7 && cLPDG <-7 === LR/RL
01243                  order=-1;
01244                  if( pRPDG == L1) nRPDG=-L2;
01245                  else             nRPDG=-L1;
01246                  if(-pLPDG == R1) nLPDG= R2;
01247                  else             nLPDG= R1;
01248                  break;
01249              }
01250              break;
01251             case 6: // ------------------------------------> QDiQ + aQaDiQ = QaQ (double)
01252              switch (af)
01253              {
01254                case 1:
01255                  order=-1;
01256                  if(-cRPDG == pL1) nLPDG= pL2;
01257                  else              nLPDG= pL1;
01258                  if( pRPDG ==  L1) nRPDG= -L2;
01259                  else              nRPDG= -L1;
01260                  break;
01261                case 2:
01262                  order= 1;
01263                  if(-cLPDG == pL1) nLPDG= pL2;
01264                  else              nLPDG= pL1;
01265                  if( pRPDG ==  R1) nRPDG= -R2;
01266                  else              nRPDG= -R1;
01267                  break;
01268                case 3:
01269                  order= 1;
01270                  if(-cRPDG == pR1) nRPDG= pR2;
01271                  else              nRPDG= pR1;
01272                  if( pLPDG ==  L1) nLPDG= -L2;
01273                  else              nLPDG= -L1;
01274                  break;
01275                case 4:
01276                  order=-1;
01277                  if(-cLPDG == pR1) nRPDG= pR2;
01278                  else              nRPDG= pR1;
01279                  if( pLPDG ==  R1) nLPDG= -R2;
01280                  else              nLPDG= -R1;
01281                  break;
01282                case 5:
01283                  order=-1;
01284                  if(-pRPDG ==  L1) nRPDG=  L2;
01285                  else              nRPDG=  L1;
01286                  if( cRPDG == pL1) nLPDG=-pL2;
01287                  else              nLPDG=-pL1;
01288                  break;
01289                case 6:
01290                  order= 1;
01291                  if(-pLPDG ==  L1) nLPDG=  L2;
01292                  else              nLPDG=  L1;
01293                  if( cRPDG == pR1) nRPDG=-pR2;
01294                  else              nRPDG=-pR1;
01295                  break;
01296                case 7:
01297                  order= 1;
01298                  if(-pRPDG ==  R1) nRPDG=  R2;
01299                  else              nRPDG=  R1;
01300                  if( cLPDG == pL1) nLPDG=-pL2;
01301                  else              nLPDG=-pL1;
01302                  break;
01303                case 8:
01304                  order=-1;
01305                  if(-pLPDG ==  R1) nLPDG=  R2;
01306                  else              nLPDG=  R1;
01307                  if( cLPDG == pR1) nRPDG=-pR2;
01308                  else              nRPDG=-pR1;
01309                  break;
01310              }
01311              break;
01312           }
01313           if(!order) G4cerr<<"-Warning-G4QInel::Constr: t="<<tf<<", a="<<af<<", cL="<<cLPDG
01314                            <<", cR="<<cRPDG<<", pL="<<pLPDG<<", pR="<<pRPDG<<G4endl;
01315           else
01316           {
01317             // With theNewHypotheticalPartons the min mass must be calculated & compared
01318             G4int LT=1;
01319             if(std::abs(nLPDG) > 7) ++LT;
01320             G4int RT=1;
01321             if(std::abs(nRPDG) > 7) ++RT;
01322             G4double minM=0.;
01323             G4bool sing=true;
01324             if(cLT==2 && cRT==2)
01325             {
01326               aLPDG=0;
01327               aRPDG=0;
01328               if(cLPDG>0)
01329               {
01330                 aLPDG=nLPDG/100;
01331                 aRPDG=(-nRPDG)/100;
01332               }
01333               else //cRPDG>0
01334               {
01335                 aRPDG=nRPDG/100;
01336                 aLPDG=(-nLPDG)/100;
01337               }
01338               G4int nL1=aLPDG/10;
01339               G4int nL2=aLPDG%10;
01340               G4int nR1=aRPDG/10;
01341               G4int nR2=aRPDG%10;
01342               if(nL1!=nR1 && nL1!=nR2 && nL2!=nR1 && nL2!=nR2) // Unreducable DiQ-aDiQ
01343               {
01344 #ifdef debug
01345                 G4cout<<"G4QIonIonCollis::Constr:aLPDG="<<aLPDG<<", aRPDG="<<aRPDG<<G4endl;
01346 #endif
01347                 sing=false;
01348                 G4QPDGCode tmp;
01349                 std::pair<G4int,G4int> pB=tmp.MakeTwoBaryons(nL1, nL2, nR1, nR2);
01350                 minM=G4QPDGCode(pB.first).GetMass()+G4QPDGCode(pB.second).GetMass();
01351               }
01352             }
01353             if(sing)
01354             {
01355               std::pair<G4int,G4int> newPair = std::make_pair(nLPDG,nRPDG);
01356               G4QContent newStQC(newPair);        // NewString QuarkContent
01357 #ifdef debug
01358               G4cout<<"G4QIn::Con: LPDG="<<nLPDG<<",RPDG="<<nRPDG<<",QC="<<newStQC<<G4endl;
01359 #endif
01360               G4int minPDG=newStQC.GetSPDGCode(); // PDG of the Lightest Hadron=String
01361               minM=G4QPDGCode(minPDG).GetMass() + mPi0; // Min SingleHadron=String Mass
01362             }
01363             // Compare this mass
01364             G4bool win=false;
01365             G4double    pSM=0.;
01366             if(pSM2>0.) pSM=std::sqrt(pSM2);
01367             if(minC && pSM2 > maxiM2)             // Up to now any positive mass is good
01368             {
01369               maxiM2=pSM2;
01370               win=true;
01371             }
01372             else if(!minC || pSM > minM)
01373             {
01374               pExcess=pSM-minM;
01375               if(minC || pExcess > excess)
01376               {
01377                 minC=false;
01378                 excess=pExcess;
01379                 win=true;
01380               }
01381             }
01382             if(win)
01383             {
01384               sst=pst;
01385               sLPDG=nLPDG;
01386               sRPDG=nRPDG;
01387               sOrd=order;
01388             }
01389           } // End of IF(new partons are created)
01390         } // End of IF(compatible partons)
01391         //} // End of positive squared mass of the fused string
01392       } // End of the LOOP over the possible partners (with the exclusive if for itself)
01393       if(sOrd)                                       // The best pStringCandidate was found
01394       {
01395         G4LorentzVector cL4M=cLeft->Get4Momentum();
01396         G4LorentzVector cR4M=cRight->Get4Momentum();
01397         G4QParton* pLeft=(*sst)->GetLeftParton();
01398         G4QParton* pRight=(*sst)->GetRightParton();
01399         G4LorentzVector pL4M=pLeft->Get4Momentum();
01400         G4LorentzVector pR4M=pRight->Get4Momentum();
01401 #ifdef debug
01402         G4cout<<"G4QIonIonCollis::Constr:cS4M="<<cS4M<<" fused/w pS4M="<<pL4M+pR4M<<G4endl;
01403 #endif
01404         if(sOrd>0)
01405         {
01406           pL4M+=cL4M;
01407           pR4M+=cR4M;
01408         }
01409         else
01410         {
01411           pL4M+=cR4M;
01412           pR4M+=cL4M;
01413         }
01414         pLeft->SetPDGCode(sLPDG);
01415         pLeft->Set4Momentum(pL4M);
01416         pRight->SetPDGCode(sRPDG);
01417         pRight->Set4Momentum(pR4M);
01418         delete (*ist);
01419         strings.erase(ist);
01420 #ifdef debug
01421         G4LorentzVector ss4M=pL4M+pR4M;
01422         G4cout<<"G4QIonIonCollision::Constr:Created,4M="<<ss4M<<",m2="<<ss4M.m2()<<G4endl;
01423 #endif
01424         if( ist != strings.begin() ) // To avoid going below begin() (for Windows)
01425         {
01426           ist--;
01427           con=false;
01428 #ifdef debug
01429           G4cout<<"G4QIonIonCollision::Construct: *** IST Decremented ***"<<G4endl;
01430 #endif
01431         }
01432         else
01433         {
01434           con=true;
01435 #ifdef debug
01436           G4cout<<"G4QIonIonCollision::Construct: *** IST Begin ***"<<G4endl;
01437 #endif
01438           break;
01439         }
01440       } // End of the IF(the best partnerString candidate was found)
01441       else
01442       {
01443 #ifdef debug
01444         G4cout<<"-Warning-G4QInel::Const: S4M="<<cS4M<<",M2="<<cSM2<<" Leave ASIS"<<G4endl;
01445 #endif
01446         ++problem;
01447         con=false;
01448       }
01449     }
01450     else con=false;
01451    } // End of loop over ist iterator
01452 #ifdef debug
01453    G4cout<<"G4QIonIonCollision::Construct: *** IST While *** , con="<<con<<G4endl;
01454 #endif
01455   } // End of "con" while 
01456 #ifdef edebug
01457   // This print has meaning only if something appear between it and the StringFragmLOOP
01458   G4LorentzVector t4M=theProjNucleus.Get4Momentum()+theTargNucleus.Get4Momentum();//NucInLS
01459   G4int rC=totChg-theProjNucleus.GetZ()-theTargNucleus.GetZ();
01460   G4int rB=totBaN-theProjNucleus.GetA()-theTargNucleus.GetA();
01461   G4int nStr=strings.size();
01462   G4cout<<"-EMCLS-G4QIn::Const: AfterSUPPRESION #ofS="<<nStr<<",tNuc4M(E=M)="<<sum<<G4endl;
01463   for(G4int i=0; i<nStr; i++)
01464   {
01465     G4LorentzVector strI4M=strings[i]->Get4Momentum();
01466     t4M+=strI4M;
01467     G4int sChg=strings[i]->GetCharge();
01468     rC-=sChg;
01469     G4int sBaN=strings[i]->GetBaryonNumber();
01470     rB-=sBaN;
01471     G4cout<<"-EMCLS-G4QIonIonCollision::Construct: St#"<<i<<", 4M="<<strI4M<<", M="
01472           <<strI4M.m()<<", C="<<sChg<<", B="<<sBaN<<G4endl;
01473   }
01474   G4cout<<"-EMCLS-G4QInel::Construct: r4M="<<t4M-totLS4M<<", rC="<<rC<<", rB="<<rB<<G4endl;
01475 #endif
01476   //
01477   // --- If a problem is foreseen then the DiQaDiQ strings should be reduced if possible --
01478   //
01479 #ifdef debug
01480     G4cout<<"G4QIonIonCollision::Construct: problem="<<problem<<G4endl;
01481 #endif
01482   if(problem)
01483   {
01484     G4int nOfStr=strings.size();
01485 #ifdef debug
01486     G4cout<<"G4QIonIonCollision::Constr:SecurityDiQaDiQReduction, #OfStr="<<nOfStr<<G4endl;
01487 #endif
01488     for (G4int astring=0; astring < nOfStr; astring++)
01489     {
01490       G4QString* curString=strings[astring];
01491       G4QParton* cLeft=curString->GetLeftParton();
01492       G4QParton* cRight=curString->GetRightParton();
01493       G4int LT=cLeft->GetType();
01494       G4int RT=cRight->GetType();
01495       G4int sPDG=cLeft->GetPDGCode();
01496       G4int nPDG=cRight->GetPDGCode();
01497       if(LT==2 && RT==2)
01498       {
01499 #ifdef debug
01500         G4cout<<"G4QIonIonCollis::Constr:TrySelfReduString,L="<<sPDG<<",R="<<nPDG<<G4endl;
01501 #endif
01502         if( cLeft->ReduceDiQADiQ(cLeft, cRight) ) // DiQ-aDiQ pair was successfully reduced
01503         {
01504           sPDG=cLeft->GetPDGCode();
01505           nPDG=cRight->GetPDGCode();
01506 #ifdef debug
01507           G4cout<<"+G4QInel::Const:#"<<astring<<" Reduced, L="<<sPDG<<", R="<<nPDG<<G4endl;
01508 #endif
01509         }
01510 #ifdef debug
01511         else G4cout<<"--*--G4QInel::Const: #"<<astring<<" DQ-aDQ reduction Failed"<<G4endl;
01512 #endif
01513       } // End of the found DiQ/aDiQ pair
01514       else if(sPDG==3 && nPDG==-3)
01515       {
01516         sPDG= 1;
01517         nPDG=-1;
01518         cLeft->SetPDGCode(sPDG);
01519         cRight->SetPDGCode(nPDG);
01520       }
01521       else if(sPDG==-3 && nPDG==3)
01522       {
01523         sPDG=-1;
01524         nPDG= 1;
01525         cLeft->SetPDGCode(sPDG);
01526         cRight->SetPDGCode(nPDG);
01527       }
01528     }
01529     SwapPartons();
01530   } // End of IF(problem)
01531 #ifdef edebug
01532   G4LorentzVector u4M=theProjNucleus.Get4Momentum()+theTargNucleus.Get4Momentum();//NucInLS
01533   G4int rCh=totChg-theProjNucleus.GetZ()-theTargNucleus.GetZ();
01534   G4int rBa=totBaN-theProjNucleus.GetA()-theTargNucleus.GetA();
01535   G4int nStri=strings.size();
01536   G4cout<<"-EMCLS-G4QIn::Const: FinalConstruct, #ofSt="<<nStri<<",tN4M(E=M)="<<t4M<<G4endl;
01537   for(G4int i=0; i<nStri; i++)
01538   {
01539     G4LorentzVector strI4M=strings[i]->Get4Momentum();
01540     u4M+=strI4M;
01541     G4int sChg=strings[i]->GetCharge();
01542     rCh-=sChg;
01543     G4int sBaN=strings[i]->GetBaryonNumber();
01544     rBa-=sBaN;
01545     G4cout<<"-EMCLS-G4QIonIonCollision::Construct: St#"<<i<<", 4M="<<strI4M<<", M="
01546           <<strI4M.m()<<", C="<<sChg<<", B="<<sBaN<<G4endl;
01547   }
01548   G4cout<<"-EMCLS-G4QInel::Construct: r4M="<<u4M-totLS4M<<",rC="<<rCh<<",rB="<<rBa<<G4endl;
01549 #endif
01550 }

G4QIonIonCollision::~G4QIonIonCollision (  ) 

Definition at line 1552 of file G4QIonIonCollision.cc.

01553 {
01554   std::for_each(strings.begin(), strings.end(), DeleteQString() );
01555 }


Member Function Documentation

G4int G4QIonIonCollision::AnnihilationOrder ( G4int  LS,
G4int  MS,
G4int  uP,
G4int  mP,
G4int  sP,
G4int  nP 
) [protected]

Definition at line 3671 of file G4QIonIonCollision.cc.

References G4cerr, G4cout, and G4endl.

Referenced by Breeder().

03673 {//                                             ^L^^ Curent ^^R^
03674   G4int Ord=0;
03675   //       Curent   Partner
03676   if      (LS==2 && MPS==2 )                 // Fuse 2 Q-aQ strings to DiQ-aDiQ
03677   {
03678 #ifdef debug
03679     G4cout<<"G4QIonIonCollision::AnnihilationOrder:QaQ/QaQ->DiQ-aDiQ, uPDG="<<uPDG<<G4endl;
03680 #endif
03681     if     ( (uPDG>0 && sPDG>0 && mPDG<0 && nPDG<0) || 
03682              (uPDG<0 && sPDG<0 && mPDG>0 && nPDG>0) ) Ord= 1; // LL/RR
03683     else if( (uPDG>0 && nPDG>0 && mPDG<0 && sPDG<0) || 
03684              (uPDG<0 && nPDG<0 && mPDG>0 && sPDG>0) ) Ord=-1; // LR/RL
03685     else G4cerr<<"-Warning-G4QIonIonCollision::AnnihilationOrder: Wrong 22 fusion, L="
03686                <<uPDG<<",R="<<mPDG<<",cL="<<sPDG<<",cR="<<nPDG<<G4endl;
03687   }
03688   else if ( LS==2 && MPS==3 )
03689   {
03690     if     (uPDG > 7)                                // Fuse pLDiQ
03691     {
03692 #ifdef debug
03693       G4cout<<"G4QIonIonCollision::AnnihOrder:pLDiQ, sPDG="<<sPDG<<", nPDG="<<nPDG<<G4endl;
03694 #endif
03695       if     ( sPDG<0 && (-sPDG==uPDG/1000 || -sPDG==(uPDG/100)%10) ) Ord= 1; // LL/RR
03696       else if( nPDG<0 && (-nPDG==uPDG/1000 || -nPDG==(uPDG/100)%10) ) Ord=-1; // LR/RL
03697       else G4cerr<<"-Warning-G4QIonIonCollision::AnnihilationOrder: Wrong pLDiQ, L="<<uPDG
03698                  <<",R="<<mPDG<<",cL="<<sPDG<<",cR="<<nPDG<<G4endl;
03699     }
03700     else if (mPDG > 7)                               // Fuse pRDiQ
03701     {
03702 #ifdef debug
03703       G4cout<<"G4QIonIonCollision::AnnihOrder:pRDiQ, sPDG="<<sPDG<<", nPDG="<<nPDG<<G4endl;
03704 #endif
03705       if     ( sPDG<0 && (-sPDG==mPDG/1000 || -sPDG==(mPDG/100)%10) ) Ord=-1; // LR/RL
03706       else if( nPDG<0 && (-nPDG==mPDG/1000 || -nPDG==(mPDG/100)%10) ) Ord= 1; // LL/RR
03707       else G4cerr<<"-Warning-G4QIonIonCollision::AnnihilationOrder: Wrong pRDiQ, L="<<uPDG
03708                  <<",R="<<mPDG<<",cL="<<sPDG<<",cR="<<nPDG<<G4endl;
03709     }
03710     else if (uPDG <-7)                               // Fuse pLaDiQ
03711     {
03712 #ifdef debug
03713       G4cout<<"G4QIonIonCollision::AnnihOrder:pLaDiQ, sPDG="<<sPDG<<",nPDG="<<nPDG<<G4endl;
03714 #endif
03715       if     ( sPDG>0 && (sPDG==(-uPDG)/1000 || sPDG==((-uPDG)/100)%10) ) Ord= 1; // LL/RR
03716       else if( nPDG>0 && (nPDG==(-uPDG)/1000 || nPDG==((-uPDG)/100)%10) ) Ord=-1; // LR/RL
03717       else G4cerr<<"-Warning-G4QIonIonCollision::AnnihilationOrder: Wrong pLaDiQ, L="<<uPDG
03718                  <<", R="<<mPDG<<", cL="<<sPDG<<", cR="<<nPDG<<G4endl;
03719     }
03720     else if (mPDG <-7)                              // Fuse pRaDiQ
03721     {
03722 #ifdef debug
03723       G4cout<<"G4QIonIonCollision::AnnihOrder:pRaDiQ, sPDG="<<sPDG<<",nPDG="<<nPDG<<G4endl;
03724 #endif
03725       if     ( sPDG>0 && (sPDG==(-mPDG)/1000 || sPDG==((-mPDG)/100)%10) ) Ord=-1; // LR/RL
03726       else if( nPDG>0 && (nPDG==(-mPDG)/1000 || nPDG==((-mPDG)/100)%10) ) Ord= 1; // LL/RR
03727       else G4cerr<<"-Warning-G4QIonIonCollision::AnnihilationOrder: Wrong pRaDiQ, L="<<uPDG
03728                  <<", R="<<mPDG<<", cL="<<sPDG<<", cR="<<nPDG<<G4endl;
03729     }
03730     else if( (sPDG<0 && (-sPDG==mPDG || -sPDG==uPDG) ) ||
03731              (nPDG<0 && (-nPDG==mPDG || -nPDG==uPDG) ) ) Ord= 2; // @@ Annihilation fusion
03732 #ifdef debug
03733     else G4cout<<"-Warning-G4QIonIonCollision::AnnihilatOrder:Wrong23StringFusion"<<G4endl;
03734     G4cout<<"G4QIonIonCollision::AnnihilationOrder: Ord="<<Ord<<",sPDG="<<sPDG<<",nPDG="
03735           <<nPDG<<", uPDG="<<uPDG<<",mPDG="<<mPDG<<G4endl;
03736 #endif
03737   }
03738   else if ( LS==3 && MPS==2 )
03739   {
03740     if     (sPDG > 7)                                // Fuse cLDiQ
03741     {
03742 #ifdef debug
03743       G4cout<<"G4QIonIonCollision::AnnihOrder:cLDiQ, uPDG="<<uPDG<<", mPDG="<<mPDG<<G4endl;
03744 #endif
03745       if     ( uPDG<0 && (-uPDG==sPDG/1000 || -uPDG==(sPDG/100)%10) ) Ord= 1; // LL/RR
03746       else if( mPDG<0 && (-mPDG==sPDG/1000 || -mPDG==(sPDG/100)%10) ) Ord=-1; // LR/RL
03747       else G4cerr<<"-Warning-G4QIonIonCollision::AnnihilationOrder: Wrong cLDiQ, L="<<uPDG
03748                  <<",R="<<mPDG<<",cL="<<sPDG<<",cR="<<nPDG<<G4endl;
03749     }
03750     else if (nPDG > 7)                               // Fuse cRDiQ
03751     {
03752 #ifdef debug
03753       G4cout<<"G4QIonIonCollision::AnnihOrder:cRDiQ, uPDG="<<uPDG<<", mPDG="<<mPDG<<G4endl;
03754 #endif
03755       if     ( uPDG<0 && (-uPDG==nPDG/1000 || -uPDG==(nPDG/100)%10) ) Ord=-1; // LR/RL
03756       else if( mPDG<0 && (-mPDG==nPDG/1000 || -mPDG==(nPDG/100)%10) ) Ord= 1; // LL/RR
03757       else G4cerr<<"-Warning-G4QIonIonCollision::AnnihilationOrder: Wrong cRDiQ, L="<<uPDG
03758                  <<",R="<<mPDG<<",cL="<<sPDG<<",cR="<<nPDG<<G4endl;
03759     }
03760     else if (sPDG <-7)                               // Fuse cLaDiQ
03761     {
03762 #ifdef debug
03763       G4cout<<"G4QIonIonCollision::AnnihOrder:cLaDiQ, uPDG="<<uPDG<<",mPDG="<<mPDG<<G4endl;
03764 #endif
03765       if     ( uPDG>0 && (uPDG==(-sPDG)/1000 || uPDG==((-sPDG)/100)%10) ) Ord= 1; // LL/RR
03766       else if( mPDG>0 && (mPDG==(-sPDG)/1000 || mPDG==((-sPDG)/100)%10) ) Ord=-1; // LR/RL
03767       else G4cerr<<"-Warning-G4QIonIonCollision::AnnihilationOrder: Wrong cLaDiQ, L="<<uPDG
03768                  <<", R="<<mPDG<<", cL="<<sPDG<<", cR="<<nPDG<<G4endl;
03769     }
03770     else if (nPDG <-7)                              // Fuse cRaDiQ
03771     {
03772 #ifdef debug
03773       G4cout<<"G4QIonIonCollision::AnnihOrder:cRaDiQ, uPDG="<<uPDG<<",mPDG="<<mPDG<<G4endl;
03774 #endif
03775       if     ( uPDG>0 && (uPDG==(-nPDG)/1000 || uPDG==((-nPDG)/100)%10) ) Ord=-1; // LR/RL
03776       else if( mPDG>0 && (mPDG==(-nPDG)/1000 || mPDG==((-nPDG)/100)%10) ) Ord= 1; // LL/RR
03777       else G4cerr<<"-Warning-G4QIonIonCollision::AnnihilationOrder: Wrong cRaDiQ, L="<<uPDG
03778                  <<", R="<<mPDG<<", cL="<<sPDG<<", cR="<<nPDG<<G4endl;
03779     }
03780     else if( (uPDG<0 && (-uPDG==sPDG || -uPDG==nPDG) ) ||
03781              (mPDG<0 && (-mPDG==sPDG || -mPDG==nPDG) ) ) Ord=2; // @@ Annihilation fusion
03782 #ifdef debug
03783     else G4cout<<"-Warning-G4QIonIonCollision::AnnihilatOrder:Wrong32StringFusion"<<G4endl;
03784     G4cout<<"G4QIonIonCollision::AnnihilationOrder: Ord="<<Ord<<",sPDG="<<sPDG<<",nPDG="
03785           <<nPDG<<", uPDG="<<uPDG<<",mPDG="<<mPDG<<G4endl;
03786 #endif
03787   }
03788   else if ( (LS==2 && MPS==4) || (LS==4 && MPS==2) )
03789   {
03790     if     (uPDG > 7)  // Annihilate pLDiQ
03791     {
03792 #ifdef debug
03793       G4cout<<"G4QIonIonCollision::AnnihilOrder:pLDiQ,sPDG="<<sPDG<<",nPDG="<<nPDG<<G4endl;
03794 #endif
03795       if     ( sPDG<0 && (-sPDG==uPDG/1000 || -sPDG==(uPDG/100)%10) &&
03796                (nPDG==(-mPDG)/1000 || nPDG==((-mPDG)/100)%10) ) Ord= 1; // LL/RR
03797       else if( nPDG<0 && (-nPDG==uPDG/1000 || -nPDG==(uPDG/100)%10) &&
03798                (sPDG==(-mPDG)/1000 || sPDG==((-mPDG)/100)%10) ) Ord=-1; // LR/RL
03799       else G4cerr<<"-Warning-G4QIonIonCollision::AnnihilationOrder: Wrong 24 pLDiQ, L="
03800                  <<uPDG<<",R="<<mPDG<<",cL="<<sPDG<<",cR="<<nPDG<<G4endl;
03801     }
03802     else if (mPDG >7) // Annihilate pRDiQ
03803     {
03804 #ifdef debug
03805       G4cout<<"G4QIonIonCollision::AnnihilOrder:PRDiQ,sPDG="<<sPDG<<",nPDG="<<nPDG<<G4endl;
03806 #endif
03807       if     ( sPDG<0 && (-sPDG==mPDG/1000 || -sPDG==(mPDG/100)%10) &&
03808                (nPDG==(-uPDG)/1000 || nPDG==((-uPDG)/100)%10) ) Ord=-1; // LR/RL
03809       else if( nPDG<0 && (-nPDG==mPDG/1000 || -nPDG==(mPDG/100)%10) &&
03810                (sPDG==(-uPDG)/1000 || sPDG==((-uPDG)/100)%10) ) Ord= 1; // LL/RR
03811       else G4cerr<<"-Warning-G4QIonIonCollision::AnnihilationOrder: Wrong 24 pLDiQ, L="
03812                  <<uPDG<<",R="<<mPDG<<",cL="<<sPDG<<",cR="<<nPDG<<G4endl;
03813     }
03814     else if (sPDG > 7)   // Annihilate cLDiQ
03815     {
03816 #ifdef debug
03817       G4cout<<"G4QIonIonCollision::AnnihilOrder:cLDiQ,uPDG="<<uPDG<<",mPDG="<<mPDG<<G4endl;
03818 #endif
03819       if     ( uPDG<0 && (-uPDG==sPDG/1000 || -uPDG==(sPDG/100)%10) &&
03820                (mPDG==(-nPDG)/1000 || mPDG==((-nPDG)/100)%10) ) Ord= 1; // LL/RR
03821       else if( mPDG<0 && (-mPDG==sPDG/1000 || -mPDG==(sPDG/100)%10) &&
03822                (uPDG==(-nPDG)/1000 || uPDG==((-nPDG)/100)%10) ) Ord=-1; // LR/RL
03823       else G4cerr<<"-Warning-G4QIonIonCollision::AnnihilationOrder: Wrong 24 cLDiQ, L="
03824                  <<uPDG<<",R="<<mPDG<<",cL="<<sPDG<<",cR="<<nPDG<<G4endl;
03825     }
03826     else if (nPDG > 7)   // Annihilate cRDiQ
03827     {
03828 #ifdef debug
03829       G4cout<<"G4QIonIonCollision::AnnihilOrder:cRDiQ,uPDG="<<uPDG<<",mPDG="<<mPDG<<G4endl;
03830 #endif
03831       if     ( uPDG<0 && (-uPDG==nPDG/1000 || -uPDG==(nPDG/100)%10) &&
03832                (mPDG==(-sPDG)/1000 || mPDG==((-sPDG)/100)%10) ) Ord=-1; // LR/RL
03833       else if( mPDG<0 && (-mPDG==nPDG/1000 || -mPDG==(nPDG/100)%10) &&
03834                (uPDG==(-sPDG)/1000 || uPDG==((-sPDG)/100)%10) ) Ord= 1; // LL/RR
03835       else G4cerr<<"-Warning-G4QIonIonCollision::AnnihilationOrder: Wrong 24 cRDiQ, L="
03836                  <<uPDG<<",R="<<mPDG<<",cL="<<sPDG<<",cR="<<nPDG<<G4endl;
03837     }
03838 #ifdef debug
03839     else G4cout<<"-Warning-G4QIonIonCollision::AnnihilOrder: Wrong 24StringFusion"<<G4endl;
03840     G4cout<<"G4QIonIonCollision::AnnihilationOrder: Ord="<<Ord<<",sPDG="<<sPDG<<",nPDG="
03841           <<nPDG<<", uPDG="<<uPDG<<",mPDG="<<mPDG<<G4endl;
03842 #endif
03843   }
03844   else if ( LS==3 && MPS==3 )
03845   {
03846     if     (uPDG > 7)  // Annihilate pLDiQ
03847     {
03848 #ifdef debug
03849       G4cout<<"G4QIonIonCollision::AnnihilOrder:pLDiQ,sPDG="<<sPDG<<",nPDG="<<nPDG<<G4endl;
03850 #endif
03851       if     ( sPDG<-7 && (-nPDG==uPDG/1000 || -nPDG==(uPDG/100)%10) &&
03852                (mPDG==(-sPDG)/1000 || mPDG==((-sPDG)/100)%10) ) Ord=-1; // LR/RL
03853       else if( nPDG<-7 && (-sPDG==uPDG/1000 || -sPDG==(uPDG/100)%10) &&
03854                (mPDG==(-nPDG)/1000 || mPDG==((-nPDG)/100)%10) ) Ord= 1; // LL/RR
03855       else G4cerr<<"-Warning-G4QIonIonCollision::AnnihilationOrder: Wrong 33 pLDiQ, L="
03856                  <<uPDG<<",R="<<mPDG<<",cL="<<sPDG<<",cR="<<nPDG<<G4endl;
03857     }
03858     else if(mPDG > 7)  // Annihilate pRDiQ
03859     {
03860 #ifdef debug
03861       G4cout<<"G4QIonIonCollision::AnnihilOrder:pRDiQ,sPDG="<<sPDG<<",nPDG="<<nPDG<<G4endl;
03862 #endif
03863       if     ( sPDG<-7 && (-nPDG==mPDG/1000 || -nPDG==(mPDG/100)%10) &&
03864                (uPDG==(-sPDG)/1000 || uPDG==((-sPDG)/100)%10) ) Ord= 1; // LL/RR
03865       else if( nPDG<-7 && (-sPDG==mPDG/1000 || -sPDG==(mPDG/100)%10) &&
03866                (uPDG==(-nPDG)/1000 || uPDG==((-nPDG)/100)%10) ) Ord=-1; // LR/RL
03867       else G4cerr<<"-Warning-G4QIonIonCollision::AnnihilationOrder: Wrong33pRDiQ, L="<<uPDG
03868                  <<",R="<<mPDG<<",cL="<<sPDG<<",cR="<<nPDG<<G4endl;
03869     }
03870     else if(sPDG > 7)  // Annihilate cLDiQ
03871     {
03872 #ifdef debug
03873       G4cout<<"G4QIonIonCollision::AnnihilOrder:cLDiQ,uPDG="<<uPDG<<",mPDG="<<mPDG<<G4endl;
03874 #endif
03875       if     ( uPDG<-7 && (-mPDG==sPDG/1000 || -mPDG==(sPDG/100)%10) &&
03876                (nPDG==(-uPDG)/1000 || nPDG==((-uPDG)/100)%10) ) Ord=-1; // LR/RL
03877       else if( mPDG<-7 && (-uPDG==sPDG/1000 || -uPDG==(sPDG/100)%10) &&
03878                (nPDG==(-mPDG)/1000 || nPDG==((-mPDG)/100)%10) ) Ord= 1; // LL/RR
03879       else G4cerr<<"-Warning-G4QIonIonCollision::AnnihilationOrder: Wrong 33 cLDiQ, L="
03880                  <<uPDG<<",R="<<mPDG<<",cL="<<sPDG<<",cR="<<nPDG<<G4endl;
03881     }
03882     else if(nPDG > 7)  // Annihilate cRDiQ
03883     {
03884 #ifdef debug
03885       G4cout<<"G4QIonIonCollision::AnnihilOrder:cRDiQ,uPDG="<<uPDG<<",mPDG="<<mPDG<<G4endl;
03886 #endif
03887       if     ( uPDG<-7 && (-mPDG==nPDG/1000 || -mPDG==(nPDG/100)%10) &&
03888                (nPDG==(-uPDG)/1000 || nPDG==((-uPDG)/100)%10) ) Ord= 1; // LL/RR
03889       else if( mPDG<-7 && (-uPDG==nPDG/1000 || -sPDG==(nPDG/100)%10) &&
03890                (sPDG==(-mPDG)/1000 || sPDG==((-mPDG)/100)%10) ) Ord=-1; // LR/RL
03891       else G4cerr<<"-Warning-G4QIonIonCollision::AnnihilationOrder: Wrong 33 cRDiQ, L="
03892                  <<uPDG<<",R="<<mPDG<<",cL="<<sPDG<<",cR="<<nPDG<<G4endl;
03893     }
03894 #ifdef debug
03895     else G4cout<<"-Warning-G4QIonIonCollision::AnnihilOrder: Wrong 33StringFusion"<<G4endl;
03896     G4cout<<"G4QIonIonCollision::AnnihilationOrder: Ord="<<Ord<<",sPDG="<<sPDG<<",nPDG="
03897           <<nPDG<<", uPDG="<<uPDG<<",mPDG="<<mPDG<<G4endl;
03898 #endif
03899   }
03900   return Ord;
03901 }

void G4QIonIonCollision::Breeder (  )  [protected]

Definition at line 1791 of file G4QIonIonCollision.cc.

References G4QContent::AddParton(), AnnihilationOrder(), DBL_MAX, G4Quasmon::DecayQHadron(), edebug, FatalException, G4QString::FragmentString(), G4cerr, G4cout, G4endl, G4Exception(), G4QParton::Get4Momentum(), G4QString::Get4Momentum(), G4QHadron::Get4Momentum(), G4QNucleus::GetA(), G4QHadron::GetBaryonNumber(), G4QString::GetBaryonNumber(), G4QHadron::GetCharge(), G4QString::GetCharge(), G4QString::GetDirection(), G4QString::GetLeftParton(), G4QPDGCode::GetMass(), G4QNucleus::GetPDG(), G4QPDGCode::GetPDGCode(), G4QHadron::GetPDGCode(), G4QParton::GetPDGCode(), G4QString::GetQC(), G4QChipolino::GetQPDG1(), G4QChipolino::GetQPDG2(), G4QString::GetRightParton(), G4QContent::GetSPDGCode(), G4QParton::GetType(), G4QNucleus::GetZ(), LT, G4QParton::ReduceDiQADiQ(), ReducePair(), G4QHadron::Set4Momentum(), G4QParton::Set4Momentum(), G4QParton::SetPDGCode(), sqr(), and SumPartonPDG().

Referenced by Fragment().

01792 { // This is the member function, which returns the resulting vector of Hadrons & Quasmons
01793   static const G4double  eps = 0.001;                              // Tolerance in MeV
01794   //
01795   // ------------ At this point the strings are fragmenting to hadrons in LS -------------
01796   //
01797 #ifdef edebug
01798   G4LorentzVector totLS4M=theProjNucleus.Get4Momentum()+theTargNucleus.Get4Momentum(); //LS
01799   G4int totChg=theProjNucleus.GetZ()+theTargNucleus.GetZ();
01800   G4int totBaN=theProjNucleus.GetA()+theTargNucleus.GetA();
01801   G4int nStri=strings.size();
01802   G4cout<<"-EMCLS-G4QIn::Breed: CHECKRecovery #ofS="<<nStri<<",N4M(E=M)="<<totLS4M<<G4endl;
01803   for(G4int i=0; i<nStri; i++)
01804   {
01805     G4LorentzVector strI4M=strings[i]->Get4Momentum();
01806     totLS4M+=strI4M;
01807     G4int sChg=strings[i]->GetCharge();
01808     totChg+=sChg;
01809     G4int sBaN=strings[i]->GetBaryonNumber();
01810     totBaN+=sBaN;
01811     G4cout<<"-EMCLS-G4QIonIonCollision::Breeder: St#"<<i<<", 4M="<<strI4M<<", M="
01812           <<strI4M.m()<<", C="<<sChg<<", B="<<sBaN<<G4endl;
01813   }
01814 #endif
01815   G4int nOfStr=strings.size();
01816 #ifdef debug
01817   G4cout<<"G4QIonIonCollision::Breeder: BeforeFragmentation, #OfStr="<<nOfStr<<G4endl;
01818 #endif
01819   G4LorentzVector ft4M(0.,0.,0.,0.);
01820   G4QContent      ftQC(0,0,0,0,0,0);
01821   G4bool          ftBad=false;
01822   for(G4int i=0; i < nOfStr; ++i)
01823   {
01824     G4LorentzVector pS4M=strings[i]->Get4Momentum(); // String 4-momentum
01825     ft4M+=pS4M;
01826     G4QContent pSQC=strings[i]->GetQC();             // String Quark Content
01827     ftQC+=pSQC;
01828     if(pS4M.m2() < 0.) ftBad=true;
01829 #ifdef debug
01830     G4cout<<">G4QIonIonCollision::Breed:1stTest,S#"<<i<<",4M="<<pS4M<<",QC="<<pSQC<<G4endl;
01831 #endif
01832   }
01833   if(ftBad)
01834   {
01835     G4Quasmon* stringQuasmon = new G4Quasmon(ftQC, ft4M);
01836 #ifdef debug
01837     G4cout<<"->G4QIonIonCollision::Breed:*TotQ*,QC="<<ftQC<<",4M="<<ft4M<<ft4M.m()<<G4endl;
01838 #endif
01839     theQuasmons.push_back(stringQuasmon);
01840     G4LorentzVector rp4M=theProjNucleus.Get4Momentum(); // Nucleus 4-momentum in LS
01841     G4int rpPDG=theProjNucleus.GetPDG();
01842     G4QHadron* resPNuc = new G4QHadron(rpPDG,rp4M);
01843     theResult->push_back(resPNuc);                  // Fill the residual projectile nucleus
01844     G4LorentzVector rt4M=theTargNucleus.Get4Momentum(); // Nucleus 4-momentum in LS
01845     G4int rtPDG=theTargNucleus.GetPDG();
01846     G4QHadron* resTNuc = new G4QHadron(rtPDG,rt4M);
01847     theResult->push_back(resTNuc);                  // Fill the residual target nucleus
01848     return;
01849   }
01850   for (G4int astring=0; astring < nOfStr; astring++)
01851   {
01852 #ifdef edebug
01853     G4LorentzVector sum=theProjNucleus.Get4Momentum()+theTargNucleus.Get4Momentum(); //inLS
01854     G4int rChg=totChg-theProjNucleus.GetZ()-theTargNucleus.GetZ();
01855     G4int rBaN=totBaN-theProjNucleus.GetA()-theTargNucleus.GetA();
01856     G4int nOfHadr=theResult->size();
01857     G4cout<<"-EMCLS-G4QIonIonCollision::Breed:#ofSt="<<nOfStr<<",#ofHad="<<nOfHadr<<G4endl;
01858     for(G4int i=astring; i<nOfStr; i++)
01859     {
01860       G4LorentzVector strI4M=strings[i]->Get4Momentum();
01861       sum+=strI4M;
01862       G4int sChg=strings[i]->GetCharge();
01863       rChg-=sChg;
01864       G4int sBaN=strings[i]->GetBaryonNumber();
01865       rBaN-=sBaN;
01866       G4cout<<"-EMCLS-G4QI::Breed:S#"<<i<<",4M="<<strI4M<<",C="<<sChg<<",B="<<sBaN<<G4endl;
01867     }
01868     for(G4int i=0; i<nOfHadr; i++)
01869     {
01870       G4LorentzVector hI4M=(*theResult)[i]->Get4Momentum();
01871       sum+=hI4M;
01872       G4int hChg=(*theResult)[i]->GetCharge();
01873       rChg-=hChg;
01874       G4int hBaN=(*theResult)[i]->GetBaryonNumber();
01875       rBaN-=hBaN;
01876       G4cout<<"-EMCLS-G4QIn::Breed: H#"<<i<<",4M="<<hI4M<<",C="<<hChg<<",B="<<hBaN<<G4endl;
01877     }
01878     G4cout<<"....-EMCLS-G4QInel::Br:r4M="<<sum-totLS4M<<",rC="<<rChg<<",rB="<<rBaN<<G4endl;
01879 #endif
01880     G4QString* curString=strings[astring];
01881     if(!curString->GetDirection()) continue;  // Historic for the dead strings: DoesNotWork
01882 #ifdef edebug
01883     G4int curStrChg = curString->GetCharge();
01884     G4int curStrBaN = curString->GetBaryonNumber();
01885 #endif
01886     G4LorentzVector curString4M = curString->Get4Momentum();
01887 #ifdef debug
01888     G4cout<<"=--=>G4QIonIonCollision::Breeder: String#"<<astring<<",s4M/m="<<curString4M
01889           <<curString4M.m()<<", LPart="<<curString->GetLeftParton()->GetPDGCode()
01890           <<", RPart="<<curString->GetRightParton()->GetPDGCode()<<G4endl;
01891 #endif
01892     G4QHadronVector* theHadrons = 0;           // Prototype of theStringFragmentationOUTPUT
01893     theHadrons=curString->FragmentString(true);// !! Fragmenting the String !!
01894     if (!theHadrons)                           // The string can not be fragmented
01895     {
01896       // First try to correct the diQ-antiDiQ strings, converting them to Q-antiQ
01897       G4QParton* cLeft=curString->GetLeftParton();
01898       G4QParton* cRight=curString->GetRightParton();
01899       G4int sPDG=cLeft->GetPDGCode();
01900       G4int nPDG=cRight->GetPDGCode();
01901       G4int LT=cLeft->GetType();
01902       G4int RT=cRight->GetType();
01903       G4int LS=LT+RT;
01904       if(LT==2 && RT==2)
01905       {
01906 #ifdef debug
01907         G4cout<<"G4QIonIonCollision::Breed:TryReduceString, L="<<sPDG<<",R="<<nPDG<<G4endl;
01908 #endif
01909         if( cLeft->ReduceDiQADiQ(cLeft, cRight) ) // DiQ-aDiQ pair was successfully reduced
01910         {
01911           LT=1;
01912           RT=1;
01913           LS=2;
01914           sPDG=cLeft->GetPDGCode();
01915           nPDG=cRight->GetPDGCode();
01916 #ifdef debug
01917           G4cout<<"G4QIonIonCollision::Breed:AfterReduction,L="<<sPDG<<",R="<<nPDG<<G4endl;
01918 #endif
01919           theHadrons=curString->FragmentString(true);
01920           cLeft=curString->GetLeftParton();
01921           cRight=curString->GetRightParton();
01922 #ifdef debug
01923           G4cout<<"G4QInel::Breed:L="<<cLeft->Get4Momentum()<<",R="<<cRight->Get4Momentum()
01924                 <<G4endl;
01925 #endif
01926         }
01927 #ifdef debug
01928         else G4cout<<"^G4QIonIonCollision::Breed: DQ-aDQ reduction to Q-aQ Failed"<<G4endl;
01929 #endif
01930       } // End of the SelfReduction
01931 #ifdef debug
01932       G4cout<<"G4QIonIonCollision::Breed:AfterRedAttempt, theH="<<theHadrons<<", L4M="
01933             <<cLeft->Get4Momentum()<<", R4M="<<cRight->Get4Momentum()<<G4endl;
01934 #endif
01935       unsigned next=astring+1;                 // The next string position
01936       if (!theHadrons)                         // The string can not be fragmented
01937       {
01938         G4int fusionDONE=0; // StringFusion didn't happen (1=Fuse L+L/R+R, -1=Fuse L+R/R+L)
01939         if(next < strings.size())             // TheString isn't theLastString can fuse
01940         {
01941           G4int fustr=0;                       // The found partner index (never can be 0)
01942           G4int swap=0;                        // working interger for swapping parton PDG
01943           G4double Vmin=DBL_MAX;               // Prototype of the found Velocity Distance 
01944           G4int dPDG=nPDG;
01945           G4int qPDG=sPDG;
01946           if(dPDG<-99 || (dPDG>0&&dPDG<7) || qPDG>99 || (qPDG<0 && qPDG>-7))
01947           {
01948             swap=qPDG;
01949             qPDG=dPDG;
01950             dPDG=swap;
01951           }
01952           if(dPDG>99) dPDG/=100;
01953           if(qPDG<-99) qPDG=-(-qPDG)/100;
01954 #ifdef debug
01955           G4cout<<"G4QIonIonCollision::Breed:TryFuseStringS, q="<<qPDG<<", a="<<dPDG
01956                 <<", n="<<next<<G4endl;
01957 #endif
01958           G4ThreeVector curV=curString4M.vect()/curString4M.e();
01959           G4int reduce=0;                      // a#of reduced Q-aQ pairs
01960           G4int restr=0;                       // To use beyon the LOOP for printing
01961           G4int MPS=0;                         // PLS for the selected string
01962           for (restr=next; restr < nOfStr; restr++)
01963           {
01964             reduce=0;
01965             G4QString* reString=strings[restr];
01966             G4QParton* Left=reString->GetLeftParton();
01967             G4QParton* Right=reString->GetRightParton();
01968             G4int uPDG=Left->GetPDGCode();
01969             G4int mPDG=Right->GetPDGCode();
01970             G4int PLT =Left->GetType();
01971             G4int PRT =Right->GetType();
01972             G4int aPDG=mPDG;
01973             G4int rPDG=uPDG;
01974             if(aPDG<-99 || (aPDG>0 && aPDG<7) || rPDG>99 || (rPDG<0 && rPDG>-7))
01975             {
01976               swap=rPDG;
01977               rPDG=aPDG;
01978               aPDG=swap;
01979             }
01980             if(aPDG > 99) aPDG/=100;
01981             if(rPDG <-99) rPDG=-(-rPDG)/100;
01982             // Try to reduce two DQ-aDQ strings
01983 #ifdef debug
01984             G4cout<<"G4Qnel::Breed: TryReduce #"<<restr<<", q="<<rPDG<<",a="<<aPDG<<G4endl;
01985 #endif
01986             if(LT==2 && RT==2 && PLT==2 && PRT==2)    // Have a chance for the reduction
01987             {
01988               G4int cQ1=(-qPDG)/10;
01989               G4int cQ2=(-qPDG)%10;
01990               G4int cA1=dPDG/10;
01991               G4int cA2=dPDG%10;
01992               G4int pQ1=(-rPDG)/10;
01993               G4int pQ2=(-rPDG)%10;
01994               G4int pA1=aPDG/10;
01995               G4int pA2=aPDG%10;
01996 #ifdef debug
01997                   G4cout<<"G4QIonIonCollision::Breeder: cQ="<<cQ1<<","<<cQ2<<", cA="<<cA1<<","
01998                     <<cA2<<", pQ="<<pQ1<<","<<pQ2<<", pA="<<pA1<<","<<pA2<<G4endl;
01999 #endif
02000               G4bool iQA = (cA1==pQ1 || cA1==pQ2 || cA2==pQ1 || cA2==pQ2);
02001               G4bool iAQ = (cQ1==pA1 || cQ1==pA2 || cQ2==pA1 || cQ2==pA2);
02002               if(iQA) reduce++;
02003               if(iAQ) reduce++;
02004               if  (reduce==2)                  // Two quark pairs can be reduced
02005               {
02006                 if(sPDG>0 && uPDG<0)           // LL/RR Reduction
02007                 {
02008                   std::pair<G4int,G4int> resLL=ReducePair(sPDG/100, (-uPDG)/100);
02009                   G4int newCL=resLL.first;
02010                   G4int newPL=resLL.second;
02011                   if(!newCL || !newPL)
02012                   {
02013                     G4cerr<<"*G4QIonIonCollision::Breed:CL="<<newCL<<",PL="<<newPL<<G4endl;
02014                     G4Exception("G4QIonIonCollision::Breeder:","72",FatalException,"2LL-");
02015                   }
02016                   std::pair<G4int,G4int> resRR=ReducePair((-nPDG)/100, mPDG/100);
02017                   G4int newCR=resRR.first;
02018                   G4int newPR=resRR.second;
02019                   if(!newCR || !newPR)
02020                   {
02021                     G4cerr<<"*G4QIonIonCollision::Breed:CR="<<newCR<<",PR="<<newPR<<G4endl;
02022                     G4Exception("G4QIonIonCollision::Breeder:","72",FatalException,"2RR-");
02023                   }
02024                   cLeft->SetPDGCode(newCL);    // Reset the left quark of curString
02025                   cRight->SetPDGCode(-newCR);  // Reset the right quark of curString
02026                   Left->SetPDGCode(-newPL);    // Reset the left quark of reString
02027                   Right->SetPDGCode(newPR);    // Reset the right quark of reString
02028                   break;                       // Break out of the reString internal LOOP
02029                 }
02030                 else if(sPDG<0 && uPDG>0)           // LL/RR Reduction
02031                 {
02032                   std::pair<G4int,G4int> resLL=ReducePair((-sPDG)/100, uPDG/100);
02033                   G4int newCL=resLL.first;
02034                   G4int newPL=resLL.second;
02035                   if(!newCL || !newPL)
02036                   {
02037                     G4cerr<<"*G4QIonIonCollision::Breed:CL="<<newCL<<",PL="<<newPL<<G4endl;
02038                     G4Exception("G4QIonIonCollision::Breeder:","72",FatalException,"2LL+");
02039                   }
02040                   std::pair<G4int,G4int> resRR=ReducePair(nPDG/100, (-mPDG)/100);
02041                   G4int newCR=resRR.first;
02042                   G4int newPR=resRR.second;
02043                   if(!newCR || !newPR)
02044                   {
02045                     G4cerr<<"*G4QIonIonCollision::Breed:CR="<<newCR<<",PR="<<newPR<<G4endl;
02046                     G4Exception("G4QIonIonCollision::Breeder:","72",FatalException,"2RR+");
02047                   }
02048                   cLeft->SetPDGCode(-newCL);   // Reset the left quark of curString
02049                   cRight->SetPDGCode(newCR);   // Reset the right quark of curString
02050                   Left->SetPDGCode(newPL);     // Reset the left quark of reString
02051                   Right->SetPDGCode(-newPR);   // Reset the right quark of reString
02052                   break;                       // Break out of the reString internal LOOP
02053                 }
02054                 else if(sPDG>0 && mPDG<0)                // LR Reduction
02055                 {
02056                   std::pair<G4int,G4int> resLL=ReducePair(sPDG/100, (-mPDG)/100);
02057                   G4int newCL=resLL.first;
02058                   G4int newPR=resLL.second;
02059                   if(!newCL || !newPR)
02060                   {
02061                     G4cerr<<"*G4QIonIonCollision::Breed:CL="<<newCL<<",PR="<<newPR<<G4endl;
02062                     G4Exception("G4QIonIonCollision::Breeder:","72",FatalException,"2-LR");
02063                   }
02064                   std::pair<G4int,G4int> resRR=ReducePair((-nPDG)/100, uPDG/100);
02065                   G4int newCR=resRR.first;
02066                   G4int newPL=resRR.second;
02067                   if(!newCR || !newPR)
02068                   {
02069                     G4cerr<<"*G4QIonIonCollision::Breed:CR="<<newCR<<",PL="<<newPL<<G4endl;
02070                     G4Exception("G4QIonIonCollision::Breeder:","72",FatalException,"2-LR");
02071                   }
02072                   cLeft->SetPDGCode(newCL);    // Reset the left quark of curString
02073                   cRight->SetPDGCode(-newCR);  // Reset the right quark of curString
02074                   Left->SetPDGCode(newPL);     // Reset the left quark of reString
02075                   Right->SetPDGCode(-newPR);   // Reset the right quark of reString
02076                   break;                       // Break out of the reString internal LOOP
02077                 }
02078                 else                           // RL Reduction
02079                 {
02080                   std::pair<G4int,G4int> resLL=ReducePair((-sPDG)/100, mPDG/100);
02081                   G4int newCL=resLL.first;
02082                   G4int newPR=resLL.second;
02083                   if(!newCL || !newPR)
02084                   {
02085                     G4cerr<<"*G4QIonIonCollision::Breed:CL="<<newCL<<",PR="<<newPR<<G4endl;
02086                     G4Exception("G4QIonIonCollision::Breeder:","72",FatalException,"2-RL");
02087                   }
02088                   std::pair<G4int,G4int> resRR=ReducePair(nPDG/100, (-uPDG)/100);
02089                   G4int newCR=resRR.first;
02090                   G4int newPL=resRR.second;
02091                   if(!newCR || !newPR)
02092                   {
02093                     G4cerr<<"*G4QIonIonCollision::Breed:CR="<<newCR<<",PL="<<newPL<<G4endl;
02094                     G4Exception("G4QIonIonCollision::Breeder:","72",FatalException,"2-RL");
02095                   }
02096                   cLeft->SetPDGCode(-newCL);   // Reset the left quark of curString
02097                   cRight->SetPDGCode(newCR);   // Reset the right quark of curString
02098                   Left->SetPDGCode(-newPL);    // Reset the left quark of reString
02099                   Right->SetPDGCode(newPR);    // Reset the right quark of reString
02100                   break;                       // Break out of the reString internal LOOP
02101                 }
02102               } // End of IF(possible reduction)
02103             }
02104             // Check StringsCanBeFused: all QaQ+QaQ(22), single QaQ+QDiQ/aQaDtQ(23/32),
02105             //                          double QaQ+DiQaDiQ(24/42), QDiQ+aDiQaQ(34/43)
02106 #ifdef debug
02107             G4cout<<"G4QInel::Breed: TryFuse/w #"<<restr<<",q="<<rPDG<<",a="<<aPDG<<G4endl;
02108 #endif
02109             G4int PLS=PLT+PRT;
02110             if( (LS==2 && PLS==2) ||                           // QaQ+QaQ always to DiQaDiQ
02111                 ( ( (LS==2 && PLS==3) || (LS==3 && PLS==2) ) &&// QaQ w QDiQ/aQaDiQ(single)
02112                   ( (aPDG> 7 && (-dPDG==aPDG/10   || -dPDG==aPDG%10) )   || // cAQ w DQ
02113                     (dPDG> 7 && (-aPDG==dPDG/10   || -aPDG==dPDG%10) )   || // AQ w cDQ
02114                     (rPDG<-7 && (qPDG==(-rPDG)/10 || qPDG==(-rPDG)%10) ) || // cQ w ADQ
02115                     (qPDG<-7 && (rPDG==(-qPDG)/10 || rPDG==(-qPDG)%10) )    // Q w cADQ
02116                     //|| (aPDG< 0 && -aPDG==qPDG) || (dPDG< 0 && -dPDG==rPDG) // aQ w Q 
02117                   )
02118                 )                 ||                           // -----------------------
02119                 ( ( LS==2 && PLS==4                          &&// QaQ w DiQaDiQ (double)
02120                     (aPDG> 7 && (-dPDG == aPDG/10 || -dPDG == aPDG%10) ) &&
02121                     (rPDG<-7 && (qPDG==(-rPDG)/10 || qPDG==(-rPDG)%10) )
02122                   )       ||
02123                   ( LS==4 && PLS==2                          &&// DiQaDiQ w QaQ (double)
02124                     (dPDG> 7 && (-aPDG == dPDG/10 || -aPDG == dPDG%10) ) &&
02125                     (qPDG<-7 && (rPDG==(-qPDG)/10 || rPDG==(-qPDG)%10) )
02126                   )
02127                 )                 ||                           // -----------------------
02128                 ( LS==3 && PLS==3                            &&// QDiQ w aDiQaQ (double)
02129                   ( (aPDG> 7 && (-dPDG == aPDG/10 || -dPDG == aPDG%10) &&
02130                      qPDG<-7 && (rPDG==(-qPDG)/10 || rPDG==(-qPDG)%10)
02131                     )       ||
02132                     (dPDG> 7 && (-aPDG == dPDG/10 || -aPDG == dPDG%10) &&
02133                      rPDG<-7 && (dPDG==(-rPDG)/10 || dPDG==(-rPDG)%10) 
02134                     )
02135                   )
02136                 )
02137               )
02138             {
02139               G4LorentzVector reString4M = reString->Get4Momentum();
02140               G4ThreeVector reV = reString4M.vect()/reString4M.e();
02141               G4double dV=(curV-reV).mag2();   // Squared difference between relVelocities
02142 #ifdef debug
02143               G4cout<<"G4QIonIonCollision::Breeder: StringCand#"<<restr<<", q="<<rPDG
02144                     <<", a="<<aPDG<<", L="<<uPDG<<", R="<<mPDG<<",dV="<<dV<<G4endl;
02145 #endif
02146               if(dV < Vmin)
02147               {
02148                 Vmin=dV;
02149                 fustr=restr;
02150                 MPS=PLS;
02151               }
02152             }
02153           }
02154           if(reduce==2)                        // String mutual reduction happened
02155           {
02156 #ifdef debug
02157             G4cout<<"-G4QIonIonCollision::Breed:Reduced #"<<astring<<" & #"<<restr<<G4endl;
02158 #endif
02159             astring--;                         // String was QCreduced using another String
02160             continue;                          // Repeat fragmentation of the same string
02161           }
02162           if(fustr)                            // The partner was found -> fuse strings
02163           {
02164 #ifdef debug
02165             G4cout<<"G4QInel::Breeder: StPartner#"<<fustr<<", LT="<<LT<<",RT="<<RT<<G4endl;
02166 #endif
02167             G4QString* fuString=strings[fustr];
02168             G4QParton* Left=fuString->GetLeftParton();
02169             G4QParton* Right=fuString->GetRightParton();
02170             G4int uPDG=Left->GetPDGCode();
02171             G4int mPDG=Right->GetPDGCode();
02172             G4int rPDG=uPDG;
02173             G4int aPDG=mPDG;
02174             if(aPDG<-99 || (aPDG>0 && aPDG<7) || rPDG>99 || (rPDG<0 && rPDG>-7))
02175             {
02176               swap=rPDG;
02177               rPDG=aPDG;
02178               aPDG=swap;
02179             }
02180             if(aPDG > 99) aPDG/=100;
02181             if(rPDG <-99) rPDG=-(-rPDG)/100;
02182             // Check that the strings can fuse (no anti-diquarks assumed)
02183             G4LorentzVector resL4M(0.,0.,0.,0.);
02184             G4LorentzVector resR4M(0.,0.,0.,0.);
02185             G4LorentzVector L4M=cLeft->Get4Momentum();
02186             G4LorentzVector R4M=cRight->Get4Momentum();
02187 #ifdef debug
02188             G4cout<<"G4QIonIonCollision::Breeder:BeforeFuDir,sL="<<sPDG<<",nR="<<nPDG
02189                   <<",uL="<<uPDG<<",mR="<<mPDG<<",L4M="<<L4M<<",R4M="<<R4M<<G4endl;
02190 #endif
02191             G4LorentzVector PL4M=Left->Get4Momentum();
02192             G4LorentzVector PR4M=Right->Get4Momentum();
02193             fusionDONE=AnnihilationOrder(LS, MPS, uPDG, mPDG, sPDG, nPDG);
02194             if     (fusionDONE>0)
02195             {
02196               if(fusionDONE>1)                             // Anihilation algorithm
02197               {
02198                 if     ( (uPDG<0 || nPDG<0) && -uPDG==nPDG ) Left->SetPDGCode(sPDG);
02199                 else if( (mPDG<0 || sPDG<0) && -mPDG==sPDG ) Right->SetPDGCode(nPDG);
02200                 else if( (uPDG<0 || sPDG<0) && -uPDG==sPDG ) Left->SetPDGCode(nPDG);
02201                 else if( (mPDG<0 || nPDG<0) && -mPDG==nPDG ) Right->SetPDGCode(sPDG);
02202               }
02203               {
02204                 Left->SetPDGCode(SumPartonPDG(uPDG, sPDG));
02205                 Right->SetPDGCode(SumPartonPDG(mPDG, nPDG));
02206               }
02207               Left->Set4Momentum(L4M+PL4M);
02208               Right->Set4Momentum(R4M+PR4M);
02209 #ifdef debug
02210               G4cout<<"G4QIonIonCollision::Breeder:LL/RR s4M="<<fuString->Get4Momentum()
02211                     <<",S="<<L4M+PL4M+R4M+PR4M<<", L="<<Left->Get4Momentum()<<", R="
02212                     <<Right->Get4Momentum()<<G4endl;
02213 #endif
02214             }
02215             else if(fusionDONE<0)
02216             {
02217               Left->SetPDGCode(SumPartonPDG(uPDG, nPDG));
02218               Left->Set4Momentum(L4M+PR4M);
02219               Right->SetPDGCode(SumPartonPDG(mPDG, sPDG));
02220               Right->Set4Momentum(R4M+PL4M);
02221 #ifdef debug
02222               G4cout<<"G4QIonIonCollision::Breeder:LR/RL s4M="<<fuString->Get4Momentum()
02223                     <<",S="<<L4M+PL4M+R4M+PR4M<<", L="<<Left->Get4Momentum()<<", R="
02224                     <<Right->Get4Momentum()<<G4endl;
02225 #endif
02226             }
02227 #ifdef debug
02228             else G4cout<<"-Warning-G4QIonIonCollision::Breeder: WrongStringFusion"<<G4endl;
02229 #endif
02230 #ifdef edebug
02231             G4cout<<"#EMC#G4QIonIonCollision::Breed:StringFused,F="<<fusionDONE<<",L="<<L4M
02232                   <<",R="<<R4M<<",pL="<<PL4M<<",pR="<<PR4M<<",nL="<<Left->Get4Momentum()
02233                   <<",nR="<<Right->Get4Momentum()<<",S="<<fuString->Get4Momentum()<<G4endl;
02234 #endif
02235             if(fusionDONE)
02236             {
02237 #ifdef debug
02238               G4cout<<"###G4QIonIonCollision::Breed: Str#"<<astring<<" fused/w Str#"<<fustr
02239                     <<"->"<<fuString->Get4Momentum()<<fuString->Get4Momentum().m()<<G4endl;
02240 #endif
02241               continue;                          // @@ killing of the curString is excluded
02242             }
02243           }
02244 #ifdef debug
02245           else
02246           {
02247 
02248             G4cerr<<"**G4QIonIonCollision::Breed:*NoPart*M="<<curString->Get4Momentum().m()
02249                   <<", F="<<fusionDONE<<", LPDG="<<curString->GetLeftParton()->GetPDGCode()
02250                   <<", RPDG="<<curString->GetRightParton()->GetPDGCode()<<G4endl;
02251           }
02252 #endif
02253         }
02254         if(!fusionDONE)                          // Fusion of theString failed, try hadrons
02255         {
02256           G4int nHadr=theResult->size();         // #of collected resulting hadrons upToNow
02257 #ifdef debug
02258           G4cout<<"+++4QFragmentation::Breeder:*TryHadr* M="<<curString->Get4Momentum().m()
02259                 <<", nH="<<nHadr<<", LPDG="<<curString->GetLeftParton()->GetPDGCode()
02260                 <<", RPDG="<<curString->GetRightParton()->GetPDGCode()<<G4endl;
02261 #endif
02262           // The only solution now is to try fusion with the secondary hadron
02263           while( (nHadr=theResult->size()) && !theHadrons)
02264           {
02265 #ifdef edebug
02266             for(G4int i=0; i<nHadr; i++)
02267             {
02268               G4LorentzVector h4M=(*theResult)[i]->Get4Momentum();
02269               G4int hPDG=(*theResult)[i]->GetPDGCode();
02270               G4QContent hQC=(*theResult)[i]->GetQC();
02271               G4cout<<"-EMC-G4QInel::Breed:H#"<<i<<",4M="<<h4M<<hQC<<",PDG="<<hPDG<<G4endl;
02272             }
02273 #endif
02274             G4int fusDONE=0;                      // String+Hadron fusion didn't happen
02275             G4int fuhad=-1;                       // The found hadron index
02276             G4int newPDG=0;                       // PDG ofTheParton afterMerging with Hadr
02277             G4int secPDG=0;                       // Second PDG oParton afterMerging w/Hadr
02278             G4double maM2=-DBL_MAX;               // Prototype of the max ResultingStringM2
02279             G4LorentzVector selH4M(0.,0.,0.,0.);  // 4-mom of the selected hadron
02280             G4QHadron* selHP=0;                   // Pointer to the used hadron for erasing
02281             G4QString* cString=strings[astring];  // Must be the last string by definition
02282             G4LorentzVector cString4M = cString->Get4Momentum();
02283             cLeft=cString->GetLeftParton();
02284             cRight=cString->GetRightParton();
02285             G4int sumT=cLeft->GetType()+cRight->GetType();
02286             sPDG=cLeft->GetPDGCode();
02287             nPDG=cRight->GetPDGCode();
02288             G4int cMax=0;                         // Both or only one cuark can merge
02289             for (G4int reh=0; reh < nHadr; reh++)
02290             {
02291               G4QHadron* curHadr=(*theResult)[reh];
02292               G4QContent curQC=curHadr->GetQC();  // Quark content of the hadron
02293               G4int partPDG1=curQC.AddParton(sPDG);
02294               G4int partPDG2=curQC.AddParton(nPDG);
02295 #ifdef debug
02296               G4cout<<"G4QIonIonCollision::Breeder: Hadron # "<<reh<<", QC="<<curQC
02297                       <<", P1="<<partPDG1<<", P2="<<partPDG2<<G4endl;
02298 #endif
02299               if(partPDG1 || partPDG2)            // Hadron can merge at least w/one parton
02300               {
02301                 G4int cCur=1;
02302                 if(sumT>3 && partPDG1 && partPDG2) cCur=2;
02303                 G4LorentzVector curHadr4M = curHadr->Get4Momentum();
02304                 G4double M2=(cString4M+curHadr4M).m2();// SquaredMass of theResultingString
02305 #ifdef debug
02306                 G4cout<<"G4QIonIonCollision::Breeder:*IN*Hadron#"<<reh<<",M2="<<M2<<G4endl;
02307 #endif
02308                 if( (sumT<4 || cCur>=cMax) && M2 > maM2) // DeepAnnihilation for DiQ-aDiQ 
02309                 {
02310                   maM2=M2;
02311                   fuhad=reh;
02312                   if(partPDG1)
02313                   {
02314                     fusDONE= 1;
02315                     newPDG=partPDG1;
02316                     if(partPDG2)
02317                     {
02318                       secPDG=partPDG2;
02319                       cMax=cCur;
02320                     }
02321                   }
02322                   else
02323                   {
02324                     fusDONE=-1;
02325                     newPDG=partPDG2;
02326                   }
02327 #ifdef debug
02328                   G4cout<<"G4QInel::Br:*Selected*,P1="<<partPDG1<<",P2="<<partPDG2<<G4endl;
02329 #endif
02330                   selH4M=curHadr4M;
02331                   selHP=curHadr;
02332                 } // End of IF(update selection)
02333               } // End of IF(HadronCanMergeWithTheString)
02334             } // End of the LOOP over Hadrons
02335 #ifdef debug
02336             G4cout<<"G4QIonIonCollision::Breeder: fuh="<<fuhad<<",fus="<<fusDONE<<G4endl;
02337 #endif
02338             if(fuhad>-1)                          // The partner was found -> fuse H&S
02339             {
02340               if     (fusDONE>0)                  // Fuse hadron with the left parton
02341               {
02342                 cLeft->SetPDGCode(newPDG);
02343                 G4LorentzVector newLeft=cLeft->Get4Momentum()+selH4M;
02344                 cLeft->Set4Momentum(newLeft);
02345                 if(secPDG && cMax>1)
02346                 {
02347 #ifdef debug
02348                   G4cout<<"G4QInel::Br:TryReduce, nPDG="<<newPDG<<",sPDG="<<secPDG<<G4endl;
02349 #endif
02350                   cLeft->ReduceDiQADiQ(cLeft, cRight);
02351                 }
02352 #ifdef debug
02353                 G4cout<<"G4QIonIonCollision::Breed: Left, s4M="<<curString->Get4Momentum()
02354                       <<", L4M="<<newLeft<<", R4M="<<cRight->Get4Momentum()<<", h4M="
02355                       <<selH4M<<", newL="<<newPDG<<", oldR="<<cRight->GetPDGCode()<<G4endl;
02356 #endif
02357               }
02358               else if(fusDONE<0)                  // Fuse hadron with the right parton
02359               {
02360                 cRight->SetPDGCode(newPDG);
02361                 G4LorentzVector newRight=cRight->Get4Momentum()+selH4M;
02362                 cRight->Set4Momentum(newRight);
02363 #ifdef debug
02364                 G4cout<<"G4QIonIonCollision::Breed: Right, s4M="<<curString->Get4Momentum()
02365                       <<", L4M="<<cLeft->Get4Momentum()<<", R4M="<<newRight<<", h4M="
02366                       <<selH4M<<", newR="<<newPDG<<", oldL="<<cLeft->GetPDGCode()<<G4endl;
02367 #endif
02368               }
02369 #ifdef debug
02370               else G4cout<<"-G4QIonIonCollision::Breed: Wrong String+HadronFusion"<<G4endl;
02371 #endif
02372 #ifdef debug
02373               if(fusDONE) G4cout<<"####G4QIonIonCollision::Breeder: String #"<<astring
02374                                 <<" is fused with Hadron #"<<fuhad
02375                                 <<", new4Mom="<<curString->Get4Momentum()
02376                                 <<", M2="<<curString->Get4Momentum().m2()
02377                                 <<", QC="<<curString->GetQC()<<G4endl;
02378 #endif
02379             }
02380             else
02381             {
02382 #ifdef debug
02383               G4cerr<<"**G4QIonIonCollision::Breed:*NoH*,M="<<curString->Get4Momentum().m()
02384                     <<", LPDG="<<curString->GetLeftParton()->GetPDGCode()
02385                     <<", RPDG="<<curString->GetRightParton()->GetPDGCode()<<G4endl;
02386               // @@ Temporary exception for the future solution
02387               //G4Exception("G4QIonIonCollision::Breed:","72",FatalException,"SHNotFused");
02388 #endif
02389               break;                           // Breake the While LOOP
02390             } // End of the namespace where both Fusion and reduction have failed
02391             // The fused hadron must be excluded from theResult
02392 #ifdef debug
02393             G4cout<<"G4QIonIonCollision::Breed: before HR, nH="<<theResult->size()<<G4endl;
02394             G4int icon=0;                              // Loop counter
02395 #endif
02396             G4QHadronVector::iterator ih;
02397 #ifdef debug
02398             G4bool found=false;
02399 #endif
02400             for(ih = theResult->begin(); ih != theResult->end(); ih++)
02401             {
02402 #ifdef debug
02403               G4cout<<"G4QInelast::Breeder:#"<<icon<<", i="<<(*ih)<<", sH="<<selHP<<G4endl;
02404               icon++;
02405 #endif
02406               if((*ih)==selHP)
02407               {
02408 #ifdef debug
02409                 G4cout<<"G4QInel::Breed: *HadronFound*, PDG="<<selHP->GetPDGCode()<<G4endl;
02410 #endif
02411                 G4LorentzVector p4M=selHP->Get4Momentum();
02412                 curString4M+=p4M;
02413 #ifdef edebug
02414                 G4int Chg=selHP->GetCharge();
02415                 G4int BaN=selHP->GetBaryonNumber();
02416                 curStrChg+=Chg;
02417                 curStrBaN+=BaN;
02418                 G4cout<<"-EMC->>>G4QIonIonCollision::Breed: S+=H, 4M="<<curString4M<<", M="
02419                       <<curString4M.m()<<", Charge="<<curStrChg<<", B="<<curStrBaN<<G4endl;
02420 #endif
02421                 delete selHP;                          // delete the Hadron
02422                 theResult->erase(ih);                  // erase the Hadron from theResult
02423 #ifdef debug
02424                 found=true;
02425 #endif
02426                 break;                                 // beak the LOOP over hadrons
02427               }
02428             } // End of the LOOP over hadrons
02429 #ifdef debug
02430             if(!found) G4cout<<"*G4QIonIonCollision::Breed:nH="<<theResult->size()<<G4endl;
02431 #endif
02432             // New attempt of the string decay
02433             theHadrons=curString->FragmentString(true);
02434 #ifdef debug
02435             G4cout<<"G4QInel::Breeder: tH="<<theHadrons<<",nH="<<theResult->size()<<G4endl;
02436 #endif
02437           } // End of the while LOOP over the fusion with hadrons
02438 #ifdef debug
02439           G4cout<<"*G4QIonIonCollision::Breed: *CanTryToDecay?* nH="<<theHadrons<<", next="
02440                 <<next<<" =? nS="<<strings.size()<<", nR="<<theResult->size()<<G4endl;
02441 #endif
02442           if(!theHadrons && next == strings.size() && !(theResult->size()))// TryToDecay
02443           {
02444             G4QContent miQC=curString->GetQC();    // QContent of the Lightest Hadron
02445             G4int miPDG=miQC.GetSPDGCode();         // PDG of the Lightest Hadron
02446             if(miPDG == 10)                        // ==> Decay Hadron-Chipolino
02447             {
02448               G4QChipolino QCh(miQC);              // define theTotalResidual as aChipolino
02449               G4QPDGCode   h1QPDG=QCh.GetQPDG1();  // QPDG of the first hadron
02450               G4double     h1M   =h1QPDG.GetMass();// Mass of the first hadron
02451               G4QPDGCode   h2QPDG=QCh.GetQPDG2();  // QPDG of the second hadron 
02452               G4double     h2M   =h2QPDG.GetMass();// Mass of the second hadron
02453               G4double     ttM   =curString4M.m(); // Real Mass of the Chipolino
02454               if(h1M+h2M<ttM+eps)                  // Two particles decay of Chipolino
02455               {
02456                 G4LorentzVector h14M(0.,0.,0.,h1M);
02457                 G4LorentzVector h24M(0.,0.,0.,h2M);
02458                 if(std::fabs(ttM-h1M-h2M)<=eps)
02459                 {
02460                   G4double part1=h1M/(h1M+h2M);
02461                   h14M=part1*curString4M;
02462                   h24M=curString4M-h14M;
02463                 }
02464                 else
02465                 {
02466                   if(!G4QHadron(curString4M).DecayIn2(h14M,h24M))
02467                   {
02468                     G4cerr<<"***G4QIonIonCollision::Breeder: tM="<<ttM<<"->h1="<<h1QPDG
02469                           <<"("<<h1M<<")+h2="<<h1QPDG<<"("<<h2M<<")="<<h1M+h2M<<G4endl;
02470                     G4Exception("G4QIonIonCollision::Breeder:","72",FatalException,"QDec");
02471                   }
02472                 }
02473                 G4QHadron* h1H = new G4QHadron(h1QPDG.GetPDGCode(),h14M);
02474                 theResult->push_back(h1H);         // (delete equivalent)  
02475 #ifdef debug
02476                 G4LorentzVector f4M=h1H->Get4Momentum();
02477                 G4int           fPD=h1H->GetPDGCode();
02478                 G4int           fCg=h1H->GetCharge();
02479                 G4int           fBN=h1H->GetBaryonNumber();
02480                 G4cout<<"-EMC->>G4QIonIonCollision::Breed:String=HadrChiPro1's filled,f4M="
02481                       <<f4M<<", fPDG="<<fPD<<", fCg="<<fCg<<", fBN="<<fBN<<G4endl;
02482 #endif
02483                 G4QHadron* h2H = new G4QHadron(h2QPDG.GetPDGCode(),h24M);
02484                 theResult->push_back(h2H);         // (delete equivalent)  
02485 #ifdef debug
02486                 G4LorentzVector s4M=h2H->Get4Momentum();
02487                 G4int           sPD=h2H->GetPDGCode();
02488                 G4int           sCg=h2H->GetCharge();
02489                 G4int           sBN=h2H->GetBaryonNumber();
02490                 G4cout<<"-EMC->>G4QIonIonCollision::Breed:String=HadrChiPro2's filled,s4M="
02491                       <<s4M<<", sPDG="<<sPD<<", sCg="<<sCg<<", sBN="<<sBN<<G4endl;
02492 #endif
02493 #ifdef edebug
02494                 G4cout<<"-EMC-..Chi..G4QIonIonCollision::Breeder: DecayCHECK, Ch4M="
02495                       <<curString4M<<", d4M="<<curString4M-h14M-h24M<<G4endl;
02496 #endif
02497                 break;                               // Go out of the main StringDecay LOOP
02498               }
02499               else
02500               {
02501                 G4Quasmon* stringQuasmon = new G4Quasmon(miQC, curString4M);// String->Quas
02502                 theQuasmons.push_back(stringQuasmon);
02503                 break;                               // Go out of the main StringDecay LOOP
02504               }
02505             }
02506             else if(miPDG)                                 // Decay Hadron as a Resonans
02507             {
02508               if     (miPDG>0 &&   miPDG%10 < 3) miPDG+=2; // Convert to Resonans
02509               else if(miPDG<0 && (-miPDG)%10< 3) miPDG-=2; // Convert to antiResonans
02510               G4Quasmon Quasm;
02511               G4QHadron* sHad = new G4QHadron(miPDG,curString4M);
02512               G4QHadronVector* tmpQHadVec=Quasm.DecayQHadron(sHad); // It deletes sHad
02513               G4int tmpN=tmpQHadVec->size();
02514 #ifdef debug
02515               G4cout<<"G4QIonIonCollision::Breeder: Decay the Last, Res#H="<<tmpN<<G4endl;
02516 #endif
02517               if(tmpN>1)
02518               {
02519                 for(G4int aH=0; aH < tmpN; aH++)
02520                 {
02521                   theResult->push_back((*tmpQHadVec)[aH]);//TheDecayProdOfHadronDecIsFilled
02522 #ifdef debug
02523                   G4QHadron*   prodH =(*tmpQHadVec)[aH];
02524                   G4LorentzVector p4M=prodH->Get4Momentum();
02525                   G4int           PDG=prodH->GetPDGCode();
02526                   G4int           Chg=prodH->GetCharge();
02527                   G4int           BaN=prodH->GetBaryonNumber();
02528                   G4cout<<"-EMC->>G4QIonIonCollis::Breed:String=Hadr,H#"<<aH<<" filled,4M="
02529                         <<p4M<<", PDG="<<PDG<<", Chg="<<Chg<<", BaN="<<BaN<<G4endl;
02530 #endif
02531                 }
02532               }
02533               else
02534               {
02535                 G4Quasmon* stringQuasmon = new G4Quasmon(miQC, curString4M);// String->Quas
02536 #ifdef debug
02537                 G4cout<<"G4QIonIonCollision::Breeder:==> to Quasm="<<miQC<<curString4M
02538                       <<", pNuc="<<theProjNucleus<<theProjNucleus.Get4Momentum()<<", tNuc="
02539                       <<theTargNucleus<<theTargNucleus.Get4Momentum()<<", NString="
02540                       <<strings.size()<<", nR="<<theResult->size()<<", nQ="
02541                       <<theQuasmons.size()<<G4endl;
02542 #endif
02543                 theQuasmons.push_back(stringQuasmon);
02544                 delete sHad;
02545                 tmpQHadVec->clear();
02546                 delete tmpQHadVec;  // WhoCallsDecayQHadron is responsible for clear&delete
02547                 break;                               // Go out of the main StringDecay LOOP
02548               }
02549               tmpQHadVec->clear();
02550               delete tmpQHadVec;  // Who calls DecayQHadron is responsible for clear&delete
02551               break;                               // Go out of the main String Decay LOOP
02552             }
02553           } // End of the DecayOfTheLast
02554         } // End of IF(String-Hadron fusion)
02555       } // End of IF(NO_Hadrons) for String-String and String-Hadron fusion
02556       // The last hope is to CORREC the string, using other strings (ForwardInLOOP)
02557 #ifdef debug
02558       G4cout<<"G4QIonIonCollision::Breeder: theH="<<theHadrons<<"?=0, next="<<next<<G4endl;
02559 #endif
02560       if(!theHadrons && next < strings.size())       // ForwardInLOOP strings exist
02561       {
02562         // @@ string can be not convertable to one hadron (2,0.0,0,2,0) - To be improved
02563         G4QContent miQC=curString->GetQC(); // QContent of the Lightest Hadron
02564         G4int miPDG=miQC.GetSPDGCode();// PDG of the Lightest Hadron
02565 #ifdef debug
02566         G4cout<<">>>G4QIonIonCollision::Breeder: SQC="<<miQC<<", miSPDG="<<miPDG<<G4endl;
02567 #endif
02568         G4double miM=0.;               // Prototype of the Mass of the Cur LightestHadron
02569         if(miPDG!=10) miM=G4QPDGCode(miPDG).GetMass(); // Mass of theCurLightestOneHadron
02570         else
02571         {
02572           G4QChipolino QCh(miQC);      // define the TotalString as a Chipolino
02573           miM=QCh.GetQPDG1().GetMass()+QCh.GetQPDG2().GetMass();//MinMass of theStringChipo
02574         }
02575         G4double cM2=curString4M.m2(); // Actual squared mass of the Cur String
02576 #ifdef debug
02577         G4cout<<">>>G4QIonIonCollision::Breeder: minMass="<<miM<<", realM2="<<cM2<<G4endl;
02578 #endif
02579         G4double   cM=0.;
02580         if(cM2>0.)
02581         {
02582           cM=std::sqrt(cM2);
02583           if(std::fabs(cM-miM) < eps)    // Convert to hadron(2 hadrons) w/o calculations
02584           {
02585             if(miPDG!=10)
02586             {
02587               G4QHadron* sHad = new G4QHadron(miPDG,curString4M);
02588               theResult->push_back(sHad);// Fill the curString as a hadron
02589 #ifdef debug
02590               G4cout<<">>>G4QIonIonCollision::Breeder:S->H="<<miPDG<<curString4M<<G4endl;
02591 #endif
02592             }
02593             else
02594             {
02595               G4QChipolino QCh(miQC);               // define TotalResidual as a Chipolino
02596               G4QPDGCode   h1QPDG=QCh.GetQPDG1();   // QPDG of the first hadron
02597               G4double     h1M   =h1QPDG.GetMass(); // Mass of the first hadron
02598               G4QPDGCode   h2QPDG=QCh.GetQPDG2();   // QPDG of the second hadron 
02599               G4double     h2M   =h2QPDG.GetMass(); // Mass of the second hadron
02600               G4double     pt1   =h1M/(h1M+h2M);    // 4-mom part of the first hadron
02601               G4LorentzVector h14M=pt1*curString4M; // 4-mom of the first hadron
02602               G4LorentzVector h24M=curString4M-h14M;// 4-mom of the second hadron
02603               G4QHadron* h1H = new G4QHadron(h1QPDG.GetPDGCode(),h14M);
02604               theResult->push_back(h1H);            // (delete equivalent)  
02605 #ifdef debug
02606               G4LorentzVector f4M=h1H->Get4Momentum();
02607               G4int           fPD=h1H->GetPDGCode();
02608               G4int           fCg=h1H->GetCharge();
02609               G4int           fBN=h1H->GetBaryonNumber();
02610               G4cout<<"-EMC->>G4QIonIonCollision::Breed:Str=2HadrAR Prod-F is filled, f4M="
02611                     <<f4M<<", fPDG="<<fPD<<", fCg="<<fCg<<", fBN="<<fBN<<G4endl;
02612 #endif
02613               G4QHadron* h2H = new G4QHadron(h2QPDG.GetPDGCode(),h24M);
02614               theResult->push_back(h2H);         // (delete equivalent)  
02615 #ifdef debug
02616               G4LorentzVector s4M=h2H->Get4Momentum();
02617               G4int           sPD=h2H->GetPDGCode();
02618               G4int           sCg=h2H->GetCharge();
02619               G4int           sBN=h2H->GetBaryonNumber();
02620               G4cout<<"-EMC->>G4QIonIonCollision::Breed:Str=2HadrAR Prod-S is filled, s4M="
02621                     <<s4M<<", sPDG="<<sPD<<", sCg="<<sCg<<", sBN="<<sBN<<G4endl;
02622 #endif
02623             }
02624             continue;                    // Continue the LOOP over the curStrings
02625           }
02626           else                           // Try to recover (+/-) to min Mass
02627           {
02628             G4ThreeVector cP=curString4M.vect(); // Momentum of the curString
02629             G4double      cE=curString4M.e();    // Energy of the curString
02630             G4ThreeVector curV=cP/cE;    // curRelativeVelocity
02631             G4double miM2=miM*miM;
02632             G4int restr=0;            // To use beyon the LOOP for printing
02633             G4int fustr=0;               // Selected String-partner (0 = NotFound)
02634             G4double selX=0.;            // Selected value of x
02635             G4double maD=-DBL_MAX;       // Maximum Free Mass
02636             G4double Vmin=DBL_MAX;       // min Velocity Distance
02637             G4LorentzVector s4M(0.,0.,0.,0.); // Selected 4-mom of the hadron
02638 #ifdef debug
02639             G4cout<<"G4QIonIonCollision::Breeder: TryRecover, cV="<<curV<<G4endl;
02640 #endif
02641             nOfStr=strings.size();
02642             for(restr=next; restr < nOfStr; ++restr) if(restr != astring)
02643             {
02644               G4QString* pString=strings[restr];
02645               G4LorentzVector p4M=pString->Get4Momentum();
02646               G4ThreeVector pP=p4M.vect();  // Momentum of the partnerString
02647               G4double      pE=p4M.e();     // Energy of the partnerString
02648               G4double D2=cE*pE-cP.dot(pP); 
02649               G4double pM2=p4M.m2();
02650               G4double dM4=pM2*(miM2-cM2);
02651               G4double D=D2*D2+dM4;
02652               G4double x=-1.;                // Bad preexpectation 
02653               if(D >= 0. && pM2>.01) x=(std::sqrt(D)-D2)/pM2; // what we should get from p
02654 #ifdef debug
02655                   else G4cout<<"G4QIonIonCollis::Breed:D="<<D<<",D2="<<D2<<",dM="<<dM4<<G4endl;
02656               G4cout<<"G4QIonIonCollision::Breed: pM2="<<pM2<<",D2="<<D2<<",x="<<x<<G4endl;
02657 #endif
02658               if(x > 0. && x < 1.)          // We are getting x part of p4M
02659               {
02660                 G4QContent pQC=pString->GetQC(); // Quark Content of The Partner
02661                 G4int pPDG=pQC.GetSPDGCode();// PDG of The Lightest Hadron for the Partner
02662                 G4double pM=0.;             // Mass of the LightestHadron
02663                 if(pPDG==10)
02664                 {
02665                   G4QChipolino QCh(pQC);    // define the TotalString as a Chipolino
02666                   pM=QCh.GetQPDG1().GetMass()+QCh.GetQPDG2().GetMass();// Mass of Chipolino
02667                 }
02668                 else pM=G4QPDGCode(pPDG).GetMass();// Mass of theLightestHadron for Partner
02669                 G4double rM=std::sqrt(pM2); // Real mass of the string-partner
02670                 G4double delta=(1.-x)*rM-pM;// @@ Minimum CM disterbance measure
02671                 if(delta > 0. && delta > maD)
02672                 {
02673                   maD=delta;
02674 #ifdef debug
02675                   G4cout<<"G4QIonIonCollision::Breed: Subtr,S#"<<restr<<",d="<<maD<<G4endl;
02676 #endif
02677                   fustr=restr;
02678                   selX=x;
02679                   s4M=p4M;
02680                 }
02681               }
02682               else if(x <= 0.)               // We are adding to p4M, so use RelVelocity
02683               {
02684                 G4ThreeVector pV=pP/pE;      // curRelativeVelocity
02685                 G4double dV=(curV-pV).mag2();// SquaredDifferenceBetweenRelVelocities
02686                 if(dV < Vmin)
02687                 {
02688 #ifdef debug
02689                   G4cout<<"G4QIonIonCollision::Breed: FreeAdd,S#"<<restr<<",x="<<x<<G4endl;
02690 #endif
02691                   Vmin=dV;
02692                   fustr=restr;
02693                   selX=x;
02694                   s4M=p4M;
02695                 }
02696               }
02697 #ifdef debug
02698               G4cout<<"G4QIonIonCollision::Breed:EndOfLOOP r="<<restr<<"<"<<nOfStr<<G4endl;
02699 #endif
02700             } // End of the LOOP over string-partners for Correction
02701 #ifdef debug
02702               G4cout<<"G4QIonIonCollision::Breeder: AfterLOOP fustr="<<fustr<<G4endl;
02703 #endif
02704             if(fustr)
02705             {
02706 #ifdef edebug
02707               G4LorentzVector sum4M=s4M+curString4M;
02708               G4cout<<"G4QIonIonCollision::Breeder: Found Sum4M="<<sum4M<<G4endl;
02709 #endif
02710               G4QString* pString=strings[fustr];
02711               curString4M+=selX*s4M;
02712               if(std::abs(miPDG)%10 > 2)                  // Decay String-Hadron-Resonance
02713               {
02714                 G4Quasmon Quasm;
02715                 G4QHadron* sHad = new G4QHadron(miPDG,curString4M);
02716                 G4QHadronVector* tmpQHadVec=Quasm.DecayQHadron(sHad); // It deletes sHad
02717 #ifdef debug
02718                 G4cout<<"G4QIonIonCollision::Breed:DecStH,nH="<<tmpQHadVec->size()<<G4endl;
02719 #endif
02720                 for(unsigned aH=0; aH < tmpQHadVec->size(); aH++)
02721                 {
02722                   theResult->push_back((*tmpQHadVec)[aH]);//TheDecayProdOfHadron is filled
02723 #ifdef debug
02724                   G4QHadron*   prodH =(*tmpQHadVec)[aH];
02725                   G4LorentzVector p4M=prodH->Get4Momentum();
02726                   G4int           PDG=prodH->GetPDGCode();
02727                   G4int           Chg=prodH->GetCharge();
02728                   G4int           BaN=prodH->GetBaryonNumber();
02729                   G4cout<<"-EMC->>G4QIonIonCollision::Breed:St=Had,pH#"<<aH<<" filled, 4M="
02730                         <<p4M<<", PDG="<<PDG<<", Chg="<<Chg<<", BaN="<<BaN<<G4endl;
02731 #endif
02732                 }
02733                 tmpQHadVec->clear();
02734                 delete tmpQHadVec;  // Who calls DecayQHadron is responsibleRorClear&Delete
02735               }
02736               else if(miPDG == 10)                   // ==> Decay Hadron-Chipolino
02737               {
02738                 G4QChipolino QCh(miQC);              // define theTotalResid as aChipolino
02739                 G4QPDGCode   h1QPDG=QCh.GetQPDG1();  // QPDG of the first hadron
02740                 G4double     h1M   =h1QPDG.GetMass();// Mass of the first hadron
02741                 G4QPDGCode   h2QPDG=QCh.GetQPDG2();  // QPDG of the second hadron 
02742                 G4double     h2M   =h2QPDG.GetMass();// Mass of the second hadron
02743                 G4double     ttM   =curString4M.m(); // Real Mass of the Chipolino
02744                 if(h1M+h2M<ttM+eps)                  // Two particles decay of Chipolino
02745                 {
02746                   G4LorentzVector h14M(0.,0.,0.,h1M);
02747                   G4LorentzVector h24M(0.,0.,0.,h2M);
02748                   if(std::fabs(ttM-h1M-h2M)<=eps)
02749                   {
02750                     G4double part1=h1M/(h1M+h2M);
02751                     h14M=part1*curString4M;
02752                     h24M=curString4M-h14M;
02753                   }
02754                   else
02755                   {
02756                     if(!G4QHadron(curString4M).DecayIn2(h14M,h24M))
02757                     {
02758                       G4cerr<<"***G4QIonIonCollision::Breeder: tM="<<ttM<<"->h1="<<h1QPDG
02759                             <<"("<<h1M<<")+h2="<<h1QPDG<<"("<<h2M<<")="<<h1M+h2M<<G4endl;
02760                       G4Exception("G4QIonIonCollision::Breeder:","72",FatalException,"CD");
02761                     }
02762                   }
02763                   G4QHadron* h1H = new G4QHadron(h1QPDG.GetPDGCode(),h14M);
02764                   theResult->push_back(h1H);        // (delete equivalent)  
02765 #ifdef debug
02766                   G4LorentzVector f4M=h1H->Get4Momentum();
02767                   G4int           fPD=h1H->GetPDGCode();
02768                   G4int           fCg=h1H->GetCharge();
02769                   G4int           fBN=h1H->GetBaryonNumber();
02770                   G4cout<<"-EMC->>>G4QIonIonCollision::Breed:Str=Hadr Prod-F's filled,f4M="
02771                         <<f4M<<", fPDG="<<fPD<<", fCg="<<fCg<<", fBN="<<fBN<<G4endl;
02772 #endif
02773                   G4QHadron* h2H = new G4QHadron(h2QPDG.GetPDGCode(),h24M);
02774                   theResult->push_back(h2H);        // (delete equivalent)  
02775 #ifdef debug
02776                   G4LorentzVector s4M=h2H->Get4Momentum();
02777                   G4int           sPD=h2H->GetPDGCode();
02778                   G4int           sCg=h2H->GetCharge();
02779                   G4int           sBN=h2H->GetBaryonNumber();
02780                   G4cout<<"-EMC->>>G4QIonIonCollision::Breed:Str=Hadr Prod-S's filled,s4M="
02781                         <<s4M<<", sPDG="<<sPD<<", sCg="<<sCg<<", sBN="<<sBN<<G4endl;
02782 #endif
02783 #ifdef edebug
02784                   G4cout<<"-EMC-Chipo.G4QIonIonCollision::Breed:DecCHECK,c4M="<<curString4M
02785                         <<", ChQC="<<miQC<<", d4M="<<curString4M-h14M-h24M<<G4endl;
02786 #endif
02787                 }
02788                 else
02789                 {
02790                   G4cerr<<"***G4QInel::Breeder: tM="<<ttM<<miQC<<"->h1="<<h1QPDG<<"(" <<h1M
02791                         <<")+h2="<<h1QPDG<<"("<<h2M<<") = "<<h1M+h2M<<G4endl;
02792                   G4Exception("G4QIonIonCollision::Breeder:","72",FatalException,"ChiDec");
02793                 }
02794               }
02795               else
02796               {
02797                 G4QHadron* sHad = new G4QHadron(miPDG,curString4M);
02798                 theResult->push_back(sHad);         // The original string-hadron is filled
02799 #ifdef debug
02800                 G4cout<<"-EMC->>>G4QIonIonCollision::Breeder:Str=Hadr Filled, 4M="
02801                       <<curString4M<<", PDG="<<miPDG<<G4endl;
02802 #endif
02803               }
02804               G4double corF=1-selX;
02805               G4QParton* Left=pString->GetLeftParton();
02806               G4QParton* Right=pString->GetRightParton();
02807               Left->Set4Momentum(corF*Left->Get4Momentum());
02808               Right->Set4Momentum(corF*Right->Get4Momentum());
02809 #ifdef edebug
02810               G4cout<<"-EMC-...Cor...G4QIonIonCollision::Breeder:CorCHECK Sum="<<sum4M
02811                     <<" =? "<<curString4M+pString->Get4Momentum()<<", M="<<miM<<" ?= "
02812                     <<curString4M.m()<<G4endl;
02813 #endif
02814 #ifdef debug
02815               G4cout<<">>>G4QIonIonCollision::Breeder:*Corrected* String->Hadr="<<miPDG
02816                     <<curString4M<<" by String #"<<fustr<<G4endl;
02817 #endif
02818               continue;                            // Continue the LOOP over the curStrings
02819             } // End of Found combination for the String on string correction
02820           } // End of the Try-to-recover String+String correction algorithm
02821         } // End of IF(CM2>0.)
02822       } // End of IF(Can try to correct by String-String)
02823 #ifdef debug
02824       else G4cerr<<"***G4QIonIonCollision::Breed:**No SSCorrection**, next="<<next<<G4endl;
02825 #endif
02826       // ------------ At this point we can reduce the 3/-3 meson to 1/-1 meson ------------
02827       G4QParton* lpcS=curString->GetLeftParton();
02828       G4QParton* rpcS=curString->GetRightParton();
02829       G4int lPDGcS=lpcS->GetPDGCode();
02830       G4int rPDGcS=rpcS->GetPDGCode();
02831       if     (lPDGcS==3 && rPDGcS==-3)
02832       {
02833         lpcS->SetPDGCode( 1);
02834         rpcS->SetPDGCode(-1);
02835       }
02836       else if(lPDGcS==-3 && rPDGcS==3)
02837       {
02838         lpcS->SetPDGCode(-1);
02839         rpcS->SetPDGCode( 1);
02840       }
02841       // -------- Now the only hope is Correction, using the already prodused Hadrons -----
02842       G4int nofRH=theResult->size();            // #of resulting Hadrons
02843 #ifdef debug
02844       G4cout<<"G4QIonIonCollision::Breeder: theH="<<theHadrons<<", #OfH="<<nofRH<<G4endl;
02845 #endif
02846       if(!theHadrons && nofRH)                  // Hadrons are existing for SH Correction
02847       {
02848 #ifdef debug
02849         G4cerr<<"!G4QIonIonCollision::Breeder:CanTrySHCor, nH="<<theResult->size()<<G4endl;
02850 #endif
02851         // @@ string can be not convertable to one hadron (2,0.0,0,2,0) - To be improved
02852         G4QContent miQC=curString->GetQC();     // QContent of the Lightest Hadron
02853         G4int miPDG=miQC.GetSPDGCode();         // PDG of the Lightest Hadron
02854         G4double miM=0.;                        // Prototype ofMass of theCurLightestHadron
02855         if(miPDG==10)                           // Mass of the Cur Lightest ChipolinoHadron
02856         {
02857           G4QChipolino QCh(miQC);               // define the TotalString as a Chipolino
02858           miM=QCh.GetQPDG1().GetMass()+QCh.GetQPDG2().GetMass();//MinMass of theStringChipo
02859         }
02860         else miM=G4QPDGCode(miPDG).GetMass();   // Mass of the Cur Lightest Hadron
02861         G4double spM=0.;                        // Mass of the selected Hadron-Partner
02862         G4ThreeVector cP=curString4M.vect();    // Momentum of the curString
02863         G4double      cE=curString4M.e();       // Energy of the curString
02864         G4ThreeVector curV=cP/cE;               // curRelativeVelocity
02865         G4int reha=0;                           // Hadron # to use beyon the LOOP
02866         G4int fuha=0;                           // Selected Hadron-partner (0 = NotFound)
02867         G4double dMmin=DBL_MAX;                 // min Excess of the mass
02868         G4LorentzVector s4M(0.,0.,0.,0.);       // Selected 4-mom of the Hadron+String
02869         G4double sM=0.;                         // Selected Mass of the Hadron+String
02870         for (reha=next; reha < nofRH; reha++)   // LOOP over already collected hadrons
02871         {
02872           G4QHadron* pHadron=(*theResult)[reha];// Pointer to the current Partner-Hadron
02873           G4LorentzVector p4M=pHadron->Get4Momentum();
02874           G4double         pM=p4M.m();          // Mass of the Partner-Hadron
02875           G4LorentzVector t4M=p4M+curString4M;  // Total momentum of the compound
02876           G4double        tM2=t4M.m2();         // Squared total mass of the compound
02877           if(tM2 >= sqr(pM+miM+eps))            // Condition of possible correction
02878           {
02879             G4double tM=std::sqrt(tM2);         // Mass of the Hadron+String compound
02880             G4double dM=tM-pM-miM;              // Excess of the compound mass
02881             if(dM < dMmin)
02882             {
02883               dMmin=dM;
02884               fuha=reha;
02885               spM=pM;
02886               s4M=t4M;
02887               sM=tM;
02888             }
02889           }
02890         } // End of the LOOP over string-partners for Correction
02891         if(fuha)                                // The hadron-partner was found
02892         { 
02893           G4QHadron* pHadron=(*theResult)[fuha];// Necessary for update
02894           G4LorentzVector mi4M(0.,0.,0.,miM);   // Prototype of the new String=Hadron
02895           if(miM+spM<sM+eps)                    // Decay into CorrectedString+theSameHadron
02896           {
02897             G4LorentzVector sp4M(0.,0.,0.,spM);
02898             if(std::fabs(sM-miM-spM)<=eps)
02899             {
02900               G4double part1=miM/(miM+spM);
02901               mi4M=part1*s4M;
02902               sp4M=s4M-mi4M;
02903             }
02904             else
02905             {
02906               if(!G4QHadron(s4M).DecayIn2(mi4M,sp4M))
02907               {
02908                 G4cerr<<"***G4QIonIonCollision::Breeder: *SH*, tM="<<sM<<"->h1=("<<miPDG
02909                       <<")"<<miM<<" + h2="<<spM<<" = "<<miM+spM<<G4endl;
02910                 G4Exception("G4QIonIonCollision::Breeder:","72",FatalException,"SHChiDec");
02911               }
02912             }
02913             pHadron->Set4Momentum(sp4M);
02914 #ifdef debug
02915             G4cout<<"-EMC->...>G4QIonIonCollision::Breed: H# "<<fuha<<" is updated, new4M="
02916                   <<sp4M<<G4endl;
02917 #endif
02918           }
02919           else
02920           {
02921             G4cerr<<"***G4QInel::Breeder: HS Failed, tM="<<sM<<"->h1M=("<<miPDG<<")"<<miM
02922                   <<"+h2M="<<spM<<" = "<<miM+spM<<G4endl;
02923             G4Exception("G4QIonIonCollision::Breeder:","72",FatalException,"HSChipoliDec");
02924           }
02925           if(std::abs(miPDG)%10 > 2)                  // Decay Hadron-Resonans
02926           {
02927             G4Quasmon Quasm;
02928             G4QHadron* sHad = new G4QHadron(miPDG,mi4M);
02929             G4QHadronVector* tmpQHadVec=Quasm.DecayQHadron(sHad); // It deletes sHad
02930 #ifdef debug
02931             G4cout<<"G4QInelast::Breeder: *HS* DecStrHad, nH="<<tmpQHadVec->size()<<G4endl;
02932 #endif
02933             for(unsigned aH=0; aH < tmpQHadVec->size(); aH++)
02934             {
02935               theResult->push_back((*tmpQHadVec)[aH]);// TheDecayProductOfTheHadronIsFilled
02936 #ifdef debug
02937               G4QHadron*   prodH =(*tmpQHadVec)[aH];
02938               G4LorentzVector p4M=prodH->Get4Momentum();
02939               G4int           PDG=prodH->GetPDGCode();
02940               G4int           Chg=prodH->GetCharge();
02941               G4int           BaN=prodH->GetBaryonNumber();
02942               G4cout<<"-EMC->>>G4QIonIonCollision::Breed:Str+Hadr PrH#"<<aH<<" filled, 4M="
02943                     <<p4M<<", PDG="<<PDG<<", Chg="<<Chg<<", BaN="<<BaN<<G4endl;
02944 #endif
02945             }
02946             tmpQHadVec->clear();
02947             delete tmpQHadVec;  // Who calls DecayQHadron is responsible for clear & delete
02948           }
02949           else if(miPDG == 10)                   // ==> Decay Hadron-Chipolino
02950           {
02951             G4QChipolino QCh(miQC);              // define the TotalResidual as a Chipolino
02952             G4QPDGCode   h1QPDG=QCh.GetQPDG1();  // QPDG of the first hadron
02953             G4double     h1M   =h1QPDG.GetMass();// Mass of the first hadron
02954             G4QPDGCode   h2QPDG=QCh.GetQPDG2();  // QPDG of the second hadron 
02955             G4double     h2M   =h2QPDG.GetMass();// Mass of the second hadron
02956             G4double     ttM   =curString4M.m(); // Real Mass of the Chipolino
02957             if(h1M+h2M<miM+eps)                  // Two particles decay of Chipolino
02958             {
02959               G4LorentzVector h14M(0.,0.,0.,h1M);
02960               G4LorentzVector h24M(0.,0.,0.,h2M);
02961               if(std::fabs(ttM-h1M-h2M)<=eps)
02962               {
02963                 G4double part1=h1M/(h1M+h2M);
02964                 h14M=part1*mi4M;
02965                 h24M=mi4M-h14M;
02966               }
02967               else
02968               {
02969                 if(!G4QHadron(mi4M).DecayIn2(h14M,h24M))
02970                 {
02971                   G4cerr<<"***G4QIonIonCollision::Breeder: HS tM="<<ttM<<"->h1="<<h1QPDG
02972                         <<"("<<h1M<<")+h2="<<h1QPDG<<"("<<h2M<<")="<<h1M+h2M<<G4endl;
02973                   G4Exception("G4QIonIonCollision::Breeder:","72",FatalException,"ChiDec");
02974                 }
02975               }
02976               G4QHadron* h1H = new G4QHadron(h1QPDG.GetPDGCode(),h14M);
02977               theResult->push_back(h1H);         // (delete equivalent)  
02978 #ifdef debug
02979               G4LorentzVector f4M=h1H->Get4Momentum();
02980               G4int           fPD=h1H->GetPDGCode();
02981               G4int           fCg=h1H->GetCharge();
02982               G4int           fBN=h1H->GetBaryonNumber();
02983               G4cout<<"-EMC->>G4QIonIonCollision::Breed: CorStrHadr Prod-1 is filled, f4M="
02984                     <<f4M<<", fPDG="<<fPD<<", fCg="<<fCg<<", fBN="<<fBN<<G4endl;
02985 #endif
02986               G4QHadron* h2H = new G4QHadron(h2QPDG.GetPDGCode(),h24M);
02987               theResult->push_back(h2H);         // (delete equivalent)  
02988 #ifdef debug
02989               G4LorentzVector n4M=h2H->Get4Momentum();
02990               G4int           nPD=h2H->GetPDGCode();
02991               G4int           nCg=h2H->GetCharge();
02992               G4int           nBN=h2H->GetBaryonNumber();
02993               G4cout<<"-EMC->>>G4QIonIonCollision::Breed:CorStrHadr Prod-2 is filled, n4M="
02994                     <<n4M<<", nPDG="<<nPD<<", nCg="<<nCg<<", nBN="<<nBN<<G4endl;
02995 #endif
02996 #ifdef edebug
02997               G4cout<<"-EMC-...HS-Chipo...G4QIonIonCollision::Breeder:DecCHECK, Ch4M="
02998                     <<mi4M<<", ChQC="<<miQC<<", d4M="<<mi4M-h14M-h24M<<G4endl;
02999 #endif
03000             }
03001           }
03002           else
03003           {
03004             G4QHadron* sHad = new G4QHadron(miPDG, mi4M);
03005             theResult->push_back(sHad);          // The original string=hadron is filled
03006 #ifdef debug
03007             G4cout<<">>..>>G4QIonIonCollision::Breeder: CorStr=Hadr is Filled, 4M="
03008                   <<curString4M<<", StPDG="<<miPDG<<G4endl;
03009 #endif
03010           }
03011 #ifdef edebug
03012           G4cout<<"-EMC-...Cor...G4QIonIonCollision::Breeder:StHadCor CHECK Sum="<<s4M
03013                 <<" =? "<<mi4M+pHadron->Get4Momentum()<<G4endl;
03014 #endif
03015 #ifdef debug
03016           G4cout<<">>>G4QIonIonCollision::Breeder:*Corrected* String+Hadr="<<miPDG
03017                 <<mi4M<<" by Hadron #"<<reha<<G4endl;
03018 #endif
03019           continue;                    // Continue the LOOP over the curStrings
03020         }
03021         else
03022         {
03023 #ifdef debug
03024           G4cout<<"-EMC->>>G4QIonIonCollision::Breeder: Str+Hadr Failed, 4M="<<curString4M
03025                 <<", PDG="<<miPDG<<G4endl;
03026 #endif
03027         }
03028         // @@@ convert string to Quasmon with curString4M
03029         G4QContent curStringQC=curString->GetQC();
03030         G4Quasmon* stringQuasmon = new G4Quasmon(curStringQC, curString4M);
03031         theQuasmons.push_back(stringQuasmon);
03032         continue;                      // Continue the LOOP over the curStrings
03033       } // End of IF(Can try the String-Hadron correction
03034     } // End of IF(NO_Hadrons) = Problem solution namespace
03035     G4Quasmon tmpQ;                                 // @@ an issue of Q to decay resonances
03036     G4int nHfin=0;
03037     if(theHadrons) nHfin=theHadrons->size();
03038     else // !! Sum Up all strings and convert them in a Quasmon (Exception for development)
03039     {
03040       G4LorentzVector ss4M(0.,0.,0.,0.);
03041       G4QContent      ssQC(0,0,0,0,0,0);
03042       G4int tnSt=strings.size();
03043       for(G4int i=astring; i < tnSt; ++i)
03044       {
03045         G4LorentzVector pS4M=strings[i]->Get4Momentum(); // String 4-momentum
03046         ss4M+=pS4M;
03047         G4QContent pSQC=strings[i]->GetQC();             // String Quark Content
03048         ssQC+=pSQC;
03049 #ifdef debug
03050         G4cout<<"=--=>G4QIonIonCollision::Breed:S#"<<i<<",4M="<<pS4M<<",QC="<<pSQC<<G4endl;
03051 #endif
03052       }
03053 #ifdef debug
03054       G4cout<<"==>G4QIonIonCollision::Breed:AllStrings are summed up in a Quasmon"<<G4endl;
03055 #endif
03056       G4Quasmon* stringQuasmon = new G4Quasmon(ssQC, ss4M);
03057       theQuasmons.push_back(stringQuasmon);
03058       break;                                   // break the LOOP over Strings
03059     }
03060 #ifdef debug
03061     G4cout<<"G4QIonIonCollision::Breeder: Trying to decay hadrons #ofHRes="<<nHfin<<G4endl;
03062 #endif
03063     for(G4int aTrack=0; aTrack<nHfin; aTrack++)
03064     {
03065       G4QHadron* curHadron=(*theHadrons)[aTrack];
03066       G4int hPDG=curHadron->GetPDGCode();
03067 #ifdef edebug
03068       G4LorentzVector curH4M=curHadron->Get4Momentum();
03069       G4int           curHCh=curHadron->GetCharge();
03070       G4int           curHBN=curHadron->GetBaryonNumber();
03071 #endif
03072 #ifdef debug
03073       G4cout<<">>...>>G4QIonIonCollision::Breeder:S#"<<astring<<",H#"<<aTrack<<",PDG="
03074             <<hPDG<<",4M="<<curHadron->Get4Momentum()<<G4endl;
03075 #endif
03076       if(std::abs(hPDG)%10 > 2)
03077       {
03078         G4QHadronVector* tmpQHadVec=tmpQ.DecayQHadron(curHadron); // It deletes curHadron
03079 #ifdef debug
03080         G4cout<<"G4QIonIonCollision::Breed:-DECAY'S DONE-,nH="<<tmpQHadVec->size()<<G4endl;
03081 #endif
03082         //G4int tmpS=tmpQHadVec->size(); // "The elegant method" (tested) is commented
03083         //theResult->resize(tmpS+theResult->size()); // Resize theQHadrons length
03084         //copy(tmpQHadVec->begin(), tmpQHadVec->end(), theResult->end()-tmpS);
03085         for(unsigned aH=0; aH < tmpQHadVec->size(); aH++)
03086         {
03087           theResult->push_back((*tmpQHadVec)[aH]);// TheDecayProduct of TheHadron is filled
03088 #ifdef edebug
03089           G4QHadron*   prodH =(*tmpQHadVec)[aH];
03090           G4LorentzVector p4M=prodH->Get4Momentum();
03091           G4int           PDG=prodH->GetPDGCode();
03092           G4int           Chg=prodH->GetCharge();
03093           G4int           BaN=prodH->GetBaryonNumber();
03094           curH4M-=p4M;
03095           curString4M-=p4M;
03096           curStrChg-=Chg;
03097           curStrBaN-=BaN;
03098           curHCh-=Chg;
03099           curHBN-=BaN;
03100           G4cout<<"-EMC->.>>G4QIonIonCollision::Breed:Str*Filled, 4M="<<p4M<<", PDG="<<PDG
03101                 <<", Chg="<<Chg<<", BaN="<<BaN<<G4endl;
03102 #endif
03103         }
03104 #ifdef edebug
03105         G4cout<<"-EMC-.G4QIn::Br:Dec,r4M="<<curH4M<<",rC="<<curHCh<<",rB="<<curHBN<<G4endl;
03106 #endif
03107         tmpQHadVec->clear();
03108         delete tmpQHadVec;  // Who calls DecayQHadron is responsible for clear & delete
03109       }
03110       else                                      // Chipolino is not checked here
03111       {
03112         theResult->push_back(curHadron);        // The original hadron is filled
03113 #ifdef edebug
03114         curString4M-=curH4M;
03115         G4int curCh=curHadron->GetCharge();
03116         G4int curBN=curHadron->GetBaryonNumber();
03117         curStrChg-=curCh;
03118         curStrBaN-=curBN;
03119         G4cout<<"-EMC->>..>>G4QIonIonCollision::Breeder: curH filled 4M="<<curH4M<<",PDG="
03120               <<curHadron->GetPDGCode()<<", Chg="<<curCh<<", BaN="<<curBN<<G4endl;
03121 #endif
03122       }
03123     }
03124     // clean up (the issues are filled to theResult)
03125     if(theHadrons) delete theHadrons;
03126 #ifdef edebug
03127     G4cout<<"-EMC-.......G4QIonIonCollision::Breeder: StringDecay CHECK, r4M="<<curString4M
03128           <<", rChg="<<curStrChg<<", rBaN="<<curStrBaN<<G4endl;
03129 #endif
03130   } // End of the main LOOP over decaying strings
03131   G4LorentzVector rp4M=theProjNucleus.Get4Momentum(); // projNucleus 4-momentum in LS
03132   G4int rpPDG=theProjNucleus.GetPDG();
03133   G4QHadron* resPNuc = new G4QHadron(rpPDG,rp4M);
03134   theResult->push_back(resPNuc);                      // Fill the residual nucleus
03135   G4LorentzVector rt4M=theTargNucleus.Get4Momentum(); // Nucleus 4-momentum in LS
03136   G4int rtPDG=theTargNucleus.GetPDG();
03137   G4QHadron* resTNuc = new G4QHadron(rtPDG,rt4M);
03138   theResult->push_back(resTNuc);                      // Fill the residual nucleus
03139 #ifdef edebug
03140   G4LorentzVector s4M(0.,0.,0.,0.);                   // Sum of the Result in LS
03141   G4int rCh=totChg;
03142   G4int rBN=totBaN;
03143   G4int nHadr=theResult->size();
03144   G4int nQuasm=theQuasmons.size();
03145   G4cout<<"-EMCLS-G4QInel::Breeder: #ofHadr="<<nHadr<<", #OfQ="<<nQuasm<<", rPA="<<rp4M.m()
03146         <<"="<<G4QNucleus(rpPDG).GetGSMass()<<", rTA="<<rt4M.m()<<"="
03147         <<G4QNucleus(rtPDG).GetGSMass()<<G4endl;
03148   for(G4int i=0; i<nHadr; i++)
03149   {
03150     G4LorentzVector hI4M=(*theResult)[i]->Get4Momentum();
03151     s4M+=hI4M;
03152     G4int hChg=(*theResult)[i]->GetCharge();
03153     rCh-=hChg;
03154     G4int hBaN=(*theResult)[i]->GetBaryonNumber();
03155     rBN-=hBaN;
03156     G4cout<<"-EMCLS-G4QIonIonCollision::Breeder: Hadron#"<<i<<", 4M="<<hI4M<<", PDG="
03157           <<(*theResult)[i]->GetPDGCode()<<", C="<<hChg<<", B="<<hBaN<<G4endl;
03158   }
03159   for(G4int i=0; i<nQuasm; i++)
03160   {
03161     G4LorentzVector hI4M=theQuasmons[i]->Get4Momentum();
03162     s4M+=hI4M;
03163     G4int hChg=theQuasmons[i]->GetCharge();
03164     rCh-=hChg;
03165     G4int hBaN=theQuasmons[i]->GetBaryonNumber();
03166     rBN-=hBaN;
03167     G4cout<<"-EMCLS-G4QIonIonCollision::Breeder: Quasmon#"<<i<<", 4M="<<hI4M<<", C="<<hChg
03168           <<", B="<<hBaN<<G4endl;
03169   }
03170   G4cout<<"-EMCLS-G4QInel::Breed: LS r4M="<<s4M-totLS4M<<",rC="<<rCh<<",rB="<<rBN<<G4endl;
03171 #endif
03172   // Now we need to coolect particles for creation of a Quasmon @@ improve !!
03173   G4int nRes=theResult->size();
03174   if(nRes>2)
03175   {
03176     G4QHadronVector::iterator ih;
03177     G4QHadronVector::iterator nih;
03178     G4QHadronVector::iterator mih;
03179     G4QHadronVector::iterator lst=theResult->end();
03180     lst--;
03181     G4double minMesEn=DBL_MAX;
03182     G4double minBarEn=DBL_MAX;
03183     G4bool nfound=false;
03184     G4bool mfound=false;
03185     for(ih = theResult->begin(); ih < theResult->end(); ++ih) if(ih != lst)
03186     {
03187       G4LorentzVector h4M=(*ih)->Get4Momentum();
03188       G4int          hPDG=(*ih)->GetPDGCode();
03189 #ifdef debug
03190       G4cout<<"%->G4QIonIonCollision::Breeder: TRY hPDG="<<hPDG<<", h4M="<<h4M<<G4endl;
03191 #endif
03192       if(hPDG>1111 && hPDG<3333)
03193       {
03194         G4double bE=h4M.e()-(*ih)->GetMass();
03195         if(bE < minBarEn)
03196         {
03197           minBarEn=bE;
03198           nih=ih;
03199           nfound=true;
03200         }
03201       }
03202       else if(hPDG>-1111)
03203       {
03204         G4double mE=h4M.e();
03205         if(mE < minMesEn)
03206         {
03207           minMesEn=mE;
03208           mih=ih;
03209           mfound=true;
03210         }
03211       }
03212     }
03213     if(nfound && mfound && minMesEn+minBarEn < maxEn)
03214     {
03215       G4QHadron* resNuc = (*theResult)[nRes-1];              // ResidualNucleus is theLastH
03216       theResult->pop_back();                                 // Fill the QHad-nucleus later
03217       G4LorentzVector q4M=(*nih)->Get4Momentum()+(*mih)->Get4Momentum();
03218       G4QContent qQC=(*nih)->GetQC()+(*mih)->GetQC();
03219       maxEn -= minBarEn+minMesEn;                            // Reduce the energy limit
03220 #ifdef debug
03221       G4cout<<"%->G4QIonIonCollision::Breeder:Exclude,4M="<<(*nih)->Get4Momentum()
03222             <<", & 4m="<<(*mih)->Get4Momentum()<<G4endl;
03223 #endif
03224       delete (*nih);
03225       delete (*mih);
03226       if(nih > mih)                                          // Exclude nucleon first
03227       {
03228         theResult->erase(nih);                               // erase Baryon from theResult
03229         theResult->erase(mih);                               // erase Meson from theResult
03230       }
03231       else                                                   // Exclude meson first
03232       {
03233         theResult->erase(mih);                               // erase Baryon from theResult
03234         theResult->erase(nih);                               // erase Meson from theResult
03235       }
03236 #ifdef debug
03237       G4cout<<"%->G4QI::Breed: BeforeLOOP, dE="<<maxEn<<", nR="<<theResult->size()<<G4endl;
03238 #endif
03239       if(maxEn > 0.)                                         // More hadrons to absorbe
03240       {
03241         for(ih = theResult->begin(); ih < theResult->end(); ih++)
03242         {
03243           G4LorentzVector h4M=(*ih)->Get4Momentum();
03244           G4int          hPDG=(*ih)->GetPDGCode();
03245           G4double dE=0.;
03246           if     (hPDG> 1111 && hPDG<3333) dE=h4M.e()-(*ih)->GetMass(); // Baryons
03247           else if(hPDG>-1111) dE=h4M.e();                    // Mesons
03248           else continue;                                     // Don't consider anti-baryons
03249           if(dE < maxEn)
03250           {
03251             maxEn-=dE;
03252             q4M+=h4M;
03253             qQC+=(*ih)->GetQC();
03254 #ifdef debug
03255             G4cout<<"%->G4QIonIonCollision::Breed:Exclude,4M="<<h4M<<",dE="<<maxEn<<G4endl;
03256 #endif
03257             delete (*ih);
03258             theResult->erase(ih);                            // erase Hadron from theResult
03259             --ih;
03260           }
03261         }
03262       }
03263       G4Quasmon* softQuasmon = new G4Quasmon(qQC, q4M);      // SoftPart -> Quasmon
03264 #ifdef debug
03265       G4cout<<"%->G4QIonIonCollision::Breed:QuasmonIsFilled,4M="<<q4M<<",QC="<<qQC<<G4endl;
03266 #endif
03267       theQuasmons.push_back(softQuasmon);
03268       theResult->push_back(resNuc);
03269     }
03270 #ifdef edebug
03271     G4LorentzVector f4M(0.,0.,0.,0.);                        // Sum of the Result in LS
03272     G4int fCh=totChg;
03273     G4int fBN=totBaN;
03274     G4int nHd=theResult->size();
03275     G4int nQm=theQuasmons.size();
03276     G4cout<<"-EMCLS-G4QIonIonCollision::Breeder:#ofHadr="<<nHd<<", #OfQ="<<nQm
03277           <<",rPA="<<rt4M.m()<<"="<<G4QNucleus(rtPDG).GetGSMass()
03278           <<",rTA="<<rt4M.m()<<"="<<G4QNucleus(rtPDG).GetGSMass()<<G4endl;
03279     for(G4int i=0; i<nHd; i++)
03280     {
03281       G4LorentzVector hI4M=(*theResult)[i]->Get4Momentum();
03282       f4M+=hI4M;
03283       G4int hChg=(*theResult)[i]->GetCharge();
03284       fCh-=hChg;
03285       G4int hBaN=(*theResult)[i]->GetBaryonNumber();
03286       fBN-=hBaN;
03287       G4cout<<"-EMCLS-G4QIonIonCollision::Breeder: Hadron#"<<i<<", 4M="<<hI4M<<", PDG="
03288             <<(*theResult)[i]->GetPDGCode()<<", C="<<hChg<<", B="<<hBaN<<G4endl;
03289     }
03290     for(G4int i=0; i<nQm; i++)
03291     {
03292       G4LorentzVector hI4M=theQuasmons[i]->Get4Momentum();
03293       f4M+=hI4M;
03294       G4int hChg=theQuasmons[i]->GetCharge();
03295       fCh-=hChg;
03296       G4int hBaN=theQuasmons[i]->GetBaryonNumber();
03297       fBN-=hBaN;
03298       G4cout<<"-EMCLS-G4QIonIonCollision::Breeder: Quasmon#"<<i<<", 4M="<<hI4M<<", C="
03299             <<hChg<<", B="<<hBaN<<G4endl;
03300     }
03301     G4cout<<"-EMCLS-G4QInel::Breed:, r4M="<<f4M-totLS4M<<", rC="<<fCh<<",rB="<<fBN<<G4endl;
03302 #endif
03303   } // End of the soft Quasmon Creation
03304   return;
03305 } // End of Breeder

G4double G4QIonIonCollision::ChooseX ( G4double  Xmin,
G4double  Xmax 
) const [protected]

Definition at line 3551 of file G4QIonIonCollision.cc.

References FatalException, G4cerr, G4cout, G4endl, G4Exception(), G4UniformRand, and sqr().

Referenced by ExciteDiffParticipants(), and ExciteSingDiffParticipants().

03552 {
03553   // choose an x between Xmin and Xmax with P(x) ~ 1/x @@ M.K. -> 1/sqrt(x)
03554   //G4double range=Xmax-Xmin;
03555   if(Xmax == Xmin) return Xmin;
03556   if( Xmin < 0. || Xmax < Xmin) 
03557   {
03558     G4cerr<<"***G4QIonIonCollision::ChooseX: Xmin="<<Xmin<<", Xmax="<<Xmax<< G4endl;
03559     G4Exception("G4QIonIonCollision::ChooseX:","72",FatalException,"Bad X or X-Range");
03560   }
03561   //G4double x;
03562   //do {x=Xmin+G4UniformRand()*range;} while ( Xmin/x < G4UniformRand() );
03563   G4double sxi=std::sqrt(Xmin);
03564   G4double x=sqr(sxi+G4UniformRand()*(std::sqrt(Xmax)-sxi));
03565 #ifdef debug
03566   G4cout<<"G4QIonIonCollision::ChooseX: DiffractiveX="<<x<<G4endl;
03567 #endif
03568   return x;
03569 } // End of ChooseX

G4bool G4QIonIonCollision::ExciteDiffParticipants ( G4QHadron aPartner,
G4QHadron bPartner 
) const [protected]

Definition at line 3308 of file G4QIonIonCollision.cc.

References ChooseX(), G4cout, G4endl, GaussianPt(), G4QHadron::Get4Momentum(), G4QHadron::GetMass(), G4QHadron::Set4Momentum(), and sqr().

Referenced by G4QIonIonCollision().

03310 {
03311   G4LorentzVector Pprojectile=projectile->Get4Momentum();
03312   G4double Mprojectile=projectile->GetMass();
03313   G4double Mprojectile2=Mprojectile*Mprojectile;
03314   G4LorentzVector Ptarget=target->Get4Momentum();
03315   G4double Mtarget=target->GetMass();
03316   G4double Mtarget2=Mtarget*Mtarget;
03317 #ifdef debug
03318   G4cout<<"G4QInel::ExciteDiffPartici: Ep="<<Pprojectile.e()<<", Et="<<Ptarget.e()<<G4endl;
03319 #endif
03320   // Transform momenta to cms and then rotate parallel to z axis;
03321   G4LorentzVector Psum=Pprojectile+Ptarget;
03322   G4LorentzRotation toCms(-Psum.boostVector()); // Boost Rotation to CMS
03323   G4LorentzVector Ptmp=toCms*Pprojectile;
03324   if(Ptmp.pz()<=0.) // "String" moving backwards in CMS, abort collision !! ? M.K.
03325   {
03326 #ifdef debug
03327     G4cout<<"G4QIonIonCollision::ExciteDiffParticipants:*1* abort Collision!! *1*"<<G4endl;
03328 #endif
03329     return false; 
03330   }         
03331   toCms.rotateZ(-Ptmp.phi());
03332   toCms.rotateY(-Ptmp.theta());
03333 #ifdef debug
03334   G4cout<<"G4QIonIonCollision::ExciteDiffParticipantts:BeforBoost Pproj="<<Pprojectile
03335         <<",Ptarg="<<Ptarget<<G4endl;
03336 #endif
03337   G4LorentzRotation toLab(toCms.inverse()); // Boost Rotation to LabSys (LS)
03338   Pprojectile.transform(toCms);
03339   Ptarget.transform(toCms);
03340 #ifdef debug
03341   G4cout<< "G4QInelast::ExciteDiffParticipantts: AfterBoost Pproj="<<Pprojectile<<",Ptarg="
03342         <<Ptarget<<", cms4M="<<Pprojectile+Ptarget<<G4endl;
03343   G4cout<<"G4QIonIonCollis::ExciteDiffParticipants:ProjX+="<<Pprojectile.plus()<<",ProjX-="
03344         <<Pprojectile.minus()<<", tX+="<< Ptarget.plus()<<",tX-="<<Ptarget.minus()<<G4endl;
03345 #endif
03346   G4LorentzVector Qmomentum(0.,0.,0.,0.);
03347   G4int whilecount=0;
03348 #ifdef debug
03349   G4cout<<"G4QIonIonCollision::ExciteDiffParticipants: Before DO"<<G4endl;
03350 #endif
03351   do
03352   {
03353     //  Generate pt  
03354     G4double maxPtSquare=sqr(Ptarget.pz());
03355 #ifdef debug
03356     G4cout<<"G4QIonIonCollision::ExciteDiffParticipants: maxPtSq="<<maxPtSquare<<G4endl;
03357     if(whilecount++>=500 && whilecount%100==0) // @@ M.K. Hardwired limits 
03358     G4cout<<"G4QIonIonCollision::ExciteDiffParticipantts: can loop, loopCount="<<whilecount
03359           <<", maxPtSquare="<<maxPtSquare<<G4endl;
03360 #endif
03361     if(whilecount>1000)                        // @@ M.K. Hardwired limits 
03362     {
03363 #ifdef debug
03364       G4cout<<"G4QIonIonCollision::ExciteDiffParticipants: *2* abort Loop!! *2*"<<G4endl;
03365 #endif
03366       return false;    //  Ignore this interaction 
03367     }
03368     Qmomentum=G4LorentzVector(GaussianPt(widthOfPtSquare,maxPtSquare),0);
03369 #ifdef debug
03370     G4cout<<"G4QIonIonCollision::ExciteDiffParticipants: generated Pt="<<Qmomentum
03371           <<", ProjPt="<<Pprojectile+Qmomentum<<", TargPt="<<Ptarget-Qmomentum<<G4endl;
03372 #endif
03373     //  Momentum transfer
03374     G4double Xmin=0.;
03375     G4double Xmax=1.;
03376     G4double Xplus =ChooseX(Xmin,Xmax);
03377     G4double Xminus=ChooseX(Xmin,Xmax);
03378 #ifdef debug
03379     G4cout<<"G4QIonIonCollision::ExciteDiffParticipant:X+="<<Xplus<<",X-="<<Xminus<<G4endl;
03380 #endif
03381     G4double pt2=Qmomentum.vect().mag2();
03382     G4double Qplus =-pt2/Xminus/Ptarget.minus();
03383     G4double Qminus= pt2/Xplus /Pprojectile.plus();
03384     Qmomentum.setPz((Qplus-Qminus)/2);
03385     Qmomentum.setE( (Qplus+Qminus)/2);
03386 #ifdef debug
03387     G4cout<<"G4QInelast::ExciteDiffParticip: Qplus="<<Qplus<<", Qminus="<<Qminus<<", pt2="
03388           <<pt2<<", Qmomentum="<<Qmomentum<<", ProjM="<<(Pprojectile+Qmomentum).mag()
03389           <<", TargM="<<(Ptarget-Qmomentum).mag()<<G4endl;
03390 #endif
03391   } while((Pprojectile+Qmomentum).mag2()<=Mprojectile2 ||
03392           (Ptarget-Qmomentum).mag2()<=Mtarget2);
03393   Pprojectile += Qmomentum;
03394   Ptarget     -= Qmomentum;
03395 #ifdef debug
03396   G4cout<<"G4QInelast::ExciteDiffParticipan: Proj(Q)="<<Pprojectile<<", Targ(Q)="<<Ptarget
03397         <<", Proj(back)="<<toLab*Pprojectile<<", Targ(bac)="<< toLab*Ptarget << G4endl;
03398 #endif
03399   // Transform back and update SplitableHadron Participant.
03400   Pprojectile.transform(toLab);
03401   Ptarget.transform(toLab);
03402 #ifdef debug
03403   G4cout<< "G4QIonIonCollision::ExciteDiffParticipants:TargetMass="<<Ptarget.mag()<<G4endl;
03404 #endif
03405   target->Set4Momentum(Ptarget);  
03406 #ifdef debug
03407   G4cout<<"G4QIonIonCollision::ExciteDiffParticipant:ProjMass="<<Pprojectile.mag()<<G4endl;
03408 #endif
03409   projectile->Set4Momentum(Pprojectile);
03410   return true;
03411 } // End of ExciteDiffParticipants

G4bool G4QIonIonCollision::ExciteSingDiffParticipants ( G4QHadron aPartner,
G4QHadron bPartner 
) const [protected]

Definition at line 3415 of file G4QIonIonCollision.cc.

References ChooseX(), G4cout, G4endl, G4UniformRand, GaussianPt(), G4QHadron::Get4Momentum(), G4QHadron::GetMass(), G4QHadron::GetMass2(), G4QHadron::Set4Momentum(), and sqr().

Referenced by G4QIonIonCollision().

03417 {
03418   G4LorentzVector Pprojectile=projectile->Get4Momentum();
03419   G4double Mprojectile=projectile->GetMass();
03420   G4double Mprojectile2=Mprojectile*Mprojectile;
03421   G4LorentzVector Ptarget=target->Get4Momentum();
03422   G4double Mtarget=target->GetMass();
03423   G4double Mtarget2=Mtarget*Mtarget;
03424 #ifdef debug
03425   G4cout<<"G4QInel::ExSingDiffPartici: Ep="<<Pprojectile.e()<<", Et="<<Ptarget.e()<<G4endl;
03426 #endif
03427   G4bool KeepProjectile= G4UniformRand() > 0.5;
03428   // Reset minMass of the non diffractive particle to its value, (minus for rounding...)
03429   if(KeepProjectile ) 
03430   {
03431 #ifdef debug
03432     G4cout<<"-1/2-G4QIonIonCollision::ExSingDiffParticipants: Projectile is fixed"<<G4endl;
03433 #endif
03434     Mprojectile2 = projectile->GetMass2()*(1.-perCent); // Isn't it too big reduction? M.K.
03435   }
03436   else
03437   {
03438 #ifdef debug
03439     G4cout<<"---1/2---G4QIonIonCollision::ExSingDiffParticipants: Target is fixed"<<G4endl;
03440 #endif
03441     Mtarget2 = target->GetMass2()*(1.-perCent); // @@ Isn't it too big reduction? M.K.
03442   }
03443   // @@ From this point it repeats the Diffractional excitation (? Use flag ?)
03444   // Transform momenta to cms and then rotate parallel to z axis;
03445   G4LorentzVector Psum=Pprojectile+Ptarget;
03446   G4LorentzRotation toCms(-Psum.boostVector()); // Boost Rotation to CMS
03447   G4LorentzVector Ptmp=toCms*Pprojectile;
03448   if(Ptmp.pz()<=0.) // "String" moving backwards in CMS, abort collision !! ? M.K.
03449   {
03450 #ifdef debug
03451     G4cout<<"G4QIonIonCollision::ExciteSingDiffParticipants: *1*abortCollision*1*"<<G4endl;
03452 #endif
03453     return false; 
03454   }         
03455   toCms.rotateZ(-Ptmp.phi());
03456   toCms.rotateY(-Ptmp.theta());
03457 #ifdef debug
03458   G4cout<<"G4QInel::ExciteSingDiffParticipantts: Be4Boost Pproj="<<Pprojectile<<", Ptarg="
03459         <<Ptarget<<G4endl;
03460 #endif
03461   G4LorentzRotation toLab(toCms.inverse()); // Boost Rotation to LabSys (LS)
03462   Pprojectile.transform(toCms);
03463   Ptarget.transform(toCms);
03464 #ifdef debug
03465   G4cout<< "G4QInelast::ExciteDiffParticipantts: AfterBoost Pproj="<<Pprojectile<<"Ptarg="
03466         <<Ptarget<<", cms4M="<<Pprojectile+Ptarget<<G4endl;
03467 
03468   G4cout<<"G4QInelast::ExciteDiffParticipantts: ProjX+="<<Pprojectile.plus()<<", ProjX-="
03469         <<Pprojectile.minus()<<", TargX+="<< Ptarget.plus()<<", TargX-="<<Ptarget.minus()
03470         <<G4endl;
03471 #endif
03472   G4LorentzVector Qmomentum(0.,0.,0.,0.);
03473   G4int whilecount=0;
03474   do
03475   {
03476     //  Generate pt  
03477     G4double maxPtSquare=sqr(Ptarget.pz());
03478     if(whilecount++>=500 && whilecount%100==0) // @@ M.K. Hardwired limits 
03479 #ifdef debug
03480     G4cout<<"G4QIonIonCollision::ExciteSingDiffParticipantts: can loop, loopCount="
03481           <<whilecount<<", maxPtSquare="<<maxPtSquare<<G4endl;
03482 #endif
03483     if(whilecount>1000)                        // @@ M.K. Hardwired limits 
03484     {
03485 #ifdef debug
03486       G4cout<<"G4QIonIonCollision::ExciteSingDiffParticipants:*2* abortLoop!! *2*"<<G4endl;
03487 #endif
03488       return false;    //  Ignore this interaction 
03489     }
03490     Qmomentum=G4LorentzVector(GaussianPt(widthOfPtSquare,maxPtSquare),0);
03491 #ifdef debug
03492     G4cout<<"G4QInelast::ExciteSingDiffParticipants: generated Pt="<<Qmomentum
03493           <<", ProjPt="<<Pprojectile+Qmomentum<<", TargPt="<<Ptarget-Qmomentum<<G4endl;
03494 #endif
03495     //  Momentum transfer
03496     G4double Xmin=0.;
03497     G4double Xmax=1.;
03498     G4double Xplus =ChooseX(Xmin,Xmax);
03499     G4double Xminus=ChooseX(Xmin,Xmax);
03500 #ifdef debug
03501     G4cout<<"G4QInel::ExciteSingDiffPartici: X-plus="<<Xplus<<", X-minus="<<Xminus<<G4endl;
03502 #endif
03503     G4double pt2=G4ThreeVector(Qmomentum.vect()).mag2();
03504     G4double Qplus =-pt2/Xminus/Ptarget.minus();
03505     G4double Qminus= pt2/Xplus /Pprojectile.plus();
03506     if (KeepProjectile)
03507       Qminus=(projectile->GetMass2()+pt2)/(Pprojectile.plus()+Qplus) - Pprojectile.minus();
03508     else Qplus=Ptarget.plus() - (target->GetMass2()+pt2)/(Ptarget.minus()-Qminus);  
03509     Qmomentum.setPz((Qplus-Qminus)/2);
03510     Qmomentum.setE( (Qplus+Qminus)/2);
03511 #ifdef debug
03512     G4cout<<"G4QInel::ExciteDiffParticip: Qplus="<<Qplus<<", Qminus="<<Qminus<<", pt2="
03513           <<pt2<<", Qmomentum="<<Qmomentum<<", ProjM="<<(Pprojectile+Qmomentum).mag()
03514           <<", TargM="<<(Ptarget-Qmomentum).mag()<<G4endl;
03515 #endif
03516     // while is different from the Double Diffractive Excitation (@@ !)
03517     //} while((Pprojectile+Qmomentum).mag2()<= Mprojectile2 ||
03518     //        (Ptarget-Qmomentum).mag2()<=Mtarget2);
03519   } while((Ptarget-Qmomentum).mag2()<=Mtarget2 ||
03520           (Pprojectile+Qmomentum).mag2()<=Mprojectile2 ||
03521           (Ptarget-Qmomentum).e() < 0. || (Pprojectile+Qmomentum).e() < 0.);
03522   Pprojectile += Qmomentum;
03523   Ptarget     -= Qmomentum;
03524 #ifdef debug
03525   G4cout<<"G4QIonIonCollision::ExciteSingDiffParticipan: Proj(Q)="<<Pprojectile<<"(E="
03526         <<Pprojectile.e()<<"), Targ(Q)="<<Ptarget<<"(E="<<Ptarget.e()
03527         <<"), Proj(back)="<<toLab*Pprojectile<<", Targ(bac)="<< toLab*Ptarget << G4endl;
03528 #endif
03529   // Transform back and update SplitableHadron Participant.
03530   Pprojectile.transform(toLab);
03531   Ptarget.transform(toLab);
03532 #ifdef debug
03533   G4cout<< "G4QIonIonCollision::ExciteSingleDiffParticipants: TgM="<<Ptarget.mag()<<G4endl;
03534 #endif
03535   target->Set4Momentum(Ptarget);  
03536 #ifdef debug
03537   G4cout<<"G4QIonIonCollision::ExciteSingleParticipants:ProjM="<<Pprojectile.mag()<<G4endl;
03538 #endif
03539   projectile->Set4Momentum(Pprojectile);
03540   return true;
03541 } // End of ExciteSingleDiffParticipants

G4QHadronVector * G4QIonIonCollision::Fragment (  ) 

Definition at line 1557 of file G4QIonIonCollision.cc.

References G4QEnvironment::AddQuasmon(), Breeder(), FatalException, G4QEnvironment::Fragment(), G4cerr, G4cout, G4endl, G4Exception(), G4QHadron::Get4Momentum(), G4QNucleus::GetA(), G4QPDGCode::GetMass(), G4QNucleus::GetPDG(), G4QHadron::GetPDGCode(), G4QHadron::GetQC(), G4QNucleus::GetQCZNS(), G4QCHIPSWorld::GetQParticle(), G4QChipolino::GetQPDG1(), G4QChipolino::GetQPDG2(), G4QContent::GetSPDGCode(), G4QNucleus::GetZ(), G4QParticle::MinMassOfFragm(), G4Quasmon::Set4Momentum(), and G4Quasmon::SetQC().

Referenced by G4QInelastic::PostStepDoIt().

01558 { // This is the member function fragmenting Strings & Quasmons (in nuclear matter)
01559 #ifdef debug
01560   G4cout<<"*******>G4QIonIonCollision::Fragment: ***Called***, Res="<<theResult<<G4endl;
01561 #endif
01562   G4int striNum=strings.size();                             // Find out if there're strings
01563   G4int hadrNum=theResult->size();                          // Find out if there're hadrons
01564 #ifdef edebug
01565   G4int nQm=theQuasmons.size();
01566   G4LorentzVector totLS4M=theProjNucleus.Get4Momentum()+theTargNucleus.Get4Momentum(); //LS
01567   G4int totChg=theProjNucleus.GetZ()+theTargNucleus.GetZ();
01568   G4int totBaN=theTargNucleus.GetA()+theTargNucleus.GetA();
01569   G4cout<<"-EMCLS-G4QIn::Fragment: CHECKRecovery, #ofS="<<striNum<<", Nuc4M(E=M)="<<totLS4M
01570         <<",#Q="<<nQm<<",#H="<<hadrNum<<G4endl;
01571   for(G4int i=0; i < striNum; i++)
01572   {
01573     G4LorentzVector strI4M=strings[i]->Get4Momentum();
01574     totLS4M+=strI4M;
01575     G4int sChg=strings[i]->GetCharge();
01576     totChg+=sChg;
01577     G4int sBaN=strings[i]->GetBaryonNumber();
01578     totBaN+=sBaN;
01579     G4cout<<"-EMCLS-G4QIonIonCollision::Fragment:String#"<<i<<", 4M="<<strI4M<<", M="
01580           <<strI4M.m()<<", C="<<sChg<<", B="<<sBaN<<G4endl;
01581   }
01582   for(G4int i=0; i < nQm; i++)
01583   {
01584     G4LorentzVector hI4M=theQuasmons[i]->Get4Momentum();
01585     totLS4M+=hI4M;
01586     G4int hChg=theQuasmons[i]->GetCharge();
01587     totChg+=hChg;
01588     G4int hBaN=theQuasmons[i]->GetBaryonNumber();
01589     totBaN+=hBaN;
01590     G4cout<<"-EMCLS-G4QIonIonCollision::Fragment: Quasmon#"<<i<<", 4M="<<hI4M<<", C="<<hChg
01591           <<", B="<<hBaN<<G4endl;
01592   }
01593   for(G4int i=0; i < hadrNum; i++)
01594   {
01595     G4LorentzVector hI4M=(*theResult)[i]->Get4Momentum();
01596     totLS4M+=hI4M;
01597     G4int hChg=(*theResult)[i]->GetCharge();
01598     totChg+=hChg;
01599     G4int hBaN=(*theResult)[i]->GetBaryonNumber();
01600     totBaN+=hBaN;
01601     G4cout<<"-EMCLS-G4QIn::Fragment:H#"<<i<<",4M="<<hI4M<<",C="<<hChg<<",B="<<hBaN<<G4endl;
01602   }
01603 #endif
01604 #ifdef debug
01605   G4cout<<"***>G4QIonIonCollision::Fragm: #OfStr="<<striNum<<", #OfRes="<<hadrNum<<G4endl;
01606 #endif
01607   if(!striNum && hadrNum)                                   // Quasi-elastic or decoupled p
01608   {
01609 #ifdef debug
01610     G4cout<<"***>G4QIonIonCollision::Fragm:**Quasi-Elastic**, #OfResult="<<hadrNum<<G4endl;
01611 #endif
01612     return theResult;
01613   }
01614   else if(striNum) Breeder();                               // Strings fragmentation
01615   else                                                      // No strings, make HadrNucleus
01616   {
01617     if(hadrNum)
01618     {
01619       for(G4int ih=0; ih<hadrNum; ih++) delete (*theResult)[ih];
01620       theResult->clear();
01621     }
01622     G4LorentzVector rp4M=theProjNucleus.Get4Momentum();     // Nucleus 4-momentum in LS
01623     G4int rpPDG=theProjNucleus.GetPDG();                    // Nuclear PDG
01624     G4QHadron* resPNuc = new G4QHadron(rpPDG,rp4M);         // Nucleus -> Hadron
01625     theResult->push_back(resPNuc);                          // Fill the residual nucleus
01626     G4LorentzVector rt4M=theTargNucleus.Get4Momentum();     // Nucleus 4-momentum in LS
01627     G4int rtPDG=theTargNucleus.GetPDG();                    // Nuclear PDG
01628     G4QHadron* resTNuc = new G4QHadron(rtPDG,rt4M);         // Nucleus -> Hadron
01629     theResult->push_back(resTNuc);                          // Fill the residual nucleus
01630   }
01631   G4int nQuas=theQuasmons.size();                           // Size of the Quasmon OUTPUT
01632   G4int theRS=theResult->size();                            // Size of Hadron Output by now
01633 #ifdef debug
01634   G4cout<<"***>G4QIonIonCollision::Fragm:beforeEnv, #OfQ="<<nQuas<<",#OfR="<<theRS<<G4endl;
01635 #endif
01636   if(nQuas && theRS)
01637   {
01638 
01639     G4QHadron* resNuc = (*theResult)[theRS-1];              // Pointer to Residual Nucleus
01640     G4LorentzVector resNuc4M = resNuc->Get4Momentum();      // 4-Momentum of the Nucleuz
01641     G4int           resNucPDG= resNuc->GetPDGCode();        // PDG Code of the Nucleus
01642     G4QNucleus      theEnv(resNucPDG);                      // NucleusHadron->NucleusAtRest
01643     delete resNuc;                                          // Delete resNucleus as aHadron
01644     theResult->pop_back();                                  // Exclude the nucleus from HV
01645     --theRS;                                                // Reduce the OUTPUT by theNucl
01646 #ifdef pdebug
01647     G4cout<<"G4QIonIonCollision::Fragm:#OfRemainingHadron="<<theRS<<", A="<<theEnv<<G4endl;
01648 #endif
01649     // Now we need to be sure that the compound nucleus is heavier than the Ground State
01650     for(G4int j=theRS-1; j>-2; --j)                         // try to reach M_compound>M_GS
01651     {
01652       G4LorentzVector qsum4M=resNuc4M;                      // Proto compound 4-momentum
01653       G4QContent qsumQC=theEnv.GetQCZNS();                  // Proto compound Quark Content
01654 #ifdef pdebug
01655       G4cout<<"G4QIonIonCollis::Fragm: rN4M"<<qsum4M<<qsum4M.m()<<",rNQC="<<qsumQC<<G4endl;
01656 #endif
01657       G4Quasmon* firstQ=0;                                  // Prototype of theFirstQuasmon
01658       G4LorentzVector first4M;                              // Proto of the FirstQuasmon 4M
01659       G4QContent firstQC;                                   // Proto of the FirstQuasmon QC
01660       for(G4int i=0; i<nQuas; ++i)                          // LOOP over Quasmons
01661       {
01662         G4Quasmon* curQuasm=theQuasmons[i];                 // current Quasmon
01663         G4LorentzVector cur4M=curQuasm->Get4Momentum();     // 4-Mom of the Quasmon
01664         G4QContent curQC=curQuasm->GetQC();                 // Quark Content of the Quasmon
01665         qsum4M+=cur4M;                                      // Add quasmon's 4-momentum
01666         qsumQC+=curQC;                                      // Add quasmon's Quark Content
01667 #ifdef pdebug
01668         G4cout<<"G4QIonIonCollis::Fragm: Q#"<<i<<", QQC="<<curQC<<", sQC="<<qsumQC<<G4endl;
01669 #endif
01670         if(!i)                                              // Remember 1-st for correction
01671         {
01672           firstQ =curQuasm;
01673           first4M=cur4M;
01674           firstQC=curQC;
01675         }
01676       }
01677       G4int miPDG=qsumQC.GetSPDGCode();                     // PDG of minM of hadron/fragm.
01678       G4double gsM=0.;                                      // Proto minM of had/frag forQC
01679       if(miPDG == 10)
01680       {
01681         G4QChipolino QCh(qsumQC);                           // define TotNuc as a Chipolino
01682         gsM=QCh.GetQPDG1().GetMass()+QCh.GetQPDG2().GetMass(); // Sum of Hadron Masses
01683         //gsM=theWorld->GetQParticle(QCh.GetQPDG1())->MinMassOfFragm() +
01684         //    theWorld->GetQParticle(QCh.GetQPDG2())->MinMassOfFragm();
01685       }
01686       // @@ it is not clear, why it does not work ?!
01687       //else if(miPDG>80000000)                             // Compound Nucleus
01688       //{
01689       //  G4QNucleus rtN(qsumQC);                           // CreatePseudoNucl for totComp
01690       //  gsM=rtN.GetGSMass();                              // MinMass of residQ+(Env-ParC)
01691       //}
01692       else if(miPDG < 80000000 && std::abs(miPDG)%10 > 2)
01693                            gsM=theWorld->GetQParticle(G4QPDGCode(miPDG))->MinMassOfFragm();
01694       else gsM=G4QPDGCode(miPDG).GetMass();                 // minM of hadron/fragm. for QC
01695       G4double reM=qsum4M.m();                              // real mass of the compound
01696 #ifdef pdebug
01697       G4cout<<"G4QIonIonCollision::Fragm: PDG="<<miPDG<<", rM="<<reM<<",GSM="<<gsM<<G4endl;
01698 #endif
01699       if(reM > gsM) break;                                  // CHIPS can be called
01700       if(j > -1)                                            // Can try to add hadrons to Q0
01701       {
01702         G4QHadron* cH = (*theResult)[j];                    // Pointer to the last Hadron
01703         G4LorentzVector h4M = cH->Get4Momentum();           // 4-Momentum of the Hadron
01704         G4QContent      hQC = cH->GetQC();                  // QC of the Hadron
01705         firstQ->Set4Momentum(first4M+h4M);                  // Update the Quasmon's 4-Mom
01706         firstQ->SetQC(firstQC+hQC);                         // Update the Quasmon's QCont
01707         delete cH;                                          // Delete the Hadron
01708         theResult->pop_back();                              // Exclude the hadron from HV
01709 #ifdef pdebug
01710         G4cout<<"G4QInel::Fragm: H#"<<j<<", hQC="<<hQC<<",hPDG="<<cH->GetPDGCode()<<G4endl;
01711 #endif
01712       }
01713       else
01714       {
01715         G4cerr<<"***G4QIonIonCollision::Fra:PDG="<<miPDG<<",M="<<reM<<",GSM="<<gsM<<G4endl;
01716         G4Exception("G4QIonIonCollision::Fragment:","27",FatalException,"Can'tRecoverGSM");
01717       }
01718     }
01719     G4double nucE=resNuc4M.e();                             // Total energy of the nuclEnv
01720     if(nucE<1.E-12) nucE=0.;                                // Computer accuracy safety
01721     G4ThreeVector   nucVel(0.,0.,0.);                       // Proto of the NucleusVelocity
01722     G4QHadronVector* output=0;                              // NucleusFragmentation Hadrons
01723     G4QEnvironment* pan= new G4QEnvironment(theEnv);        // ---> DELETED --->----------+
01724 #ifdef pdebug
01725     G4cout<<"G4QIonIonCollision::Fragm: nucE="<<nucE<<", nQ="<<nQuas<<G4endl; //          |
01726 #endif
01727     if(nucE) nucVel=resNuc4M.vect()/nucE;                   // The NucleusVelocity        |
01728     for(G4int i=0; i<nQuas; ++i)                            // LOOP over Quasmons         |
01729     {                                                       //                            |
01730       G4Quasmon* curQuasm=theQuasmons[i];                   // current Quasmon            |
01731       if(nucE) curQuasm->Boost(-nucVel);                    // Boost it to CMS of Nucleus |
01732       pan->AddQuasmon(curQuasm);                            // Fill the predefined Quasmon|
01733 #ifdef pdebug
01734       G4LorentzVector cQ4M=curQuasm->Get4Momentum();        // Just for printing          |
01735       G4cout<<"G4QIonIonCollision::Fragment: Quasmon# "<<i<<" added, 4M="<<cQ4M<<G4endl;//|
01736 #endif
01737     }                                                       //                            |
01738     try                                                     //                            |
01739     {                                                       //                            |
01740       delete output;                                        //                            |
01741       output = pan->Fragment();// DESTROYED after theHadrons are transferred to theResult |
01742     }                                                       //                          | |
01743     catch (G4QException& error)                             //                          | |
01744     {                                                       //                          | |
01745       G4cerr<<"***G4QIonIonCollision::Fragment: G4QE Exception is catched"<<G4endl; //  | |
01746       // G4Exception("G4QIonIonCollision::Fragment:","27",FatalException,"CHIPSCrash");//  | |
01747       G4Exception("G4QIonIonCollision::Fragment()", "HAD_CHPS_0027",
01748                   FatalException, "CHIPSCrash");
01749     }                                                       //                          | |
01750     delete pan;                              // Delete the Nuclear Environment <-----<--+-+
01751     if(output)                               // Output exists                           |
01752     {                                        //                                         |
01753       G4int nOut=output->size();             // #ofHadrons in the Nuclear Fragmentation |
01754       for(G4int j=0; j<nOut; j++)            // LOOP over Hadrons transferring to LS    |
01755       {                                      //                                         |
01756         G4QHadron* curHadron=(*output)[j];   // Hadron from the nucleus fragmentation   |
01757         if(nucE) curHadron->Boost(nucVel);   // Boost it back to Laboratory System      |
01758         theResult->push_back(curHadron);     // Transfer it to the result               |
01759       }                                      //                                         |
01760       delete output;                         // Delete the OUTPUT <-----<-----<-----<---+
01761     }
01762   }
01763   else if(!striNum) G4cout<<"-Warning-G4QIonIonCollision::Fragment:NothingWasDone"<<G4endl;
01764 #ifdef debug
01765   G4cout<<"=--=>G4QIonIonCollision::Fragment: Final #OfResult="<<theResult->size()<<G4endl;
01766 #endif
01767     G4int nQ =theQuasmons.size();
01768     if(nQ) theQuasmons.clear();                              // @@ Not necesary ?
01769 #ifdef edebug
01770     G4LorentzVector f4M(0.,0.,0.,0.);                        // Sum of the Result in LS
01771     G4int fCh=totChg;
01772     G4int fBN=totBaN;
01773     G4int nHd=theResult->size();
01774     G4cout<<"-EMCLS-G4QIonIonCollision::Fragment: #ofHadr="<<nHd<<",#OfQuasm="<<nQ<<G4endl;
01775     for(G4int i=0; i<nHd; i++)
01776     {
01777       G4LorentzVector hI4M=(*theResult)[i]->Get4Momentum();
01778       f4M+=hI4M;
01779       G4int hChg=(*theResult)[i]->GetCharge();
01780       fCh-=hChg;
01781       G4int hBaN=(*theResult)[i]->GetBaryonNumber();
01782       fBN-=hBaN;
01783       G4cout<<"-EMCLS-G4QIonIonCollision::Fragment: Hadron#"<<i<<", 4M="<<hI4M<<", PDG="
01784             <<(*theResult)[i]->GetPDGCode()<<", C="<<hChg<<", B="<<hBaN<<G4endl;
01785     }
01786     G4cout<<"-EMCLS-G4QInel::Fragm:, r4M="<<f4M-totLS4M<<", rC="<<fCh<<",rB="<<fBN<<G4endl;
01787 #endif
01788   return theResult;
01789 } // End of fragmentation

G4ThreeVector G4QIonIonCollision::GaussianPt ( G4double  widthSquare,
G4double  maxPtSquare 
) const [protected]

Definition at line 3572 of file G4QIonIonCollision.cc.

References G4cout, G4endl, G4UniformRand, and sqr().

Referenced by ExciteDiffParticipants(), and ExciteSingDiffParticipants().

03573 {
03574 #ifdef debug
03575   G4cout<<"G4QIonIonCollision::GaussianPt: w^2="<<widthSq<<",maxPt2="<<maxPtSquare<<G4endl;
03576 #endif
03577   G4double pt2=0.;
03578   G4double rm=maxPtSquare/widthSq;                      // Negative
03579   if(rm>-.01) pt2=widthSq*(std::sqrt(1.-G4UniformRand()*(1.-sqr(1.+rm)))-1.);
03580   else        pt2=widthSq*std::log(1.-G4UniformRand()*(1.-std::exp(rm)));
03581   pt2=std::sqrt(pt2);                                   // It is not pt2, but pt now
03582   G4double phi=G4UniformRand()*twopi;
03583   return G4ThreeVector(pt2*std::cos(phi),pt2*std::sin(phi),0.);    
03584 } // End of GaussianPt

G4bool G4QIonIonCollision::IsSingleDiffractive (  )  [inline, protected]

Definition at line 76 of file G4QIonIonCollision.hh.

References G4UniformRand.

Referenced by G4QIonIonCollision().

00076 {G4bool r=false; if(G4UniformRand()<1.) r=true; return r;}

std::pair< G4int, G4int > G4QIonIonCollision::ReducePair ( G4int  P1,
G4int  P2 
) const [protected]

Definition at line 3651 of file G4QIonIonCollision.cc.

References G4cout, and G4endl.

Referenced by Breeder().

03652 {
03653 #ifdef debug
03654   G4cout<<"G4QIonIonCollision::ReducePair: **Called** P1="<<P1<<", P2="<<P2<<G4endl;
03655 #endif
03656   G4int P11=P1/10;
03657   G4int P12=P1%10;
03658   G4int P21=P2/10;
03659   G4int P22=P2%10;
03660   if     (P11==P21) return std::make_pair(P12,P22);
03661   else if(P11==P22) return std::make_pair(P12,P21);
03662   else if(P12==P21) return std::make_pair(P11,P22);
03663   else if(P12==P22) return std::make_pair(P11,P21);
03664   //#ifdef debug
03665   G4cout<<"-Warning-G4QIonIonCollision::ReducePair:**Failed**,P1="<<P1<<",P2="<<P2<<G4endl;
03666   //#endif
03667   return std::make_pair(0,0);                         // Reduction failed
03668 }

void G4QIonIonCollision::SetParameters ( G4int  nC,
G4double  strTens,
G4double  tubeDens,
G4double  SigPt 
) [static]

Definition at line 3543 of file G4QIonIonCollision.cc.

03544 {
03545   nCutMax            = nC;             // max number of pomeron cuts
03546   stringTension      = stT;            // string tension for absorbed energy
03547   tubeDensity        = tbD;            // Flux Tube Density of nuclear nucleons
03548   widthOfPtSquare    = -2*SigPt*SigPt; // width^2 of pt for string excitation
03549 }

G4int G4QIonIonCollision::SumPartonPDG ( G4int  PDG1,
G4int  PFG2 
) const [protected]

Definition at line 3586 of file G4QIonIonCollision.cc.

References FatalException, G4cerr, G4endl, and G4Exception().

Referenced by Breeder().

03587 {
03588   if      (PDG1 < 7 && PDG1 > 0 && PDG2 < 7 && PDG2 > 0) // Sum up two Q in DiQ (S=0)
03589   {
03590     if(PDG1 > PDG2) return PDG1*1000+PDG2*100+1;
03591     else            return PDG2*1000+PDG1*100+1;
03592   }
03593   else if (PDG1 >-7 && PDG1 < 0 && PDG2 >-7 && PDG2 < 0) // Sum up two AQ in ADiQ (S=0)
03594   {
03595     if(-PDG1 > -PDG2) return PDG1*1000+PDG2*100-1;
03596     else              return PDG2*1000+PDG1*100-1;
03597   }
03598   else if (PDG1 <-99 && PDG2 < 7 && PDG2 > 0) // Sum up Q and ADiQ in AQ
03599   {
03600     G4int PDG=-PDG1/100;
03601     if(PDG2==PDG/10) return -PDG%10;
03602     if(PDG2==PDG%10) return -PDG/10;
03603     else
03604     {
03605       G4cerr<<"***4QFragmentation::SumPartonPDG: PDG1="<<PDG1<<", PDG2="<<PDG2<<G4endl;
03606       G4Exception("G4QIonIonCollision::SumPartonPDG:","72",FatalException,"Q&ADiQnoMatch");
03607     }
03608   }
03609   else if (PDG2 <-99 && PDG1 < 7 && PDG1 > 0) // Sum up ADiQ and Q in AQ
03610   {
03611     G4int PDG=-PDG2/100;
03612     if(PDG1==PDG/10) return -PDG%10;
03613     if(PDG1==PDG%10) return -PDG/10;
03614     else
03615     {
03616       G4cerr<<"***4QFragmentation::SumPartonPDG: PDG1="<<PDG1<<", PDG2="<<PDG2<<G4endl;
03617       G4Exception("G4QIonIonCollision::SumPartonPDG:","72",FatalException,"ADiQ&QnoMatch");
03618     }
03619   }
03620   else if (PDG1 > 99 && PDG2 >-7 && PDG2 < 0) // Sum up DiQ and AQ in Q
03621   {
03622     G4int PDG=PDG1/100;
03623     if(PDG2==-PDG/10) return PDG%10;
03624     if(PDG2==-PDG%10) return PDG/10;
03625     else
03626     {
03627       G4cerr<<"***4QFragmentation::SumPartonPDG: PDG1="<<PDG1<<", PDG2="<<PDG2<<G4endl;
03628       G4Exception("G4QIonIonCollision::SumPartonPDG:","72",FatalException,"DiQ&AQnoMatch");
03629     }
03630   }
03631   else if (PDG2 > 99 && PDG1 >-7 && PDG1 < 0) // Sum up AQ and DiQ in Q
03632   {
03633     G4int PDG=PDG2/100;
03634     if(PDG1==-PDG/10) return PDG%10;
03635     if(PDG1==-PDG%10) return PDG/10;
03636     else
03637     {
03638       G4cerr<<"***G4QIonIonCollision::SumPartonPDG: PDG1="<<PDG1<<", PDG2="<<PDG2<<G4endl;
03639       G4Exception("G4QIonIonCollision::SumPartonPDG:","72",FatalException,"AQ&DiQnoMatch");
03640     }
03641   }
03642   else
03643   {
03644     G4cerr<<"***G4QIonIonCollision::SumPartonPDG: PDG1="<<PDG1<<", PDG2="<<PDG2<<G4endl;
03645     G4Exception("G4QIonIonCollision::SumPartonPDG:","72",FatalException,"Can'tSumUpParts");
03646   }
03647   return 0;
03648 } // End of SumPartonPDG

void G4QIonIonCollision::SwapPartons (  )  [protected]

Definition at line 3903 of file G4QIonIonCollision.cc.

References DBL_MAX, G4cout, G4endl, G4QParton::Get4Momentum(), G4QParton::GetPDGCode(), and G4QParton::GetType().

Referenced by G4QIonIonCollision().

03904 {
03905   static const G4double baryM=800.;                  // Mass excess for baryons
03906   G4QStringVector::iterator ist;
03907   for(ist = strings.begin(); ist < strings.end(); ist++)
03908   {
03909     G4LorentzVector cS4M=(*ist)->Get4Momentum();
03910     G4double cSM2=cS4M.m2();                         // Squared mass of the String
03911     if(cSM2<0.1)                                     // Parton Swapping is needed
03912     {
03913       G4QParton* cLeft=(*ist)->GetLeftParton();      // Current String Left Parton 
03914       G4QParton* cRight=(*ist)->GetRightParton();    // Current String Right Parton 
03915       G4int cLPDG=cLeft->GetPDGCode();
03916       G4int cRPDG=cRight->GetPDGCode();
03917       G4LorentzVector cL4M=cLeft->Get4Momentum();
03918       G4LorentzVector cR4M=cRight->Get4Momentum();
03919       G4int cLT=cLeft->GetType();
03920       G4int cRT=cRight->GetType();
03921       G4QStringVector::iterator sst;                 // Selected partner string
03922       G4QStringVector::iterator pst;                 // LOOP iterator
03923       G4double maxM=-DBL_MAX;                        // Swapping providing the max mass
03924       G4int    sDir=0;                               // Selected direction of swapping
03925       for(pst = strings.begin(); pst < strings.end(); pst++) if(pst != ist)
03926       {
03927         G4QParton* pLeft=(*pst)->GetLeftParton();    // Partner String Left Parton 
03928         G4QParton* pRight=(*pst)->GetRightParton();  // Partner String Right Parton 
03929         G4int pLPDG=pLeft->GetPDGCode();
03930         G4int pRPDG=pRight->GetPDGCode();
03931         G4LorentzVector pL4M=pLeft->Get4Momentum();
03932         G4LorentzVector pR4M=pRight->Get4Momentum();
03933         G4int pLT=pLeft->GetType();
03934         G4int pRT=pRight->GetType();
03935         G4double LM=0.;
03936         G4double RM=0.;
03937         if(((cLPDG<-7 || (cLPDG>0 && cLPDG< 7) ) && (pLPDG<-7 || (pLPDG>0 && pLPDG< 7) ))||
03938            ((cLPDG> 7 || (cLPDG<0 && cLPDG>-7) ) && (pLPDG> 7 || (pLPDG<0 && cLPDG>-7) )))
03939         {
03940           G4double pLM2=(cL4M+pR4M).m2();                      // new partner M2
03941           G4double cLM2=(cR4M+pL4M).m2();                      // new partner M2
03942           if(pLM2>0. && cLM2>0.)
03943           {
03944             G4double pLM=std::sqrt(pLM2);
03945             if(cLT+pRT==3) pLM-=baryM;
03946             G4double cLM=std::sqrt(cLM2);
03947             if(cRT+pLT==3) cLM-=baryM;
03948             LM=std::min(pLM2,cLM2);
03949           }
03950         }
03951         if(((cRPDG<-7 || (cRPDG>0 && cRPDG< 7) ) && (pRPDG<-7 || (pRPDG>0 && pRPDG< 7) ))||
03952            ((cRPDG> 7 || (cRPDG<0 && cRPDG>-7) ) && (pRPDG> 7 || (pRPDG<0 && cRPDG>-7) )) )
03953         {
03954           G4double pRM2=(cR4M+pL4M).m2();                      // new partner M2
03955           G4double cRM2=(cL4M+pR4M).m2();                      // new partner M2
03956           if(pRM2>0. && cRM2>0.)
03957           {
03958             G4double pRM=std::sqrt(pRM2);
03959             if(cRT+pLT==3) pRM-=baryM;
03960             G4double cRM=std::sqrt(cRM2);
03961             if(cLT+pRT==3) cRM-=baryM;
03962             RM=std::min(pRM,cRM);
03963           }
03964         }
03965         G4int dir=0;
03966         G4double sM=0.;
03967         if( LM && LM > RM )
03968         {
03969           dir= 1;
03970           sM=LM;
03971         }
03972         else if(RM)
03973         {
03974           dir=-1;
03975           sM=RM;
03976         }
03977         if(sM > maxM)
03978         {
03979           sst=pst;
03980           maxM=sM;
03981           sDir=dir;
03982         }
03983       }
03984       if(sDir)
03985       {
03986         G4QParton* pLeft=(*sst)->GetLeftParton();    // Partner String Left Parton 
03987         G4QParton* pRight=(*sst)->GetRightParton();  // Partner String Right Parton 
03988         G4QParton* swap=pLeft;                       // Prototype remember the partner Left
03989         if(sDir>0)                                   // swap left partons
03990         {
03991           (*sst)->SetLeftParton(cLeft);
03992           (*ist)->SetLeftParton(swap);
03993         }
03994         else
03995         {
03996           swap=pRight;
03997           (*sst)->SetRightParton(cRight);
03998           (*ist)->SetRightParton(swap);
03999         }
04000       }
04001 #ifdef debug
04002       else G4cout<<"***G4QIonIonCollision::SwapPartons:**Failed**,cLPDG="<<cLPDG<<",cRPDG="
04003                  <<cRPDG<<",-->cM2="<<cSM2<<G4endl;
04004 #endif
04005       
04006     }
04007   }
04008 }


The documentation for this class was generated from the following files:
Generated on Mon May 27 17:53:09 2013 for Geant4 by  doxygen 1.4.7