source: trunk/source/processes/hadronic/models/parton_string/qgsm/src/G4QGSMSplitableHadron.cc @ 1200

Last change on this file since 1200 was 819, checked in by garnier, 16 years ago

import all except CVS
File size: 16.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 "G4QGSMSplitableHadron.hh"
27#include "G4ParticleTable.hh"
28#include "G4PionPlus.hh"
29#include "G4PionMinus.hh"
30#include "G4Gamma.hh"
31#include "G4PionZero.hh"
32#include "G4KaonPlus.hh"
33#include "G4KaonMinus.hh"
34
35// based on prototype by Maxim Komogorov
36// Splitting into methods, and centralizing of model parameters HPW Feb 1999
37// restructuring HPW Feb 1999
38// fixing bug in the sampling of 'x', HPW Feb 1999
39// fixing bug in sampling pz, HPW Feb 1999.
40// Code now also good for p-nucleus scattering (before only p-p), HPW Feb 1999.
41// Using Parton more directly, HPW Feb 1999.
42// Shortening the algorithm for sampling x, HPW Feb 1999.
43// sampling of x replaced by formula, taking X_min into account in the correlated sampling. HPW, Feb 1999.
44// logic much clearer now. HPW Feb 1999
45// Removed the ordering problem. No Direction needed in selection of valence quark types. HPW Mar'99.
46// Fixing p-t distributions for scattering of nuclei.
47// Separating out parameters.
48
49void G4QGSMSplitableHadron::InitParameters()
50{
51  // changing rapidity distribution for all
52  alpha = -0.5; // Note that this number is still assumed in the algorithm
53                // needs to be generalized.
54  // changing rapidity distribution for projectile like
55  beta = 2.5;// Note that this number is still assumed in the algorithm
56                // needs to be generalized.
57  theMinPz = 0.5*G4PionMinus::PionMinus()->GetPDGMass(); 
58//  theMinPz = 0.1*G4PionMinus::PionMinus()->GetPDGMass();
59//  theMinPz = G4PionMinus::PionMinus()->GetPDGMass();
60  // as low as possible, otherwise, we have unphysical boundary conditions in the sampling.
61  StrangeSuppress = 0.48;
62  sigmaPt = 0.*GeV; // widens eta slightly, if increased to 1.7,
63                    // but Maxim's algorithm breaks energy conservation
64                    // to be revised.
65  widthOfPtSquare = 0.01*GeV*GeV;
66  Direction = FALSE;
67  minTransverseMass = 1*keV;
68} 
69
70G4QGSMSplitableHadron::G4QGSMSplitableHadron() 
71{
72  InitParameters();
73}
74
75G4QGSMSplitableHadron::G4QGSMSplitableHadron(const G4ReactionProduct & aPrimary, G4bool aDirection)
76 :G4VSplitableHadron(aPrimary)
77{
78  InitParameters();
79  Direction = aDirection;
80}
81
82 
83G4QGSMSplitableHadron::G4QGSMSplitableHadron(const G4ReactionProduct & aPrimary)
84      :  G4VSplitableHadron(aPrimary)
85{
86  InitParameters();
87}
88
89G4QGSMSplitableHadron::G4QGSMSplitableHadron(const G4Nucleon & aNucleon)
90      :  G4VSplitableHadron(aNucleon)
91{
92  InitParameters();
93}
94
95G4QGSMSplitableHadron::G4QGSMSplitableHadron(const G4Nucleon & aNucleon, G4bool aDirection)
96      :  G4VSplitableHadron(aNucleon)
97{
98  InitParameters();
99  Direction = aDirection;
100}
101
102G4QGSMSplitableHadron::~G4QGSMSplitableHadron(){}
103
104const G4QGSMSplitableHadron & G4QGSMSplitableHadron::operator=(const G4QGSMSplitableHadron &)
105{
106  throw G4HadronicException(__FILE__, __LINE__, "G4QGSMSplitableHadron::operator= meant to not be accessable");
107  return *this;
108}
109
110
111//**************************************************************************************************************************
112   
113void G4QGSMSplitableHadron::SplitUp()
114{ 
115  if (IsSplit()) return;
116  Splitting();
117  if (Color.size()!=0) return;
118  if (GetSoftCollisionCount() == 0)
119  {
120    DiffractiveSplitUp();
121  }
122  else
123  {
124    SoftSplitUp();
125  }
126}
127   
128void G4QGSMSplitableHadron::DiffractiveSplitUp()
129{   
130  // take the particle definitions and get the partons HPW
131  G4Parton * Left = NULL;
132  G4Parton * Right = NULL;
133  GetValenceQuarkFlavors(GetDefinition(), Left, Right);
134  Left->SetPosition(GetPosition());
135  Right->SetPosition(GetPosition());
136 
137  G4LorentzVector HadronMom = Get4Momentum();
138  //std::cout << "DSU 1 - "<<HadronMom<<std::endl;
139
140  // momenta of string ends
141  G4double pt2 = HadronMom.perp2();
142  G4double transverseMass2 = HadronMom.plus()*HadronMom.minus();
143  G4double maxAvailMomentum2 = sqr(std::sqrt(transverseMass2) - std::sqrt(pt2));
144  G4ThreeVector pt(minTransverseMass, minTransverseMass, 0);
145  if(maxAvailMomentum2/widthOfPtSquare>0.01) pt = GaussianPt(widthOfPtSquare, maxAvailMomentum2);
146  //std::cout << "DSU 1.1 - "<< maxAvailMomentum2<< pt <<std::endl;
147
148  G4LorentzVector LeftMom(pt, 0.);
149  G4LorentzVector RightMom;
150  RightMom.setPx(HadronMom.px() - pt.x());
151  RightMom.setPy(HadronMom.py() - pt.y());
152  //std::cout << "DSU 2 - "<<RightMom<<" "<< LeftMom <<std::endl;
153
154  G4double Local1 = HadronMom.minus() + (RightMom.perp2() - LeftMom.perp2())/HadronMom.plus();
155  G4double Local2 = std::sqrt(std::max(0., sqr(Local1) - 4.*RightMom.perp2()*HadronMom.minus()/HadronMom.plus()));
156  //std::cout << "DSU 3 - "<< Local1 <<" "<< Local2 <<std::endl;
157  if (Direction) Local2 = -Local2;
158  G4double RightMinus   = 0.5*(Local1 + Local2);
159  G4double LeftMinus = HadronMom.minus() - RightMinus;
160  //std::cout << "DSU 4 - "<< RightMinus <<" "<< LeftMinus << " "<<HadronMom.minus() <<std::endl;
161
162  G4double LeftPlus  = LeftMom.perp2()/LeftMinus;
163  G4double RightPlus = HadronMom.plus() - LeftPlus;
164  //std::cout << "DSU 5 - "<< RightPlus <<" "<< LeftPlus <<std::endl;
165  LeftMom.setPz(0.5*(LeftPlus - LeftMinus));
166  LeftMom.setE (0.5*(LeftPlus + LeftMinus));
167  RightMom.setPz(0.5*(RightPlus - RightMinus));
168  RightMom.setE (0.5*(RightPlus + RightMinus));
169  //std::cout << "DSU 6 - "<< LeftMom <<" "<< RightMom <<std::endl;
170  Left->Set4Momentum(LeftMom);
171  Right->Set4Momentum(RightMom);
172  Color.push_back(Left);
173  AntiColor.push_back(Right);
174}
175
176
177void G4QGSMSplitableHadron::SoftSplitUp()
178{
179   //... sample transversal momenta for sea and valence quarks
180   G4double phi, pts;
181   G4double SumPy = 0.;
182   G4double SumPx = 0.;
183   G4ThreeVector Pos    = GetPosition();
184   G4int nSeaPair = GetSoftCollisionCount()-1; 
185   
186   // here the condition,to ensure viability of splitting, also in cases
187   // where difractive excitation occured together with soft scattering.
188   //   G4double LightConeMomentum = (Direction)? Get4Momentum().plus() : Get4Momentum().minus();
189   //   G4double Xmin = theMinPz/LightConeMomentum;
190   G4double Xmin = theMinPz/( Get4Momentum().e() - GetDefinition()->GetPDGMass() );
191   while(Xmin>=1-(2*nSeaPair+1)*Xmin) Xmin*=0.95;
192
193   G4int aSeaPair;
194   for (aSeaPair = 0; aSeaPair < nSeaPair; aSeaPair++)
195   {
196     //  choose quark flavour, d:u:s = 1:1:(1/StrangeSuppress-2)
197
198     G4int aPDGCode = 1 + (G4int)(G4UniformRand()/StrangeSuppress); 
199
200     //  BuildSeaQuark() determines quark spin, isospin and colour
201     //  via parton-constructor G4Parton(aPDGCode)
202
203     G4Parton * aParton = BuildSeaQuark(false, aPDGCode, nSeaPair);
204
205//              G4cerr << "G4QGSMSplitableHadron::SoftSplitUp()" << G4endl;
206
207//              G4cerr << "Parton 1: "
208//                     << " PDGcode: "  << aPDGCode
209//                     << " - Name: "   << aParton->GetDefinition()->GetParticleName()
210//                     << " - Type: "   << aParton->GetDefinition()->GetParticleType()
211//                     << " - Spin-3: " << aParton->GetSpinZ()
212//                     << " - Colour: " << aParton->GetColour() << G4endl;
213
214     // save colour a spin-3 for anti-quark
215
216     G4int firstPartonColour = aParton->GetColour();
217     G4double firstPartonSpinZ = aParton->GetSpinZ();
218
219     SumPx += aParton->Get4Momentum().px();
220     SumPy += aParton->Get4Momentum().py();
221     Color.push_back(aParton);
222
223     // create anti-quark
224
225     aParton = BuildSeaQuark(true, aPDGCode, nSeaPair);
226     aParton->SetSpinZ(-firstPartonSpinZ);
227     aParton->SetColour(-firstPartonColour);
228
229//              G4cerr << "Parton 2: "
230//                     << " PDGcode: "  << -aPDGCode
231//                     << " - Name: "   << aParton->GetDefinition()->GetParticleName()
232//                     << " - Type: "   << aParton->GetDefinition()->GetParticleType()
233//                     << " - Spin-3: " << aParton->GetSpinZ()
234//                     << " - Colour: " << aParton->GetColour() << G4endl;
235//              G4cerr << "------------" << G4endl;
236
237     SumPx += aParton->Get4Momentum().px();
238     SumPy += aParton->Get4Momentum().py();
239     AntiColor.push_back(aParton);
240   }
241   // Valence quark   
242   G4Parton* pColorParton = NULL;   
243   G4Parton* pAntiColorParton = NULL;   
244   GetValenceQuarkFlavors(GetDefinition(), pColorParton, pAntiColorParton);
245   G4int ColorEncoding = pColorParton->GetPDGcode();
246   G4int AntiColorEncoding = pAntiColorParton->GetPDGcode();
247   
248   pts   =  sigmaPt*std::sqrt(-std::log(G4UniformRand()));
249   phi   = 2.*pi*G4UniformRand();
250   G4double Px = pts*std::cos(phi);
251   G4double Py = pts*std::sin(phi);
252   SumPx += Px;
253   SumPy += Py;
254
255   if (ColorEncoding < 0) // use particle definition
256   {
257     G4LorentzVector ColorMom(-SumPx, -SumPy, 0, 0);
258     pColorParton->Set4Momentum(ColorMom);
259     G4LorentzVector AntiColorMom(Px, Py, 0, 0);
260     pAntiColorParton->Set4Momentum(AntiColorMom);
261   }
262   else
263   {
264     G4LorentzVector ColorMom(Px, Py, 0, 0);
265     pColorParton->Set4Momentum(ColorMom);
266     G4LorentzVector AntiColorMom(-SumPx, -SumPy, 0, 0);
267     pAntiColorParton->Set4Momentum(AntiColorMom);
268   }
269   Color.push_back(pColorParton);
270   AntiColor.push_back(pAntiColorParton);
271
272   // Sample X
273   G4int nAttempt = 0;
274   G4double SumX = 0;
275   G4double aBeta = beta;
276   G4double ColorX, AntiColorX;
277   G4double HPWtest = 0;
278   if (GetDefinition() == G4PionMinus::PionMinusDefinition()) aBeta = 1.;       
279   if (GetDefinition() == G4Gamma::GammaDefinition()) aBeta = 1.;       
280   if (GetDefinition() == G4PionPlus::PionPlusDefinition()) aBeta = 1.;     
281   if (GetDefinition() == G4PionZero::PionZeroDefinition()) aBeta = 1.;       
282   if (GetDefinition() == G4KaonPlus::KaonPlusDefinition()) aBeta = 0.;       
283   if (GetDefinition() == G4KaonMinus::KaonMinusDefinition()) aBeta = 0.;       
284   do
285   {
286     SumX = 0;
287     nAttempt++;
288     G4int    NumberOfUnsampledSeaQuarks = 2*nSeaPair;
289     G4double beta1 = beta;
290     if (std::abs(ColorEncoding) <= 1000 && std::abs(AntiColorEncoding) <= 1000) beta1 = 1.; //...  in a meson
291     ColorX = SampleX(Xmin, NumberOfUnsampledSeaQuarks, 2*nSeaPair, aBeta);
292     HPWtest = ColorX;
293     while (ColorX < Xmin || ColorX > 1.|| 1. -  ColorX <= Xmin) {;} 
294     Color.back()->SetX(SumX = ColorX);// this is the valenz quark.
295     for(G4int aPair = 0; aPair < nSeaPair; aPair++) 
296     {
297       NumberOfUnsampledSeaQuarks--;
298       ColorX = SampleX(Xmin, NumberOfUnsampledSeaQuarks, 2*nSeaPair, aBeta);
299       Color[aPair]->SetX(ColorX);
300       SumX += ColorX; 
301       NumberOfUnsampledSeaQuarks--;
302       AntiColorX = SampleX(Xmin, NumberOfUnsampledSeaQuarks, 2*nSeaPair, aBeta);
303       AntiColor[aPair]->SetX(AntiColorX); // the 'sea' partons
304       SumX += AntiColorX;
305       if (1. - SumX <= Xmin)  break;
306     }
307   } 
308   while (1. - SumX <= Xmin); 
309   (*(AntiColor.end()-1))->SetX(1. - SumX); // the di-quark takes the rest, then go to momentum
310   /// and here is the bug ;-) @@@@@@@@@@@@@
311   if(getenv("debug_QGSMSplitableHadron") )G4cout << "particle energy at split = "<<Get4Momentum().t()<<G4endl;
312   G4double lightCone = ((!Direction) ? Get4Momentum().minus() : Get4Momentum().plus());
313    G4double lightCone2 = ((!Direction) ? Get4Momentum().plus() : Get4Momentum().minus());
314  // lightCone -= 0.5*Get4Momentum().m();
315   // hpw testing @@@@@ lightCone = 2.*Get4Momentum().t();
316   if(getenv("debug_QGSMSplitableHadron") )G4cout << "Light cone = "<<lightCone<<G4endl;
317   for(aSeaPair = 0; aSeaPair < nSeaPair+1; aSeaPair++) 
318   {
319     G4Parton* aParton = Color[aSeaPair];
320     aParton->DefineMomentumInZ(lightCone, lightCone2, Direction);
321
322     aParton = AntiColor[aSeaPair]; 
323     aParton->DefineMomentumInZ(lightCone, lightCone2, Direction);
324   } 
325//--DEBUG--   cout <<G4endl<<"XSAMPLE "<<HPWtest<<G4endl;
326   return;
327} 
328
329void G4QGSMSplitableHadron::GetValenceQuarkFlavors(const G4ParticleDefinition * aPart, G4Parton *& Parton1, G4Parton *& Parton2)
330{
331   // Note! convention aEnd = q or (qq)bar and bEnd = qbar or qq.
332  G4int aEnd;
333  G4int bEnd;
334  G4int HadronEncoding = aPart->GetPDGEncoding();
335  if (aPart->GetBaryonNumber() == 0) 
336  {
337    theMesonSplitter.SplitMeson(HadronEncoding, &aEnd, &bEnd);
338  }
339  else
340  {
341    theBaryonSplitter.SplitBarion(HadronEncoding, &aEnd, &bEnd);
342  } 
343
344  Parton1 = new G4Parton(aEnd);
345  Parton1->SetPosition(GetPosition());
346
347//      G4cerr << "G4QGSMSplitableHadron::GetValenceQuarkFlavors()" << G4endl;
348//      G4cerr << "Parton 1: "
349//             << " PDGcode: "  << aEnd
350//             << " - Name: "   << Parton1->GetDefinition()->GetParticleName()
351//             << " - Type: "   << Parton1->GetDefinition()->GetParticleType()
352//             << " - Spin-3: " << Parton1->GetSpinZ()
353//             << " - Colour: " << Parton1->GetColour() << G4endl;
354
355  Parton2 = new G4Parton(bEnd);
356  Parton2->SetPosition(GetPosition());
357
358//      G4cerr << "Parton 2: "
359//             << " PDGcode: "  << bEnd
360//             << " - Name: "   << Parton2->GetDefinition()->GetParticleName()
361//             << " - Type: "   << Parton2->GetDefinition()->GetParticleType()
362//             << " - Spin-3: " << Parton2->GetSpinZ()
363//             << " - Colour: " << Parton2->GetColour() << G4endl;
364//      G4cerr << "... now checking for color and spin conservation - yielding: " << G4endl;
365
366  // colour of parton 1 choosen at random by G4Parton(aEnd)
367  // colour of parton 2 is the opposite:
368
369  Parton2->SetColour(-(Parton1->GetColour()));
370
371        // isospin-3 of both partons is handled by G4Parton(PDGCode)
372
373        // spin-3 of parton 1 and 2 choosen at random by G4Parton(aEnd)
374        // spin-3 of parton 2 may be constrained by spin of original particle:
375
376  if ( std::abs(Parton1->GetSpinZ() + Parton2->GetSpinZ()) > aPart->GetPDGSpin()) 
377  {
378                Parton2->SetSpinZ(-(Parton2->GetSpinZ()));   
379  } 
380
381//      G4cerr << "Parton 2: "
382//             << " PDGcode: "  << bEnd
383//             << " - Name: "   << Parton2->GetDefinition()->GetParticleName()
384//             << " - Type: "   << Parton2->GetDefinition()->GetParticleType()
385//             << " - Spin-3: " << Parton2->GetSpinZ()
386//             << " - Colour: " << Parton2->GetColour() << G4endl;
387//      G4cerr << "------------" << G4endl;
388
389}
390 
391   
392G4ThreeVector G4QGSMSplitableHadron::GaussianPt(G4double widthSquare, G4double maxPtSquare)
393{
394  G4double R;
395  while((R = -widthSquare*std::log(G4UniformRand())) > maxPtSquare) {;}
396  R = std::sqrt(R);
397  G4double phi = twopi*G4UniformRand();
398  return G4ThreeVector (R*std::cos(phi), R*std::sin(phi), 0.);   
399}
400
401G4Parton * G4QGSMSplitableHadron::
402BuildSeaQuark(G4bool isAntiQuark, G4int aPDGCode, G4int /* nSeaPair*/)
403{
404  if (isAntiQuark) aPDGCode*=-1;
405  G4Parton* result = new G4Parton(aPDGCode);   
406  result->SetPosition(GetPosition());
407  G4ThreeVector aPtVector = GaussianPt(sigmaPt, DBL_MAX);
408  G4LorentzVector a4Momentum(aPtVector, 0);
409  result->Set4Momentum(a4Momentum);
410  return result;
411}
412
413G4double G4QGSMSplitableHadron::
414SampleX(G4double anXmin, G4int nSea, G4int totalSea, G4double aBeta)
415{
416  G4double result;
417  G4double x1, x2;
418  G4double ymax = 0;
419  for(G4int ii=1; ii<100; ii++)
420  {
421    G4double y = std::pow(1./G4double(ii), alpha);
422    y *= std::pow( std::pow(1-anXmin-totalSea*anXmin, alpha+1) - std::pow(anXmin, alpha+1), nSea);
423    y *= std::pow(1-anXmin-totalSea*anXmin, aBeta+1) - std::pow(anXmin, aBeta+1);
424    if(y>ymax) ymax = y;
425  }
426  G4double y;
427  G4double xMax=1-(totalSea+1)*anXmin;
428  if(anXmin > xMax) 
429  {
430    G4cout << "anXmin = "<<anXmin<<" nSea = "<<nSea<<" totalSea = "<< totalSea<<G4endl;
431    throw G4HadronicException(__FILE__, __LINE__, "G4QGSMSplitableHadron - Fatal: Cannot sample parton densities under these constraints.");
432  }
433  do
434  {
435    x1 = CLHEP::RandFlat::shoot(anXmin, xMax);
436    y = std::pow(x1, alpha);
437    y *= std::pow( std::pow(1-x1-totalSea*anXmin, alpha+1) - std::pow(anXmin, alpha+1), nSea);
438    y *= std::pow(1-x1-totalSea*anXmin, aBeta+1) - std::pow(anXmin, aBeta+1); 
439    x2 = ymax*G4UniformRand();
440  }
441  while(x2>y);
442  result = x1;
443  return result; 
444}
Note: See TracBrowser for help on using the repository browser.