G4QIonIonCollision.cc

Go to the documentation of this file.
00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00026 //
00027 // $Id$
00028 //
00029 // -----------------------------------------------------------------------------
00030 //      GEANT 4 class header file
00031 //
00032 //                 History: 
00033 //     Created by Mikhail Kossov, October 2006
00034 //     CHIPS QGS fragmentation class 
00035 //     For comparison similar member functions can be found in the G4 classes:
00036 //     G4QGSParticipants
00037 //     G4QGSModels
00038 //     G4ExcitedStringDecay
00039 // -----------------------------------------------------------------------------
00040 // Short description: CHIPS QG string fragmentation class
00041 // Rhe key member function is Scatter, making the interaction (see G4QCollision)
00042 // -----------------------------------------------------------------------------
00043 //#define debug
00044 //#define edebug
00045 //#define pdebug
00046 //#define sdebug
00047 //#define ppdebug
00048 
00049 #include "G4QIonIonCollision.hh"
00050 #include "G4PhysicalConstants.hh"
00051 #include "G4SystemOfUnits.hh"
00052 
00053 // Promoting model parameters from local variables class properties @@(? M.K.)
00054 
00055 // Definition of static parameters
00056 G4int    G4QIonIonCollision::nCutMax=7; 
00057 G4double G4QIonIonCollision::stringTension=1.*GeV/fermi;
00058 G4double G4QIonIonCollision::tubeDensity  =1./fermi;
00059 // Parameters of diffractional fragmentation
00060 G4double G4QIonIonCollision::widthOfPtSquare=-0.75*GeV*GeV; // ptWidth2 forStringExcitation
00061 
00062 G4QIonIonCollision::G4QIonIonCollision(G4QNucleus &pNucleus, const G4QNucleus &tNucleus)
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 }
01551 
01552 G4QIonIonCollision::~G4QIonIonCollision()
01553 {
01554   std::for_each(strings.begin(), strings.end(), DeleteQString() );
01555 }
01556 
01557 G4QHadronVector* G4QIonIonCollision::Fragment()
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
01790 
01791 void G4QIonIonCollision::Breeder()
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
03306 
03307 // Excite double diffractive string
03308 G4bool G4QIonIonCollision::ExciteDiffParticipants(G4QHadron* projectile,
03309                                                 G4QHadron* target) const
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
03412 
03413 
03414 // Excite single diffractive string
03415 G4bool G4QIonIonCollision::ExciteSingDiffParticipants(G4QHadron* projectile,
03416                                                     G4QHadron* target) const
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
03542 
03543 void G4QIonIonCollision::SetParameters(G4int nC,G4double stT, G4double tbD, G4double SigPt)
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 }
03550 
03551 G4double G4QIonIonCollision::ChooseX(G4double Xmin, G4double Xmax) const
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
03570 
03571 // Add CHIPS exponential Pt distribution (see Fragmentation)
03572 G4ThreeVector G4QIonIonCollision::GaussianPt(G4double widthSq, G4double maxPtSquare) const
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
03585 
03586 G4int G4QIonIonCollision::SumPartonPDG(G4int PDG1, G4int PDG2) const
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
03649 
03650 // Reduces quark pairs (unsigned 2 digits) to quark singles (unsigned)
03651 std::pair<G4int,G4int> G4QIonIonCollision::ReducePair(G4int P1, G4int P2) const
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 }
03669 
03670 // Select LL/RR (1) or LR/RL (-1) annihilation order (0, if the annihilation is impossible)
03671 G4int G4QIonIonCollision::AnnihilationOrder(G4int LS, G4int MPS, G4int uPDG, G4int mPDG,
03672                                           G4int sPDG, G4int nPDG) // ^L^^ Partner^^R^
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 }
03902 
03903 void G4QIonIonCollision::SwapPartons() // Swap string partons, if a string has negative M2
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 }

Generated on Mon May 27 17:49:35 2013 for Geant4 by  doxygen 1.4.7