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

Last change on this file since 1199 was 819, checked in by garnier, 17 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.