source: trunk/source/processes/hadronic/models/parton_string/hadronization/src/G4QGSMFragmentation.cc@ 1315

Last change on this file since 1315 was 1228, checked in by garnier, 16 years ago

update geant4.9.3 tag
File size: 13.7 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//
27// $Id: G4QGSMFragmentation.cc,v 1.9 2008/06/23 08:35:55 vuzhinsk Exp $
28// GEANT4 tag $Name: geant4-09-03 $
29//
30// -----------------------------------------------------------------------------
31// GEANT 4 class implementation file
32//
33// History: first implementation, Maxim Komogorov, 10-Jul-1998
34// -----------------------------------------------------------------------------
35#include "G4QGSMFragmentation.hh"
36#include "G4FragmentingString.hh"
37#include "G4DiQuarks.hh"
38#include "G4Quarks.hh"
39
40#include "Randomize.hh"
41#include "G4ios.hh"
42
43// Class G4QGSMFragmentation
44//****************************************************************************************
45
46G4QGSMFragmentation::G4QGSMFragmentation() :
47arho(0.5), aphi(0.), an(-0.5), ala(-0.75), aksi(-1.), alft(0.5)
48 {
49 }
50
51G4QGSMFragmentation::G4QGSMFragmentation(const G4QGSMFragmentation &) : G4VLongitudinalStringDecay(),
52arho(0.5), aphi(0.), an(-0.5), ala(-0.75), aksi(-1.), alft(0.5)
53 {
54 }
55
56G4QGSMFragmentation::~G4QGSMFragmentation()
57 {
58 }
59
60//****************************************************************************************
61
62const G4QGSMFragmentation & G4QGSMFragmentation::operator=(const G4QGSMFragmentation &)
63 {
64 throw G4HadronicException(__FILE__, __LINE__, "G4QGSMFragmentation::operator= meant to not be accessable");
65 return *this;
66 }
67
68int G4QGSMFragmentation::operator==(const G4QGSMFragmentation &right) const
69 {
70 return !memcmp(this, &right, sizeof(G4QGSMFragmentation));
71 }
72
73int G4QGSMFragmentation::operator!=(const G4QGSMFragmentation &right) const
74 {
75 return memcmp(this, &right, sizeof(G4QGSMFragmentation));
76 }
77
78//****************************************************************************************
79//----------------------------------------------------------------------------------------------------------
80
81G4KineticTrackVector* G4QGSMFragmentation::FragmentString(const G4ExcitedString& theString)
82{
83// Can no longer modify Parameters for Fragmentation.
84 PastInitPhase=true;
85
86// check if string has enough mass to fragment...
87 G4KineticTrackVector * LeftVector=LightFragmentationTest(&theString);
88 if ( LeftVector != 0 ) return LeftVector;
89
90 LeftVector = new G4KineticTrackVector;
91 G4KineticTrackVector * RightVector=new G4KineticTrackVector;
92
93// this should work but its only a semi deep copy. %GF G4ExcitedString theStringInCMS(theString);
94 G4ExcitedString *theStringInCMS=CPExcited(theString);
95 G4LorentzRotation toCms=theStringInCMS->TransformToAlignedCms();
96
97 G4bool success=false, inner_sucess=true;
98 G4int attempt=0;
99 while ( !success && attempt++ < StringLoopInterrupt )
100 {
101 G4FragmentingString *currentString=new G4FragmentingString(*theStringInCMS);
102
103 std::for_each(LeftVector->begin(), LeftVector->end(), DeleteKineticTrack());
104 LeftVector->clear();
105 std::for_each(RightVector->begin(), RightVector->end(), DeleteKineticTrack());
106 RightVector->clear();
107
108 inner_sucess=true; // set false on failure..
109 while (! StopFragmenting(currentString) )
110 { // Split current string into hadron + new string
111 G4FragmentingString *newString=0; // used as output from SplitUp...
112 G4KineticTrack * Hadron=Splitup(currentString,newString);
113 if ( Hadron != 0 && IsFragmentable(newString))
114 {
115 if ( currentString->GetDecayDirection() > 0 )
116 LeftVector->push_back(Hadron);
117 else
118 RightVector->push_back(Hadron);
119 delete currentString;
120 currentString=newString;
121 } else {
122 // abandon ... start from the beginning
123 if (newString) delete newString; // Uzhi restore 20.06.08
124 if (Hadron) delete Hadron;
125 inner_sucess=false;
126 break;
127 }
128 }
129 // Split current string into 2 final Hadrons
130 if ( inner_sucess &&
131 SplitLast(currentString,LeftVector, RightVector) )
132 {
133 success=true;
134 }
135 delete currentString;
136 }
137
138 delete theStringInCMS;
139
140 if ( ! success )
141 {
142 std::for_each(LeftVector->begin(), LeftVector->end(), DeleteKineticTrack());
143 LeftVector->clear();
144 std::for_each(RightVector->begin(), RightVector->end(), DeleteKineticTrack());
145 delete RightVector;
146 return LeftVector;
147 }
148
149 // Join Left- and RightVector into LeftVector in correct order.
150 while(!RightVector->empty())
151 {
152 LeftVector->push_back(RightVector->back());
153 RightVector->erase(RightVector->end()-1);
154 }
155 delete RightVector;
156
157 CalculateHadronTimePosition(theString.Get4Momentum().mag(), LeftVector);
158
159 G4LorentzRotation toObserverFrame(toCms.inverse());
160
161 for(size_t C1 = 0; C1 < LeftVector->size(); C1++)
162 {
163 G4KineticTrack* Hadron = LeftVector->operator[](C1);
164 G4LorentzVector Momentum = Hadron->Get4Momentum();
165 Momentum = toObserverFrame*Momentum;
166 Hadron->Set4Momentum(Momentum);
167 G4LorentzVector Coordinate(Hadron->GetPosition(), Hadron->GetFormationTime());
168 Momentum = toObserverFrame*Coordinate;
169 Hadron->SetFormationTime(Momentum.e());
170 G4ThreeVector aPosition(Momentum.vect());
171 Hadron->SetPosition(theString.GetPosition()+aPosition);
172 }
173 return LeftVector;
174
175
176
177}
178
179//----------------------------------------------------------------------------------------------------------
180
181G4double G4QGSMFragmentation::GetLightConeZ(G4double zmin, G4double zmax, G4int PartonEncoding, G4ParticleDefinition* , G4double , G4double )
182{
183 G4double z;
184 G4double theA(0), d1, d2, yf;
185 G4int absCode = std::abs( PartonEncoding );
186 if (absCode < 10)
187 {
188 if(absCode == 1 || absCode == 2) theA = arho;
189 else if(absCode == 3) theA = aphi;
190 else throw G4HadronicException(__FILE__, __LINE__, "Unknown PDGencoding in G4QGSMFragmentation::G4LightConeZ");
191
192 do
193 {
194 z = zmin + G4UniformRand() * (zmax - zmin);
195 d1 = (1. - z);
196 d2 = (alft - theA);
197 yf = std::pow(d1, d2);
198 }
199 while (G4UniformRand() > yf);
200 }
201 else
202 {
203 if(absCode == 1103 || absCode == 2101 ||
204 absCode == 2203 || absCode == 2103)
205 {
206 d2 = (alft - (2.*an - arho));
207 }
208 else if(absCode == 3101 || absCode == 3103 ||
209 absCode == 3201 || absCode == 3203)
210 {
211 d2 = (alft - (2.*ala - arho));
212 }
213 else
214 {
215 d2 = (alft - (2.*aksi - arho));
216 }
217
218 do
219 {
220 z = zmin + G4UniformRand() * (zmax - zmin);
221 d1 = (1. - z);
222 yf = std::pow(d1, d2);
223 }
224 while (G4UniformRand() > yf);
225 }
226 return z;
227}
228//-----------------------------------------------------------------------------------------
229
230G4LorentzVector * G4QGSMFragmentation::SplitEandP(G4ParticleDefinition * pHadron,
231 G4FragmentingString * string, // Uzhi
232 G4FragmentingString * ) // Uzhi
233{
234 G4double HadronMass = pHadron->GetPDGMass();
235
236 // calculate and assign hadron transverse momentum component HadronPx andHadronPy
237 G4ThreeVector thePt;
238 thePt=SampleQuarkPt();
239 G4ThreeVector HadronPt = thePt +string->DecayPt();
240 HadronPt.setZ(0);
241 //... sample z to define hadron longitudinal momentum and energy
242 //... but first check the available phase space
243 G4double DecayQuarkMass2 = sqr(string->GetDecayParton()->GetPDGMass());
244 G4double HadronMass2T = sqr(HadronMass) + HadronPt.mag2();
245 if (DecayQuarkMass2 + HadronMass2T >= SmoothParam*(string->Mass2()) )
246 return 0; // have to start all over!
247
248 //... then compute allowed z region z_min <= z <= z_max
249
250 G4double zMin = HadronMass2T/(string->Mass2());
251 G4double zMax = 1. - DecayQuarkMass2/(string->Mass2());
252 if (zMin >= zMax) return 0; // have to start all over!
253
254 G4double z = GetLightConeZ(zMin, zMax,
255 string->GetDecayParton()->GetPDGEncoding(), pHadron,
256 HadronPt.x(), HadronPt.y());
257
258 //... now compute hadron longitudinal momentum and energy
259 // longitudinal hadron momentum component HadronPz
260
261 HadronPt.setZ(0.5* string->GetDecayDirection() *
262 (z * string->LightConeDecay() -
263 HadronMass2T/(z * string->LightConeDecay())));
264 G4double HadronE = 0.5* (z * string->LightConeDecay() +
265 HadronMass2T/(z * string->LightConeDecay()));
266
267 G4LorentzVector * a4Momentum= new G4LorentzVector(HadronPt,HadronE);
268
269 return a4Momentum;
270}
271
272
273//-----------------------------------------------------------------------------------------
274
275G4bool G4QGSMFragmentation::SplitLast(G4FragmentingString * string,
276 G4KineticTrackVector * LeftVector,
277 G4KineticTrackVector * RightVector)
278{
279 //... perform last cluster decay
280 G4ThreeVector ClusterVel =string->Get4Momentum().boostVector();
281 G4double ResidualMass =string->Mass();
282 G4double ClusterMassCut = ClusterMass;
283 G4int cClusterInterrupt = 0;
284 G4ParticleDefinition * LeftHadron, * RightHadron;
285 do
286 {
287 if (cClusterInterrupt++ >= ClusterLoopInterrupt)
288 {
289 return false;
290 }
291 G4ParticleDefinition * quark = NULL;
292 string->SetLeftPartonStable(); // to query quark contents..
293 if (string->DecayIsQuark() && string->StableIsQuark() )
294 {
295 //... there are quarks on cluster ends
296 LeftHadron= QuarkSplitup(string->GetLeftParton(), quark);
297 } else {
298 //... there is a Diquark on cluster ends
299 G4int IsParticle;
300 if ( string->StableIsQuark() ) {
301 IsParticle=(string->GetLeftParton()->GetPDGEncoding()>0) ? -1 : +1;
302 } else {
303 IsParticle=(string->GetLeftParton()->GetPDGEncoding()>0) ? +1 : -1;
304 }
305 pDefPair QuarkPair = CreatePartonPair(IsParticle,false); // no diquarks wanted
306 quark = QuarkPair.second;
307 LeftHadron=hadronizer->Build(QuarkPair.first, string->GetLeftParton());
308 }
309 RightHadron = hadronizer->Build(string->GetRightParton(), quark);
310
311 //... repeat procedure, if mass of cluster is too low to produce hadrons
312 //... ClusterMassCut = 0.15*GeV model parameter
313 if ( quark->GetParticleSubType()== "quark" ) {ClusterMassCut = 0.;}
314 else {ClusterMassCut = ClusterMass;}
315 }
316 while (ResidualMass <= LeftHadron->GetPDGMass() + RightHadron->GetPDGMass() + ClusterMassCut);
317
318 //... compute hadron momenta and energies
319 G4LorentzVector LeftMom, RightMom;
320 G4ThreeVector Pos;
321 Sample4Momentum(&LeftMom, LeftHadron->GetPDGMass(), &RightMom, RightHadron->GetPDGMass(), ResidualMass);
322 LeftMom.boost(ClusterVel);
323 RightMom.boost(ClusterVel);
324 LeftVector->push_back(new G4KineticTrack(LeftHadron, 0, Pos, LeftMom));
325 RightVector->push_back(new G4KineticTrack(RightHadron, 0, Pos, RightMom));
326
327 return true;
328
329}
330
331//----------------------------------------------------------------------------------------------------------
332
333G4bool G4QGSMFragmentation::IsFragmentable(const G4FragmentingString * const string)
334{
335 return sqr(FragmentationMass(string)+MassCut) <
336 string->Mass2();
337}
338
339//----------------------------------------------------------------------------------------------------------
340
341G4bool G4QGSMFragmentation::StopFragmenting(const G4FragmentingString * const string)
342{
343 return
344 sqr(FragmentationMass(string,&G4HadronBuilder::BuildHighSpin)+MassCut) >
345 string->Get4Momentum().mag2();
346}
347
348//----------------------------------------------------------------------------------------------------------
349
350//----------------------------------------------------------------------------------------------------------
351
352void G4QGSMFragmentation::Sample4Momentum(G4LorentzVector* Mom, G4double Mass, G4LorentzVector* AntiMom, G4double AntiMass, G4double InitialMass)
353 {
354 G4double r_val = sqr(InitialMass*InitialMass - Mass*Mass - AntiMass*AntiMass) - sqr(2.*Mass*AntiMass);
355 G4double Pabs = (r_val > 0.)? std::sqrt(r_val)/(2.*InitialMass) : 0;
356
357 //... sample unit vector
358 G4double pz = 1. - 2.*G4UniformRand();
359 G4double st = std::sqrt(1. - pz * pz)*Pabs;
360 G4double phi = 2.*pi*G4UniformRand();
361 G4double px = st*std::cos(phi);
362 G4double py = st*std::sin(phi);
363 pz *= Pabs;
364
365 Mom->setPx(px); Mom->setPy(py); Mom->setPz(pz);
366 Mom->setE(std::sqrt(Pabs*Pabs + Mass*Mass));
367
368 AntiMom->setPx(-px); AntiMom->setPy(-py); AntiMom->setPz(-pz);
369 AntiMom->setE (std::sqrt(Pabs*Pabs + AntiMass*AntiMass));
370 }
371
372
373//*********************************************************************************************
Note: See TracBrowser for help on using the repository browser.