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

Last change on this file since 880 was 819, checked in by garnier, 17 years ago

import all except CVS
File size: 13.5 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.6 2007/04/24 14:55:23 gunter Exp $
28// GEANT4 tag $Name: geant4-09-01-patch-02 $
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;
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)
232{
233 G4double HadronMass = pHadron->GetPDGMass();
234
235 // calculate and assign hadron transverse momentum component HadronPx andHadronPy
236 G4ThreeVector thePt;
237 thePt=SampleQuarkPt();
238 G4ThreeVector HadronPt = thePt +string->DecayPt();
239 HadronPt.setZ(0);
240 //... sample z to define hadron longitudinal momentum and energy
241 //... but first check the available phase space
242 G4double DecayQuarkMass2 = sqr(string->GetDecayParton()->GetPDGMass());
243 G4double HadronMass2T = sqr(HadronMass) + HadronPt.mag2();
244 if (DecayQuarkMass2 + HadronMass2T >= SmoothParam*(string->Mass2()) )
245 return 0; // have to start all over!
246
247 //... then compute allowed z region z_min <= z <= z_max
248
249 G4double zMin = HadronMass2T/(string->Mass2());
250 G4double zMax = 1. - DecayQuarkMass2/(string->Mass2());
251 if (zMin >= zMax) return 0; // have to start all over!
252
253 G4double z = GetLightConeZ(zMin, zMax,
254 string->GetDecayParton()->GetPDGEncoding(), pHadron,
255 HadronPt.x(), HadronPt.y());
256
257 //... now compute hadron longitudinal momentum and energy
258 // longitudinal hadron momentum component HadronPz
259
260 HadronPt.setZ(0.5* string->GetDecayDirection() *
261 (z * string->LightConeDecay() -
262 HadronMass2T/(z * string->LightConeDecay())));
263 G4double HadronE = 0.5* (z * string->LightConeDecay() +
264 HadronMass2T/(z * string->LightConeDecay()));
265
266 G4LorentzVector * a4Momentum= new G4LorentzVector(HadronPt,HadronE);
267
268 return a4Momentum;
269}
270
271
272//-----------------------------------------------------------------------------------------
273
274G4bool G4QGSMFragmentation::SplitLast(G4FragmentingString * string,
275 G4KineticTrackVector * LeftVector,
276 G4KineticTrackVector * RightVector)
277{
278 //... perform last cluster decay
279 G4ThreeVector ClusterVel =string->Get4Momentum().boostVector();
280 G4double ResidualMass =string->Mass();
281 G4double ClusterMassCut = ClusterMass;
282 G4int cClusterInterrupt = 0;
283 G4ParticleDefinition * LeftHadron, * RightHadron;
284 do
285 {
286 if (cClusterInterrupt++ >= ClusterLoopInterrupt)
287 {
288 return false;
289 }
290 G4ParticleDefinition * quark = NULL;
291 string->SetLeftPartonStable(); // to query quark contents..
292 if (string->DecayIsQuark() && string->StableIsQuark() )
293 {
294 //... there are quarks on cluster ends
295 LeftHadron= QuarkSplitup(string->GetLeftParton(), quark);
296 } else {
297 //... there is a Diquark on cluster ends
298 G4int IsParticle;
299 if ( string->StableIsQuark() ) {
300 IsParticle=(string->GetLeftParton()->GetPDGEncoding()>0) ? -1 : +1;
301 } else {
302 IsParticle=(string->GetLeftParton()->GetPDGEncoding()>0) ? +1 : -1;
303 }
304 pDefPair QuarkPair = CreatePartonPair(IsParticle,false); // no diquarks wanted
305 quark = QuarkPair.second;
306 LeftHadron=hadronizer->Build(QuarkPair.first, string->GetLeftParton());
307 }
308 RightHadron = hadronizer->Build(string->GetRightParton(), quark);
309
310 //... repeat procedure, if mass of cluster is too low to produce hadrons
311 //... ClusterMassCut = 0.15*GeV model parameter
312 if ( quark->GetParticleSubType()== "quark" ) {ClusterMassCut = 0.;}
313 else {ClusterMassCut = ClusterMass;}
314 }
315 while (ResidualMass <= LeftHadron->GetPDGMass() + RightHadron->GetPDGMass() + ClusterMassCut);
316
317 //... compute hadron momenta and energies
318 G4LorentzVector LeftMom, RightMom;
319 G4ThreeVector Pos;
320 Sample4Momentum(&LeftMom, LeftHadron->GetPDGMass(), &RightMom, RightHadron->GetPDGMass(), ResidualMass);
321 LeftMom.boost(ClusterVel);
322 RightMom.boost(ClusterVel);
323 LeftVector->push_back(new G4KineticTrack(LeftHadron, 0, Pos, LeftMom));
324 RightVector->push_back(new G4KineticTrack(RightHadron, 0, Pos, RightMom));
325
326 return true;
327
328}
329
330//----------------------------------------------------------------------------------------------------------
331
332G4bool G4QGSMFragmentation::IsFragmentable(const G4FragmentingString * const string)
333{
334 return sqr(FragmentationMass(string)+MassCut) <
335 string->Mass2();
336}
337
338//----------------------------------------------------------------------------------------------------------
339
340G4bool G4QGSMFragmentation::StopFragmenting(const G4FragmentingString * const string)
341{
342 return
343 sqr(FragmentationMass(string,&G4HadronBuilder::BuildHighSpin)+MassCut) >
344 string->Get4Momentum().mag2();
345}
346
347//----------------------------------------------------------------------------------------------------------
348
349//----------------------------------------------------------------------------------------------------------
350
351void G4QGSMFragmentation::Sample4Momentum(G4LorentzVector* Mom, G4double Mass, G4LorentzVector* AntiMom, G4double AntiMass, G4double InitialMass)
352 {
353 G4double r_val = sqr(InitialMass*InitialMass - Mass*Mass - AntiMass*AntiMass) - sqr(2.*Mass*AntiMass);
354 G4double Pabs = (r_val > 0.)? std::sqrt(r_val)/(2.*InitialMass) : 0;
355
356 //... sample unit vector
357 G4double pz = 1. - 2.*G4UniformRand();
358 G4double st = std::sqrt(1. - pz * pz)*Pabs;
359 G4double phi = 2.*pi*G4UniformRand();
360 G4double px = st*std::cos(phi);
361 G4double py = st*std::sin(phi);
362 pz *= Pabs;
363
364 Mom->setPx(px); Mom->setPy(py); Mom->setPz(pz);
365 Mom->setE(std::sqrt(Pabs*Pabs + Mass*Mass));
366
367 AntiMom->setPx(-px); AntiMom->setPy(-py); AntiMom->setPz(-pz);
368 AntiMom->setE (std::sqrt(Pabs*Pabs + AntiMass*AntiMass));
369 }
370
371
372//*********************************************************************************************
Note: See TracBrowser for help on using the repository browser.