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

Last change on this file since 1348 was 1055, checked in by garnier, 15 years ago

maj sur la beta de geant 4.9.3
File size: 12.8 KB
Line 
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),ThresholdParameter(0.45*GeV),
39                                         nCutMax(7),ThresholdParameter(0.000*GeV),
40                                         QGSMThreshold(3*GeV),theNucleonRadius(1.5*fermi)
41                                         
42{
43}
44
45G4QGSParticipants::G4QGSParticipants(const G4QGSParticipants &right)
46: G4VParticipants(), nCutMax(right.nCutMax),ThresholdParameter(right.ThresholdParameter),
47  QGSMThreshold(right.QGSMThreshold),theNucleonRadius(right.theNucleonRadius)
48{
49}
50
51
52G4QGSParticipants::~G4QGSParticipants()
53{
54}
55
56void G4QGSParticipants::BuildInteractions(const G4ReactionProduct  &thePrimary) 
57{
58 
59  // Find the collisions and collition conditions
60  G4VSplitableHadron* aProjectile = SelectInteractions(thePrimary);
61
62  // now build the parton pairs. HPW
63  SplitHadrons();
64 
65  // soft collisions first HPW, ordering is vital
66  PerformSoftCollisions();
67   
68  // the rest is diffractive HPW
69  PerformDiffractiveCollisions();
70 
71  // clean-up, if necessary
72  std::for_each(theInteractions.begin(), theInteractions.end(), DeleteInteractionContent());
73  theInteractions.clear();
74  std::for_each(theTargets.begin(), theTargets.end(), DeleteSplitableHadron());
75  theTargets.clear();
76  delete aProjectile;
77}
78
79G4VSplitableHadron* G4QGSParticipants::SelectInteractions(const G4ReactionProduct  &thePrimary)
80{
81  G4VSplitableHadron* aProjectile = new G4QGSMSplitableHadron(thePrimary, TRUE); // @@@ check the TRUE
82  G4PomeronCrossSection theProbability(thePrimary.GetDefinition()); // @@@ should be data member
83  G4double outerRadius = theNucleus->GetOuterRadius();
84  // Check reaction threshold
85
86  theNucleus->StartLoop();
87  G4Nucleon * pNucleon = theNucleus->GetNextNucleon();
88  G4LorentzVector aPrimaryMomentum(thePrimary.GetMomentum(), thePrimary.GetTotalEnergy());
89//--DEBUG--G4cout << " qgspart- " << aPrimaryMomentum << " # " << aPrimaryMomentum.mag()
90//--DEBUG--      << pNucleon->Get4Momentum() << " # " << (pNucleon->Get4Momentum()).mag()<< G4endl;
91  G4double s = (aPrimaryMomentum + pNucleon->Get4Momentum()).mag2();
92  G4double ThresholdMass = thePrimary.GetMass() + pNucleon->GetDefinition()->GetPDGMass(); 
93  ModelMode = SOFT;
94  if (sqr(ThresholdMass + ThresholdParameter) > s)
95  {
96    ModelMode = DIFFRACTIVE;
97    //throw G4HadronicException(__FILE__, __LINE__, "Initial energy is too low. The 4-vectors of the input are inconsistant with the particle masses.");
98  }
99  if (sqr(ThresholdMass + QGSMThreshold) > s) // thus only diffractive in cascade!
100  {
101    ModelMode = DIFFRACTIVE;
102  }
103 
104  // first find the collisions HPW
105  std::for_each(theInteractions.begin(), theInteractions.end(), DeleteInteractionContent());
106  theInteractions.clear();
107  G4int totalCuts = 0;
108  G4double impactUsed = 0;
109
110  #ifdef debug_QGS
111  G4double eK = thePrimary.GetKineticEnergy()/GeV;
112  #endif
113#ifdef debug_G4QGSParticipants
114  G4LorentzVector intNuclMom;
115#endif
116  while(theInteractions.size() == 0)
117  {
118    // choose random impact parameter HPW
119    std::pair<G4double, G4double> theImpactParameter;
120    theImpactParameter = theNucleus->ChooseImpactXandY(outerRadius+theNucleonRadius);
121    G4double impactX = theImpactParameter.first; 
122    G4double impactY = theImpactParameter.second;
123   
124    // loop over nuclei to find collissions HPW
125    theNucleus->StartLoop();
126    G4int nucleonCount = 0; // debug
127    G4QGSParticipants_NPart = 0;
128#ifdef debug_G4QGSParticipants
129    intNuclMom=aPrimaryMomentum;
130#endif
131    while( (pNucleon = theNucleus->GetNextNucleon()) )
132    {
133      if(totalCuts>1.5*thePrimary.GetKineticEnergy()/GeV) 
134      {
135        break;
136      }
137       nucleonCount++; // debug
138      // Needs to be moved to Probability class @@@
139      G4double s = (aPrimaryMomentum + pNucleon->Get4Momentum()).mag2();
140      G4double Distance2 = sqr(impactX - pNucleon->GetPosition().x()) +
141                           sqr(impactY - pNucleon->GetPosition().y());
142      G4double Probability = theProbability.GetInelasticProbability(s, Distance2); 
143      // test for inelastic collision
144      G4double rndNumber = G4UniformRand();
145//      ModelMode = DIFFRACTIVE;
146      if (Probability > rndNumber)
147      {
148//--DEBUG--        cout << "DEBUG p="<< Probability<<" r="<<rndNumber<<" d="<<std::sqrt(Distance2)<<G4endl;
149//--DEBUG--         G4cout << " qgspart+  " << aPrimaryMomentum << " # " << aPrimaryMomentum.mag()
150//--DEBUG--              << pNucleon->Get4Momentum() << " # " << (pNucleon->Get4Momentum()).mag()<< G4endl;
151
152#ifdef debug_G4QGSParticipants
153       intNuclMom += pNucleon->Get4Momentum();
154#endif
155       G4QGSMSplitableHadron* aTarget = new G4QGSMSplitableHadron(*pNucleon);
156        G4QGSParticipants_NPart ++;
157        theTargets.push_back(aTarget);
158        pNucleon->Hit(aTarget);
159        if ((theProbability.GetDiffractiveProbability(s, Distance2)/Probability > G4UniformRand() 
160             &&(ModelMode==SOFT)) || (ModelMode==DIFFRACTIVE ))
161        { 
162          // diffractive interaction occurs
163          if(IsSingleDiffractive())
164          {
165            theSingleDiffExcitation.ExciteParticipants(aProjectile, aTarget);
166          }
167          else
168          {
169            theDiffExcitaton.ExciteParticipants(aProjectile, aTarget);
170          }
171          G4InteractionContent * aInteraction = new G4InteractionContent(aProjectile);
172          aInteraction->SetTarget(aTarget); 
173          theInteractions.push_back(aInteraction);
174          aInteraction->SetNumberOfDiffractiveCollisions(1);
175          totalCuts += 1;
176        }
177        else
178        {
179          // nondiffractive soft interaction occurs
180          // sample nCut+1 (cut Pomerons) pairs of strings can be produced
181          G4int nCut;
182          G4double * running = new G4double[nCutMax];
183          running[0] = 0;
184          for(nCut = 0; nCut < nCutMax; nCut++)
185          {
186            running[nCut] = theProbability.GetCutPomeronProbability(s, Distance2, nCut + 1);
187            if(nCut!=0) running[nCut] += running[nCut-1];
188          }
189          G4double random = running[nCutMax-1]*G4UniformRand();
190          for(nCut = 0; nCut < nCutMax; nCut++)
191          {
192             if(running[nCut] > random) break; 
193          }
194          delete [] running;
195            nCut = 0;
196          aTarget->IncrementCollisionCount(nCut+1);
197          aProjectile->IncrementCollisionCount(nCut+1);
198          G4InteractionContent * aInteraction = new G4InteractionContent(aProjectile);
199          aInteraction->SetTarget(aTarget);
200          aInteraction->SetNumberOfSoftCollisions(nCut+1);
201          theInteractions.push_back(aInteraction);
202          totalCuts += nCut+1;
203          impactUsed=Distance2;
204       }
205      }
206    }
207//--DEBUG--  g4cout << G4endl<<"NUCLEONCOUNT "<<nucleonCount<<G4endl;
208//--DEBUG--  G4cout << " Interact 4-Vect " << intNuclMom << G4endl;
209  }
210//--DEBUG--  cout << G4endl<<"CUTDEBUG "<< totalCuts <<G4endl;
211//--DEBUG--  cout << "Impact Parameter used = "<<impactUsed<<G4endl;
212  return aProjectile;
213}
214
215void G4QGSParticipants::PerformDiffractiveCollisions()
216{
217  // remove the "G4PartonPair::PROJECTILE", etc., which are not necessary. @@@
218  unsigned int i;
219  for(i = 0; i < theInteractions.size(); i++) 
220  {
221    G4InteractionContent* anIniteraction = theInteractions[i];
222    G4VSplitableHadron* aProjectile = anIniteraction->GetProjectile();
223    G4Parton* aParton = aProjectile->GetNextParton();
224    G4PartonPair * aPartonPair;
225    // projectile first HPW
226    if (aParton)
227    {
228      aPartonPair = new G4PartonPair(aParton, aProjectile->GetNextAntiParton(), 
229                                     G4PartonPair::DIFFRACTIVE, 
230                                     G4PartonPair::PROJECTILE);
231#ifdef debug_G4QGSPart_PDiffColl
232        G4cout << "DiffPair Pro " << aPartonPair->GetParton1()->GetPDGcode() << " " 
233                              << aPartonPair->GetParton1()->Get4Momentum() << " "
234                              << aPartonPair->GetParton1()->GetX() << " " << G4endl;
235        G4cout << "         " << aPartonPair->GetParton2()->GetPDGcode() << " " 
236                              << aPartonPair->GetParton2()->Get4Momentum() << " "
237                              << aPartonPair->GetParton2()->GetX() << " " << G4endl;
238#endif
239      thePartonPairs.push_back(aPartonPair);
240    }
241    // then target HPW
242    G4VSplitableHadron* aTarget = anIniteraction->GetTarget();
243    aParton = aTarget->GetNextParton();
244    if (aParton)
245    {
246      aPartonPair = new G4PartonPair(aParton, aTarget->GetNextAntiParton(), 
247                                     G4PartonPair::DIFFRACTIVE, 
248                                     G4PartonPair::TARGET);
249#ifdef debug_G4QGSPart_PDiffColl
250        G4cout << "DiffPair Tgt " << aPartonPair->GetParton1()->GetPDGcode() << " " 
251                              << aPartonPair->GetParton1()->Get4Momentum() << " "
252                              << aPartonPair->GetParton1()->GetX() << " " << G4endl;
253        G4cout << "         " << aPartonPair->GetParton2()->GetPDGcode() << " " 
254                              << aPartonPair->GetParton2()->Get4Momentum() << " "
255                              << aPartonPair->GetParton2()->GetX() << " " << G4endl;
256#endif
257      thePartonPairs.push_back(aPartonPair);
258    }
259  }
260}
261
262void G4QGSParticipants::PerformSoftCollisions()
263{
264  std::vector<G4InteractionContent*>::iterator i;
265  G4LorentzVector str4Mom;
266  i = theInteractions.begin();
267  while ( i != theInteractions.end() )   
268  {
269    G4InteractionContent* anIniteraction = *i;
270    G4PartonPair * aPair = NULL;
271    if (anIniteraction->GetNumberOfSoftCollisions())
272    { 
273      G4VSplitableHadron* pProjectile = anIniteraction->GetProjectile();
274      G4VSplitableHadron* pTarget     = anIniteraction->GetTarget();
275      for (G4int j = 0; j < anIniteraction->GetNumberOfSoftCollisions(); j++)
276      {
277        aPair = new G4PartonPair(pTarget->GetNextParton(), pProjectile->GetNextAntiParton(), 
278                                 G4PartonPair::SOFT, G4PartonPair::TARGET);
279#ifdef debug_G4QGSPart_PSoftColl
280        G4cout << "SoftPair " << aPair->GetParton1()->GetPDGcode() << " " 
281                              << aPair->GetParton1()->Get4Momentum() << " "
282                              << aPair->GetParton1()->GetX() << " " << G4endl;
283        G4cout << "         " << aPair->GetParton2()->GetPDGcode() << " " 
284                              << aPair->GetParton2()->Get4Momentum() << " "
285                              << aPair->GetParton2()->GetX() << " " << G4endl;
286#endif
287#ifdef debug_G4QGSParticipants
288        str4Mom += aPair->GetParton1()->Get4Momentum();
289        str4Mom += aPair->GetParton2()->Get4Momentum();
290#endif
291        thePartonPairs.push_back(aPair);
292        aPair = new G4PartonPair(pProjectile->GetNextParton(), pTarget->GetNextAntiParton(), 
293                                 G4PartonPair::SOFT, G4PartonPair::PROJECTILE);
294#ifdef debug_G4QGSPart_PSoftColl
295        G4cout << "SoftPair " << aPair->GetParton1()->GetPDGcode() << " " 
296                              << aPair->GetParton1()->Get4Momentum() << " "
297                              << aPair->GetParton1()->GetX() << " " << G4endl;
298        G4cout << "         " << aPair->GetParton2()->GetPDGcode() << " " 
299                              << aPair->GetParton2()->Get4Momentum() << " "
300                              << aPair->GetParton2()->GetX() << " " << G4endl;
301#endif
302#ifdef debug_G4QGSParticipants
303        str4Mom += aPair->GetParton1()->Get4Momentum();
304        str4Mom += aPair->GetParton2()->Get4Momentum();
305#endif
306        thePartonPairs.push_back(aPair);
307      } 
308      delete *i;
309      i=theInteractions.erase(i);    // i now points to the next interaction
310    } else {   
311      i++;
312    } 
313  }
314#ifdef debug_G4QGSPart_PSoftColl
315  G4cout << " string 4 mom " << str4Mom << G4endl;
316#endif
317}
Note: See TracBrowser for help on using the repository browser.