G4QGSParticipants.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 #include <utility>
00027 
00028 #include "G4QGSParticipants.hh"
00029 #include "globals.hh"
00030 #include "G4SystemOfUnits.hh"
00031 #include "G4LorentzVector.hh"
00032 
00033 // Class G4QGSParticipants 
00034 
00035 // HPW Feb 1999
00036 // Promoting model parameters from local variables class properties
00037 G4int G4QGSParticipants_NPart = 0;
00038 
00039 G4QGSParticipants::G4QGSParticipants() : theDiffExcitaton(), //0.7*GeV, 250*MeV, 250*MeV),
00040                 ModelMode(SOFT),
00041                 //nCutMax(7),ThresholdParameter(0.45*GeV),
00042                 nCutMax(7),ThresholdParameter(0.000*GeV),
00043                 QGSMThreshold(3*GeV),theNucleonRadius(1.5*fermi)
00044 
00045 {
00046 }
00047 
00048 G4QGSParticipants::G4QGSParticipants(const G4QGSParticipants &right)
00049 : G4VParticipants(),ModelMode(right.ModelMode), nCutMax(right.nCutMax),
00050   ThresholdParameter(right.ThresholdParameter), QGSMThreshold(right.QGSMThreshold),
00051   theNucleonRadius(right.theNucleonRadius)
00052 {
00053 }
00054 
00055 
00056 G4QGSParticipants::~G4QGSParticipants()
00057 {
00058 }
00059 
00060 void G4QGSParticipants::BuildInteractions(const G4ReactionProduct  &thePrimary) 
00061 {
00062 
00063         // Find the collisions and collition conditions
00064         G4VSplitableHadron* aProjectile = SelectInteractions(thePrimary);
00065 
00066         // now build the parton pairs. HPW
00067         SplitHadrons();
00068 
00069         // soft collisions first HPW, ordering is vital
00070         PerformSoftCollisions();
00071 
00072         // the rest is diffractive HPW
00073         PerformDiffractiveCollisions();
00074 
00075         // clean-up, if necessary
00076         std::for_each(theInteractions.begin(), theInteractions.end(), DeleteInteractionContent());
00077         theInteractions.clear();
00078         std::for_each(theTargets.begin(), theTargets.end(), DeleteSplitableHadron());
00079         theTargets.clear();
00080         delete aProjectile;
00081 }
00082 
00083 G4VSplitableHadron* G4QGSParticipants::SelectInteractions(const G4ReactionProduct  &thePrimary)
00084 {
00085         G4VSplitableHadron* aProjectile = new G4QGSMSplitableHadron(thePrimary, TRUE); // @@@ check the TRUE
00086         G4PomeronCrossSection theProbability(thePrimary.GetDefinition()); // @@@ should be data member
00087         G4double outerRadius = theNucleus->GetOuterRadius();
00088         // Check reaction threshold
00089 
00090         theNucleus->StartLoop();
00091         G4Nucleon * pNucleon = theNucleus->GetNextNucleon();
00092         G4LorentzVector aPrimaryMomentum(thePrimary.GetMomentum(), thePrimary.GetTotalEnergy());
00093         //--DEBUG--G4cout << " qgspart- " << aPrimaryMomentum << " # " << aPrimaryMomentum.mag()
00094         //--DEBUG--      << pNucleon->Get4Momentum() << " # " << (pNucleon->Get4Momentum()).mag()<< G4endl;
00095         G4double s_nucleus = (aPrimaryMomentum + pNucleon->Get4Momentum()).mag2();
00096         G4double ThresholdMass = thePrimary.GetMass() + pNucleon->GetDefinition()->GetPDGMass();
00097         ModelMode = SOFT;
00098         if (sqr(ThresholdMass + ThresholdParameter) > s_nucleus)
00099         {
00100                 ModelMode = DIFFRACTIVE;
00101                 //throw G4HadronicException(__FILE__, __LINE__, "Initial energy is too low. The 4-vectors of the input are inconsistant with the particle masses.");
00102         }
00103         if (sqr(ThresholdMass + QGSMThreshold) > s_nucleus) // thus only diffractive in cascade!
00104         {
00105                 ModelMode = DIFFRACTIVE;
00106         }
00107 
00108         // first find the collisions HPW
00109         std::for_each(theInteractions.begin(), theInteractions.end(), DeleteInteractionContent());
00110         theInteractions.clear();
00111         G4int totalCuts = 0;
00112 
00113         #ifdef debug_QGS
00114                 G4double eK = thePrimary.GetKineticEnergy()/GeV;
00115         #endif
00116         #ifdef debug_G4QGSParticipants
00117                 G4double impactUsed = 0;
00118                 G4LorentzVector intNuclMom;
00119         #endif
00120 
00121         while(theInteractions.size() == 0)
00122         {
00123                 // choose random impact parameter HPW
00124                 std::pair<G4double, G4double> theImpactParameter;
00125                 theImpactParameter = theNucleus->ChooseImpactXandY(outerRadius+theNucleonRadius);
00126                 G4double impactX = theImpactParameter.first;
00127                 G4double impactY = theImpactParameter.second;
00128 
00129                 // loop over nucleons to find collisions
00130                 theNucleus->StartLoop();
00131                 G4int nucleonCount = 0; // debug
00132                 G4QGSParticipants_NPart = 0;
00133                 #ifdef debug_G4QGSParticipants
00134                         intNuclMom=aPrimaryMomentum;
00135                 #endif
00136                 while( (pNucleon = theNucleus->GetNextNucleon()) )
00137                 {
00138                         if(totalCuts>1.5*thePrimary.GetKineticEnergy()/GeV)
00139                         {
00140                                 break;
00141                         }
00142                         nucleonCount++; // debug
00143                         // Needs to be moved to Probability class @@@
00144                         G4LorentzVector nucleonMomentum=pNucleon->Get4Momentum();
00145                         nucleonMomentum.setE(nucleonMomentum.e()-pNucleon->GetBindingEnergy());
00146                         G4double s_nucleon = (aPrimaryMomentum + nucleonMomentum).mag2();
00147                         G4double Distance2 = sqr(impactX - pNucleon->GetPosition().x()) +
00148                                         sqr(impactY - pNucleon->GetPosition().y());
00149                         G4double Probability = theProbability.GetInelasticProbability(s_nucleon, Distance2);
00150                         // test for inelastic collision
00151                         G4double rndNumber = G4UniformRand();
00152                         //      ModelMode = DIFFRACTIVE;
00153                         if (Probability > rndNumber)
00154                         {
00155                         #ifdef debug_G4QGSParticipants
00156                                 G4cout << "DEBUG p="<< Probability<<" r="<<rndNumber<<" d="<<std::sqrt(Distance2)<<G4endl;
00157                                 G4cout << " qgspart+  " << aPrimaryMomentum << " # " << aPrimaryMomentum.mag()
00158                           << pNucleon->Get4Momentum() << " # " << (pNucleon->Get4Momentum()).mag()<< G4endl;
00159                                 intNuclMom += nucleonMomentum;
00160                         #endif
00161                                 pNucleon->SetMomentum(nucleonMomentum);
00162                                 G4QGSMSplitableHadron* aTarget = new G4QGSMSplitableHadron(*pNucleon);
00163                                 G4QGSParticipants_NPart ++;
00164                                 theTargets.push_back(aTarget);
00165                                 pNucleon->Hit(aTarget);
00166                                 if ((theProbability.GetDiffractiveProbability(s_nucleon, Distance2)/Probability > G4UniformRand()
00167                                                 &&(ModelMode==SOFT)) || (ModelMode==DIFFRACTIVE ))
00168                                 {
00169                                         // diffractive interaction occurs
00170                                         if(IsSingleDiffractive())
00171                                         {
00172                                                 theSingleDiffExcitation.ExciteParticipants(aProjectile, aTarget);
00173                                         }
00174                                         else
00175                                         {
00176                                                 theDiffExcitaton.ExciteParticipants(aProjectile, aTarget);
00177                                         }
00178                                         G4InteractionContent * aInteraction = new G4InteractionContent(aProjectile);
00179                                         aInteraction->SetTarget(aTarget);
00180                                         theInteractions.push_back(aInteraction);
00181                                         aInteraction->SetNumberOfDiffractiveCollisions(1);
00182                                         totalCuts += 1;
00183                                 }
00184                                 else
00185                                 {
00186                                         // nondiffractive soft interaction occurs
00187                                         // sample nCut+1 (cut Pomerons) pairs of strings can be produced
00188                                         G4int nCut;
00189                                         G4double * running = new G4double[nCutMax];
00190                                         running[0] = 0;
00191                                         for(nCut = 0; nCut < nCutMax; nCut++)
00192                                         {
00193                                                 running[nCut] = theProbability.GetCutPomeronProbability(s_nucleon, Distance2, nCut + 1);
00194                                                 if(nCut!=0) running[nCut] += running[nCut-1];
00195                                         }
00196                                         G4double random = running[nCutMax-1]*G4UniformRand();
00197                                         for(nCut = 0; nCut < nCutMax; nCut++)
00198                                         {
00199                                                 if(running[nCut] > random) break;
00200                                         }
00201                                         delete [] running;
00202                                         nCut = 0;
00203                                         aTarget->IncrementCollisionCount(nCut+1);
00204                                         aProjectile->IncrementCollisionCount(nCut+1);
00205                                         G4InteractionContent * aInteraction = new G4InteractionContent(aProjectile);
00206                                         aInteraction->SetTarget(aTarget);
00207                                         aInteraction->SetNumberOfSoftCollisions(nCut+1);
00208                                         theInteractions.push_back(aInteraction);
00209                                         totalCuts += nCut+1;
00210                                         #ifdef debug_G4QGSParticipants
00211                                                 impactUsed=Distance2;
00212                                         #endif
00213                                 }
00214                         }
00215                 }
00216         #ifdef debug_G4QGSParticipants
00217                         G4cout << G4endl<<"NUCLEONCOUNT "<<nucleonCount<<G4endl;
00218                         G4cout << " Interact 4-Vect " << intNuclMom << G4endl;
00219         #endif
00220         }
00221         #ifdef debug_G4QGSParticipants
00222                 G4cout << G4endl<<"CUTDEBUG "<< totalCuts <<G4endl;
00223                 G4cout << "Impact Parameter used = "<<impactUsed<<G4endl;
00224         #endif
00225         return aProjectile;
00226 }
00227 
00228 void G4QGSParticipants::PerformDiffractiveCollisions()
00229 {
00230         // remove the "G4PartonPair::PROJECTILE", etc., which are not necessary. @@@
00231         unsigned int i;
00232         for(i = 0; i < theInteractions.size(); i++)
00233         {
00234                 G4InteractionContent* anIniteraction = theInteractions[i];
00235                 G4VSplitableHadron* aProjectile = anIniteraction->GetProjectile();
00236                 G4Parton* aParton = aProjectile->GetNextParton();
00237                 G4PartonPair * aPartonPair;
00238                 // projectile first HPW
00239                 if (aParton)
00240                 {
00241                         aPartonPair = new G4PartonPair(aParton, aProjectile->GetNextAntiParton(),
00242                                         G4PartonPair::DIFFRACTIVE,
00243                                         G4PartonPair::PROJECTILE);
00244                         #ifdef debug_G4QGSPart_PDiffColl
00245                         G4cout << "DiffPair Pro " << aPartonPair->GetParton1()->GetPDGcode() << " "
00246                                         << aPartonPair->GetParton1()->Get4Momentum() << " "
00247                                         << aPartonPair->GetParton1()->GetX() << " " << G4endl;
00248                         G4cout << "         " << aPartonPair->GetParton2()->GetPDGcode() << " "
00249                                         << aPartonPair->GetParton2()->Get4Momentum() << " "
00250                                         << aPartonPair->GetParton2()->GetX() << " " << G4endl;
00251                         #endif
00252                         thePartonPairs.push_back(aPartonPair);
00253                 }
00254                 // then target HPW
00255                 G4VSplitableHadron* aTarget = anIniteraction->GetTarget();
00256                 aParton = aTarget->GetNextParton();
00257                 if (aParton)
00258                 {
00259                         aPartonPair = new G4PartonPair(aParton, aTarget->GetNextAntiParton(),
00260                                         G4PartonPair::DIFFRACTIVE,
00261                                         G4PartonPair::TARGET);
00262                                         #ifdef debug_G4QGSPart_PDiffColl
00263                                         G4cout << "DiffPair Tgt " << aPartonPair->GetParton1()->GetPDGcode() << " "
00264                                                         << aPartonPair->GetParton1()->Get4Momentum() << " "
00265                                                         << aPartonPair->GetParton1()->GetX() << " " << G4endl;
00266                                         G4cout << "         " << aPartonPair->GetParton2()->GetPDGcode() << " "
00267                                                         << aPartonPair->GetParton2()->Get4Momentum() << " "
00268                                                         << aPartonPair->GetParton2()->GetX() << " " << G4endl;
00269                                         #endif
00270                                         thePartonPairs.push_back(aPartonPair);
00271                 }
00272         }
00273 }
00274 
00275 void G4QGSParticipants::PerformSoftCollisions()
00276 {
00277         std::vector<G4InteractionContent*>::iterator i;
00278         G4LorentzVector str4Mom;
00279         i = theInteractions.begin();
00280         while ( i != theInteractions.end() )
00281         {
00282                 G4InteractionContent* anIniteraction = *i;
00283                 G4PartonPair * aPair = NULL;
00284                 if (anIniteraction->GetNumberOfSoftCollisions())
00285                 {
00286                         G4VSplitableHadron* pProjectile = anIniteraction->GetProjectile();
00287                         G4VSplitableHadron* pTarget     = anIniteraction->GetTarget();
00288                         for (G4int j = 0; j < anIniteraction->GetNumberOfSoftCollisions(); j++)
00289                         {
00290                                 aPair = new G4PartonPair(pTarget->GetNextParton(), pProjectile->GetNextAntiParton(),
00291                                                 G4PartonPair::SOFT, G4PartonPair::TARGET);
00292                                 #ifdef debug_G4QGSPart_PSoftColl
00293                                 G4cout << "SoftPair " << aPair->GetParton1()->GetPDGcode() << " "
00294                                                 << aPair->GetParton1()->Get4Momentum() << " "
00295                                                 << aPair->GetParton1()->GetX() << " " << G4endl;
00296                                 G4cout << "         " << aPair->GetParton2()->GetPDGcode() << " "
00297                                                 << aPair->GetParton2()->Get4Momentum() << " "
00298                                                 << aPair->GetParton2()->GetX() << " " << G4endl;
00299                                 #endif
00300                                 #ifdef debug_G4QGSParticipants
00301                                 str4Mom += aPair->GetParton1()->Get4Momentum();
00302                                 str4Mom += aPair->GetParton2()->Get4Momentum();
00303                                 #endif
00304                                 thePartonPairs.push_back(aPair);
00305                                 aPair = new G4PartonPair(pProjectile->GetNextParton(), pTarget->GetNextAntiParton(),
00306                                                 G4PartonPair::SOFT, G4PartonPair::PROJECTILE);
00307                                 #ifdef debug_G4QGSPart_PSoftColl
00308                                 G4cout << "SoftPair " << aPair->GetParton1()->GetPDGcode() << " "
00309                                                 << aPair->GetParton1()->Get4Momentum() << " "
00310                                                 << aPair->GetParton1()->GetX() << " " << G4endl;
00311                                 G4cout << "         " << aPair->GetParton2()->GetPDGcode() << " "
00312                                                 << aPair->GetParton2()->Get4Momentum() << " "
00313                                                 << aPair->GetParton2()->GetX() << " " << G4endl;
00314                                 #endif
00315                                 #ifdef debug_G4QGSParticipants
00316                                 str4Mom += aPair->GetParton1()->Get4Momentum();
00317                                 str4Mom += aPair->GetParton2()->Get4Momentum();
00318                                 #endif
00319                                 thePartonPairs.push_back(aPair);
00320                         }
00321                         delete *i;
00322                         i=theInteractions.erase(i);    // i now points to the next interaction
00323                 } else {
00324                         i++;
00325                 }
00326         }
00327         #ifdef debug_G4QGSPart_PSoftColl
00328                 G4cout << " string 4 mom " << str4Mom << G4endl;
00329         #endif
00330 }

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