source: trunk/source/processes/hadronic/models/parton_string/qgsm/src/G4QGSParticipants.cc @ 978

Last change on this file since 978 was 962, checked in by garnier, 15 years ago

update processes
File size: 12.7 KB
RevLine 
[819]1//
2// ********************************************************************
3// * License and Disclaimer                                           *
4// *                                                                  *
5// * The  Geant4 software  is  copyright of the Copyright Holders  of *
6// * the Geant4 Collaboration.  It is provided  under  the terms  and *
7// * conditions of the Geant4 Software License,  included in the file *
8// * LICENSE and available at  http://cern.ch/geant4/license .  These *
9// * include a list of copyright holders.                             *
10// *                                                                  *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work  make  any representation or  warranty, express or implied, *
14// * regarding  this  software system or assume any liability for its *
15// * use.  Please see the license in the file  LICENSE  and URL above *
16// * for the full disclaimer and the limitation of liability.         *
17// *                                                                  *
18// * This  code  implementation is the result of  the  scientific and *
19// * technical work of the GEANT4 collaboration.                      *
20// * By using,  copying,  modifying or  distributing the software (or *
21// * any work based  on the software)  you  agree  to acknowledge its *
22// * use  in  resulting  scientific  publications,  and indicate your *
23// * acceptance of all terms of the Geant4 Software license.          *
24// ********************************************************************
25//
26#include "globals.hh"
27#include "G4QGSParticipants.hh"
28#include "G4LorentzVector.hh"
29#include <utility>
30
31// Class G4QGSParticipants
32
33// HPW Feb 1999
34// Promoting model parameters from local variables class properties
35G4int G4QGSParticipants_NPart = 0;
36
37G4QGSParticipants::G4QGSParticipants() : theDiffExcitaton(), //0.7*GeV, 250*MeV, 250*MeV),
38                                         nCutMax(7),ThersholdParameter(0.45*GeV),
39                                         QGSMThershold(3*GeV),theNucleonRadius(1.5*fermi)
40                                         
41{
42}
43
44G4QGSParticipants::G4QGSParticipants(const G4QGSParticipants &right)
45: G4VParticipants(), nCutMax(right.nCutMax),ThersholdParameter(right.ThersholdParameter),
46  QGSMThershold(right.QGSMThershold),theNucleonRadius(right.theNucleonRadius)
47{
48}
49
50
51G4QGSParticipants::~G4QGSParticipants()
52{
53}
54
55void G4QGSParticipants::BuildInteractions(const G4ReactionProduct  &thePrimary) 
56{
57 
58  // Find the collisions and collition conditions
59  G4VSplitableHadron* aProjectile = SelectInteractions(thePrimary);
60
61  // now build the parton pairs. HPW
62  SplitHadrons();
63 
64  // soft collisions first HPW, ordering is vital
65  PerformSoftCollisions();
66   
67  // the rest is diffractive HPW
68  PerformDiffractiveCollisions();
69 
70  // clean-up, if necessary
71  std::for_each(theInteractions.begin(), theInteractions.end(), DeleteInteractionContent());
72  theInteractions.clear();
73  std::for_each(theTargets.begin(), theTargets.end(), DeleteSplitableHadron());
74  theTargets.clear();
75  delete aProjectile;
76}
77
78G4VSplitableHadron* G4QGSParticipants::SelectInteractions(const G4ReactionProduct  &thePrimary)
79{
80  G4VSplitableHadron* aProjectile = new G4QGSMSplitableHadron(thePrimary, TRUE); // @@@ check the TRUE
81  G4PomeronCrossSection theProbability(thePrimary.GetDefinition()); // @@@ should be data member
82  G4double outerRadius = theNucleus->GetOuterRadius();
83  // Check reaction threshold
84
85  theNucleus->StartLoop();
86  G4Nucleon * pNucleon = theNucleus->GetNextNucleon();
87  G4LorentzVector aPrimaryMomentum(thePrimary.GetMomentum(), thePrimary.GetTotalEnergy());
88//--DEBUG--G4cout << " qgspart- " << aPrimaryMomentum << " # " << aPrimaryMomentum.mag()
89//--DEBUG--      << pNucleon->Get4Momentum() << " # " << (pNucleon->Get4Momentum()).mag()<< G4endl;
90  G4double s = (aPrimaryMomentum + pNucleon->Get4Momentum()).mag2();
91  G4double ThresholdMass = thePrimary.GetMass() + pNucleon->GetDefinition()->GetPDGMass(); 
92  ModelMode = SOFT;
93  if (sqr(ThresholdMass + ThersholdParameter) > s)
94  {
95    throw G4HadronicException(__FILE__, __LINE__, "Initial energy is too low. The 4-vectors of the input are inconsistant with the particle masses.");
96  }
97  if (sqr(ThresholdMass + QGSMThershold) > s) // thus only diffractive in cascade!
98  {
99    ModelMode = DIFFRACTIVE;
100  }
101 
102  // first find the collisions HPW
103  std::for_each(theInteractions.begin(), theInteractions.end(), DeleteInteractionContent());
104  theInteractions.clear();
105  G4int totalCuts = 0;
106  G4double impactUsed = 0;
107
108  #ifdef debug_QGS
109  G4double eK = thePrimary.GetKineticEnergy()/GeV;
110  #endif
111#ifdef debug_G4QGSParticipants
112  G4LorentzVector intNuclMom;
113#endif
114  while(theInteractions.size() == 0)
115  {
116    // choose random impact parameter HPW
117    std::pair<G4double, G4double> theImpactParameter;
118    theImpactParameter = theNucleus->ChooseImpactXandY(outerRadius+theNucleonRadius);
119    G4double impactX = theImpactParameter.first; 
120    G4double impactY = theImpactParameter.second;
121   
122    // loop over nuclei to find collissions HPW
123    theNucleus->StartLoop();
124    G4int nucleonCount = 0; // debug
125    G4QGSParticipants_NPart = 0;
126#ifdef debug_G4QGSParticipants
127    intNuclMom=aPrimaryMomentum;
128#endif
129    while( (pNucleon = theNucleus->GetNextNucleon()) )
130    {
131      if(totalCuts>1.5*thePrimary.GetKineticEnergy()/GeV) 
132      {
133        break;
134      }
135       nucleonCount++; // debug
136      // Needs to be moved to Probability class @@@
137      G4double s = (aPrimaryMomentum + pNucleon->Get4Momentum()).mag2();
138      G4double Distance2 = sqr(impactX - pNucleon->GetPosition().x()) +
139                           sqr(impactY - pNucleon->GetPosition().y());
140      G4double Probability = theProbability.GetInelasticProbability(s, Distance2); 
141      // test for inelastic collision
142      G4double rndNumber = G4UniformRand();
143//      ModelMode = DIFFRACTIVE;
144      if (Probability > rndNumber)
145      {
146//--DEBUG--        cout << "DEBUG p="<< Probability<<" r="<<rndNumber<<" d="<<std::sqrt(Distance2)<<G4endl;
147//--DEBUG--         G4cout << " qgspart+  " << aPrimaryMomentum << " # " << aPrimaryMomentum.mag()
148//--DEBUG--              << pNucleon->Get4Momentum() << " # " << (pNucleon->Get4Momentum()).mag()<< G4endl;
149
150#ifdef debug_G4QGSParticipants
151       intNuclMom += pNucleon->Get4Momentum();
152#endif
153       G4QGSMSplitableHadron* aTarget = new G4QGSMSplitableHadron(*pNucleon);
154        G4QGSParticipants_NPart ++;
155        theTargets.push_back(aTarget);
156        pNucleon->Hit(aTarget);
157        if ((theProbability.GetDiffractiveProbability(s, Distance2)/Probability > G4UniformRand() 
158             &&(ModelMode==SOFT)) || (ModelMode==DIFFRACTIVE ))
159        { 
160          // diffractive interaction occurs
161          if(IsSingleDiffractive())
162          {
163            theSingleDiffExcitation.ExciteParticipants(aProjectile, aTarget);
164          }
165          else
166          {
167            theDiffExcitaton.ExciteParticipants(aProjectile, aTarget);
168          }
169          G4InteractionContent * aInteraction = new G4InteractionContent(aProjectile);
170          aInteraction->SetTarget(aTarget); 
171          theInteractions.push_back(aInteraction);
172          aInteraction->SetNumberOfDiffractiveCollisions(1);
173          totalCuts += 1;
174        }
175        else
176        {
177          // nondiffractive soft interaction occurs
178          // sample nCut+1 (cut Pomerons) pairs of strings can be produced
179          G4int nCut;
180          G4double * running = new G4double[nCutMax];
181          running[0] = 0;
182          for(nCut = 0; nCut < nCutMax; nCut++)
183          {
184            running[nCut] = theProbability.GetCutPomeronProbability(s, Distance2, nCut + 1);
185            if(nCut!=0) running[nCut] += running[nCut-1];
186          }
187          G4double random = running[nCutMax-1]*G4UniformRand();
188          for(nCut = 0; nCut < nCutMax; nCut++)
189          {
190             if(running[nCut] > random) break; 
191          }
192          delete [] running;
193            nCut = 0;
194          aTarget->IncrementCollisionCount(nCut+1);
195          aProjectile->IncrementCollisionCount(nCut+1);
196          G4InteractionContent * aInteraction = new G4InteractionContent(aProjectile);
197          aInteraction->SetTarget(aTarget);
198          aInteraction->SetNumberOfSoftCollisions(nCut+1);
199          theInteractions.push_back(aInteraction);
200          totalCuts += nCut+1;
201          impactUsed=Distance2;
202       }
203      }
204    }
205//--DEBUG--  g4cout << G4endl<<"NUCLEONCOUNT "<<nucleonCount<<G4endl;
206//--DEBUG--  G4cout << " Interact 4-Vect " << intNuclMom << G4endl;
207  }
208//--DEBUG--  cout << G4endl<<"CUTDEBUG "<< totalCuts <<G4endl;
209//--DEBUG--  cout << "Impact Parameter used = "<<impactUsed<<G4endl;
210  return aProjectile;
211}
212
213void G4QGSParticipants::PerformDiffractiveCollisions()
214{
215  // remove the "G4PartonPair::PROJECTILE", etc., which are not necessary. @@@
216  unsigned int i;
217  for(i = 0; i < theInteractions.size(); i++) 
218  {
219    G4InteractionContent* anIniteraction = theInteractions[i];
220    G4VSplitableHadron* aProjectile = anIniteraction->GetProjectile();
221    G4Parton* aParton = aProjectile->GetNextParton();
222    G4PartonPair * aPartonPair;
223    // projectile first HPW
224    if (aParton)
225    {
226      aPartonPair = new G4PartonPair(aParton, aProjectile->GetNextAntiParton(), 
227                                     G4PartonPair::DIFFRACTIVE, 
228                                     G4PartonPair::PROJECTILE);
229#ifdef debug_G4QGSPart_PDiffColl
230        G4cout << "DiffPair Pro " << aPartonPair->GetParton1()->GetPDGcode() << " " 
231                              << aPartonPair->GetParton1()->Get4Momentum() << " "
232                              << aPartonPair->GetParton1()->GetX() << " " << G4endl;
233        G4cout << "         " << aPartonPair->GetParton2()->GetPDGcode() << " " 
234                              << aPartonPair->GetParton2()->Get4Momentum() << " "
235                              << aPartonPair->GetParton2()->GetX() << " " << G4endl;
236#endif
237      thePartonPairs.push_back(aPartonPair);
238    }
239    // then target HPW
240    G4VSplitableHadron* aTarget = anIniteraction->GetTarget();
241    aParton = aTarget->GetNextParton();
242    if (aParton)
243    {
244      aPartonPair = new G4PartonPair(aParton, aTarget->GetNextAntiParton(), 
245                                     G4PartonPair::DIFFRACTIVE, 
246                                     G4PartonPair::TARGET);
247#ifdef debug_G4QGSPart_PDiffColl
248        G4cout << "DiffPair Tgt " << aPartonPair->GetParton1()->GetPDGcode() << " " 
249                              << aPartonPair->GetParton1()->Get4Momentum() << " "
250                              << aPartonPair->GetParton1()->GetX() << " " << G4endl;
251        G4cout << "         " << aPartonPair->GetParton2()->GetPDGcode() << " " 
252                              << aPartonPair->GetParton2()->Get4Momentum() << " "
253                              << aPartonPair->GetParton2()->GetX() << " " << G4endl;
254#endif
255      thePartonPairs.push_back(aPartonPair);
256    }
257  }
258}
259
260void G4QGSParticipants::PerformSoftCollisions()
261{
262  std::vector<G4InteractionContent*>::iterator i;
263  G4LorentzVector str4Mom;
[962]264  i = theInteractions.begin();
265  while ( i != theInteractions.end() )   
[819]266  {
267    G4InteractionContent* anIniteraction = *i;
268    G4PartonPair * aPair = NULL;
269    if (anIniteraction->GetNumberOfSoftCollisions())
270    { 
271      G4VSplitableHadron* pProjectile = anIniteraction->GetProjectile();
272      G4VSplitableHadron* pTarget     = anIniteraction->GetTarget();
273      for (G4int j = 0; j < anIniteraction->GetNumberOfSoftCollisions(); j++)
274      {
275        aPair = new G4PartonPair(pTarget->GetNextParton(), pProjectile->GetNextAntiParton(), 
276                                 G4PartonPair::SOFT, G4PartonPair::TARGET);
277#ifdef debug_G4QGSPart_PSoftColl
278        G4cout << "SoftPair " << aPair->GetParton1()->GetPDGcode() << " " 
279                              << aPair->GetParton1()->Get4Momentum() << " "
280                              << aPair->GetParton1()->GetX() << " " << G4endl;
281        G4cout << "         " << aPair->GetParton2()->GetPDGcode() << " " 
282                              << aPair->GetParton2()->Get4Momentum() << " "
283                              << aPair->GetParton2()->GetX() << " " << G4endl;
284#endif
285#ifdef debug_G4QGSParticipants
286        str4Mom += aPair->GetParton1()->Get4Momentum();
287        str4Mom += aPair->GetParton2()->Get4Momentum();
288#endif
289        thePartonPairs.push_back(aPair);
290        aPair = new G4PartonPair(pProjectile->GetNextParton(), pTarget->GetNextAntiParton(), 
291                                 G4PartonPair::SOFT, G4PartonPair::PROJECTILE);
292#ifdef debug_G4QGSPart_PSoftColl
293        G4cout << "SoftPair " << aPair->GetParton1()->GetPDGcode() << " " 
294                              << aPair->GetParton1()->Get4Momentum() << " "
295                              << aPair->GetParton1()->GetX() << " " << G4endl;
296        G4cout << "         " << aPair->GetParton2()->GetPDGcode() << " " 
297                              << aPair->GetParton2()->Get4Momentum() << " "
298                              << aPair->GetParton2()->GetX() << " " << G4endl;
299#endif
300#ifdef debug_G4QGSParticipants
301        str4Mom += aPair->GetParton1()->Get4Momentum();
302        str4Mom += aPair->GetParton2()->Get4Momentum();
303#endif
304        thePartonPairs.push_back(aPair);
305      } 
306      delete *i;
[962]307      i=theInteractions.erase(i);    // i now points to the next interaction
308    } else {   
309      i++;
310    } 
[819]311  }
312#ifdef debug_G4QGSPart_PSoftColl
313  G4cout << " string 4 mom " << str4Mom << G4endl;
314#endif
315}
Note: See TracBrowser for help on using the repository browser.