| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // 26 //
>> 27 // $Id: G4LundStringFragmentation.cc,v 1.23 2010-09-22 12:36:37 vuzhinsk Exp $
>> 28 // GEANT4 tag $Name: not supported by cvs2svn $ 1.8
27 // 29 //
28 // ------------------------------------------- 30 // -----------------------------------------------------------------------------
29 // GEANT 4 class implementation file 31 // GEANT 4 class implementation file
30 // 32 //
31 // History: first implementation, Maxim K 33 // History: first implementation, Maxim Komogorov, 10-Jul-1998
32 // ------------------------------------------- 34 // -----------------------------------------------------------------------------
33 #include "G4LundStringFragmentation.hh" 35 #include "G4LundStringFragmentation.hh"
34 #include "G4PhysicalConstants.hh" <<
35 #include "G4SystemOfUnits.hh" <<
36 #include "Randomize.hh" <<
37 #include "G4FragmentingString.hh" 36 #include "G4FragmentingString.hh"
38 #include "G4DiQuarks.hh" 37 #include "G4DiQuarks.hh"
39 #include "G4Quarks.hh" 38 #include "G4Quarks.hh"
40 #include "G4HadronicParameters.hh" <<
41 #include "G4Exp.hh" <<
42 #include "G4Pow.hh" <<
43 39
44 //#define debug_LUNDfragmentation << 40 #include "Randomize.hh"
45 41
46 // Class G4LundStringFragmentation 42 // Class G4LundStringFragmentation
47 //******************************************** 43 //*************************************************************************************
48 44
49 G4LundStringFragmentation::G4LundStringFragmen 45 G4LundStringFragmentation::G4LundStringFragmentation()
50 : G4VLongitudinalStringDecay("LundStringFrag << 46 {
51 { << 47 // ------ For estimation of a minimal string mass ---------------
52 SetMassCut(210.*MeV); // Mpi + Delta << 48 Mass_of_light_quark =140.*MeV;
53 // For ProduceOneH << 49 Mass_of_heavy_quark =500.*MeV;
54 // that no one pi- << 50 Mass_of_string_junction=720.*MeV;
55 SigmaQT = 0.435 * GeV; << 51 // ------ An estimated minimal string mass ----------------------
56 Tmt = 190.0 * MeV; << 52 MinimalStringMass = 0.;
57 << 53 MinimalStringMass2 = 0.;
58 SetStringTensionParameter(1.*GeV/fermi); << 54 // ------ Minimal invariant mass used at a string fragmentation -
59 SetDiquarkBreakProbability(0.3); << 55 WminLUND = 0.45*GeV; //0.23*GeV; // Uzhi 0.7 -> 0.23 3.8.10 //0.8 1.5
60 << 56 // ------ Smooth parameter used at a string fragmentation for ---
61 SetStrangenessSuppression((1.0 - 0.12)/2.0 << 57 // ------ smearinr sharp mass cut-off ---------------------------
62 SetDiquarkSuppression(0.07); << 58 SmoothParam = 0.2;
63 << 59
64 // Check if charmed and bottom hadrons are << 60 // SetStringTensionParameter(0.25);
65 // set the non-zero probabilities for c-cb << 61 SetStringTensionParameter(1.);
66 // else set them to 0.0. If these probabil << 62
67 // hadrons can't/can be created during the << 63 SetDiquarkBreakProbability(0.05); // Vova Aug. 22
68 // (i.e. not heavy) projectile hadron nucl << 64
69 if ( G4HadronicParameters::Instance()->Ena << 65 // For treating of small string decays
70 SetProbCCbar(0.0002); // According to O << 66 for(G4int i=0; i<3; i++)
71 SetProbBBbar(5.0e-5); // According to O << 67 { for(G4int j=0; j<3; j++)
72 } else { << 68 { for(G4int k=0; k<6; k++)
73 SetProbCCbar(0.0); << 69 { Meson[i][j][k]=0; MesonWeight[i][j][k]=0.;
74 SetProbBBbar(0.0); << 70 }
75 } << 71 }
76 << 72 }
77 SetMinMasses(); // For treating of small << 73 //--------------------------
78 } << 74 Meson[0][0][0]=111; // dbar-d Pi0
>> 75 MesonWeight[0][0][0]=pspin_meson*(1.-scalarMesonMix[0]);
79 76
80 //-------------------------------------------- << 77 Meson[0][0][1]=221; // dbar-d Eta
>> 78 MesonWeight[0][0][1]=pspin_meson*(scalarMesonMix[0]-scalarMesonMix[1]);
81 79
82 G4KineticTrackVector* G4LundStringFragmentatio << 80 Meson[0][0][2]=331; // dbar-d EtaPrime
83 { << 81 MesonWeight[0][0][2]=pspin_meson*(scalarMesonMix[1]);
84 // Can no longer modify Parameters for Fragm <<
85 82
86 PastInitPhase=true; << 83 Meson[0][0][3]=113; // dbar-d Rho0
>> 84 MesonWeight[0][0][3]=(1.-pspin_meson)*(1.-vectorMesonMix[0]);
87 85
88 G4FragmentingString aString(theString); << 86 Meson[0][0][4]=223; // dbar-d Omega
89 SetMinimalStringMass(&aString); << 87 MesonWeight[0][0][4]=(1.-pspin_meson)*(vectorMesonMix[0]);
>> 88 //--------------------------
90 89
91 #ifdef debug_LUNDfragmentation << 90 Meson[0][1][0]=211; // dbar-u Pi+
92 G4cout<<G4endl<<"LUND StringFragmentat << 91 MesonWeight[0][1][0]=pspin_meson;
93 G4cout<<G4endl<<"LUND StringFragm: Str <<
94 <<theString.Get4Momentum <<
95 <<"4Mom "<<theString.Get <<
96 <<"--------------------- <<
97 G4cout<<"String ends Direct "<<theStri <<
98 <<theStri <<
99 <<theStri <<
100 G4cout<<"Left mom "<<theString.GetLef <<
101 G4cout<<"Right mom "<<theString.GetRig <<
102 G4cout<<"Check for Fragmentation "<<G4 <<
103 #endif <<
104 92
105 G4KineticTrackVector * LeftVector(0); << 93 Meson[0][1][1]=213; // dbar-u Rho+
>> 94 MesonWeight[0][1][1]=(1.-pspin_meson);
>> 95 //--------------------------
106 96
107 if (!aString.IsAFourQuarkString() && !IsItFr << 97 Meson[0][2][0]=311; // dbar-s K0bar
108 { << 98 MesonWeight[0][2][0]=pspin_meson;
109 #ifdef debug_LUNDfragmentation <<
110 G4cout<<"Non fragmentable - th <<
111 #endif <<
112 // SetMassCut(210.*MeV); // F <<
113 // t <<
114 <<
115 G4double Mcut = GetMassCut(); <<
116 SetMassCut(10000.*MeV); <<
117 LeftVector=ProduceOneHadron(&theString); <<
118 SetMassCut(Mcut); <<
119 99
120 if ( LeftVector ) << 100 Meson[0][2][1]=313; // dbar-s K*0bar
121 { << 101 MesonWeight[0][2][1]=(1.-pspin_meson);
122 if ( LeftVector->size() > 0) << 102 //--------------------------
123 { << 103 //--------------------------
124 LeftVector->operator[](0)->SetForm << 104 Meson[1][0][0]=211; // ubar-d Pi-
125 LeftVector->operator[](0)->SetPosi << 105 MesonWeight[1][0][0]=pspin_meson;
126 } <<
127 if (LeftVector->size() > 1) <<
128 { <<
129 // 2 hadrons created from qq-qqbar <<
130 LeftVector->operator[](1)->SetFormationT <<
131 LeftVector->operator[](1)->SetPosition(t <<
132 } <<
133 } <<
134 return LeftVector; <<
135 } <<
136 106
137 #ifdef debug_LUNDfragmentation << 107 Meson[1][0][1]=213; // ubar-d Rho-
138 G4cout<<"The string will be fragmented << 108 MesonWeight[1][0][1]=(1.-pspin_meson);
139 #endif << 109 //--------------------------
140 110
141 // The string can fragment. At least two par << 111 Meson[1][1][0]=111; // ubar-u Pi0
142 LeftVector =new G4KineticTrackVec << 112 MesonWeight[1][1][0]=pspin_meson*(1.-scalarMesonMix[0]);
143 G4KineticTrackVector * RightVector=new G4Kin <<
144 113
145 G4bool success = Loop_toFragmentString << 114 Meson[1][1][1]=221; // ubar-u Eta
>> 115 MesonWeight[1][1][1]=pspin_meson*(scalarMesonMix[0]-scalarMesonMix[1]);
146 116
147 if ( ! success ) << 117 Meson[1][1][2]=331; // ubar-u EtaPrime
148 { << 118 MesonWeight[1][1][2]=pspin_meson*(scalarMesonMix[1]);
149 std::for_each(LeftVector->begin(), LeftVec <<
150 LeftVector->clear(); <<
151 std::for_each(RightVector->begin(), RightV <<
152 delete RightVector; <<
153 return LeftVector; <<
154 } <<
155 119
156 // Join Left- and RightVector into LeftVecto << 120 Meson[1][1][3]=113; // ubar-u Rho0
157 while (!RightVector->empty()) << 121 MesonWeight[1][1][3]=(1.-pspin_meson)*(1.-vectorMesonMix[0]);
158 { <<
159 LeftVector->push_back(RightVector->back()) <<
160 RightVector->erase(RightVector->end()-1); <<
161 } <<
162 delete RightVector; <<
163 122
164 return LeftVector; << 123 Meson[1][1][4]=223; // ubar-u Omega
165 } << 124 MesonWeight[1][1][4]=(1.-pspin_meson)*(scalarMesonMix[0]);
>> 125 //--------------------------
166 126
167 //-------------------------------------------- << 127 Meson[1][2][0]=321; // ubar-s K-
>> 128 MesonWeight[1][2][0]=pspin_meson;
168 129
169 G4bool G4LundStringFragmentation::IsItFragment << 130 Meson[1][2][1]=323; // ubar-s K*-bar -
170 { << 131 MesonWeight[1][2][1]=(1.-pspin_meson);
171 SetMinimalStringMass(string); << 132 //--------------------------
172 //G4cout<<"MinM StrM "<<MinimalStringM << 133 //--------------------------
173 134
174 return std::abs(MinimalStringMass) < string- << 135 Meson[2][0][0]=311; // sbar-d K0
>> 136 MesonWeight[2][0][0]=pspin_meson;
175 137
176 //MinimalStringMass is negative and la << 138 Meson[2][0][1]=313; // sbar-d K*0
177 } << 139 MesonWeight[2][0][1]=(1.-pspin_meson);
>> 140 //--------------------------
178 141
179 //-------------------------------------------- << 142 Meson[2][1][0]=321; // sbar-u K+
>> 143 MesonWeight[2][1][0]=pspin_meson;
180 144
181 G4bool G4LundStringFragmentation::Loop_toFragm << 145 Meson[2][1][1]=323; // sbar-u K*+
182 << 146 MesonWeight[2][1][1]=(1.-pspin_meson);
183 << 147 //--------------------------
184 { <<
185 #ifdef debug_LUNDfragmentation <<
186 G4cout<<"Loop_toFrag "<<theString.GetL <<
187 <<theString.GetL <<
188 <<" "<<theString.GetR <<
189 <<theString.GetR <<
190 <<"Direction "<<theString.GetD <<
191 #endif <<
192 148
193 G4LorentzRotation toCmsI, toObserverFr << 149 Meson[2][2][0]=221; // sbar-s Eta
>> 150 MesonWeight[2][2][0]=pspin_meson*(1.-scalarMesonMix[5]);
194 151
195 G4bool final_success=false; << 152 Meson[2][2][1]=331; // sbar-s EtaPrime
196 G4bool inner_success=true; << 153 MesonWeight[2][2][1]=pspin_meson*(1.-scalarMesonMix[5]);
197 154
198 G4int attempt=0; << 155 Meson[2][2][3]=333; // sbar-s EtaPrime
>> 156 MesonWeight[2][2][3]=(1.-pspin_meson)*(vectorMesonMix[5]);
>> 157 //--------------------------
199 158
200 while ( ! final_success && attempt++ < Strin << 159 for(G4int i=0; i<3; i++)
201 { // If the string fragmentation does << 160 { for(G4int j=0; j<3; j++)
202 // repeat the fragmentation. << 161 { for(G4int k=0; k<3; k++)
>> 162 { for(G4int l=0; l<4; l++)
>> 163 { Baryon[i][j][k][l]=0; BaryonWeight[i][j][k][l]=0.;}
>> 164 }
>> 165 }
>> 166 }
203 167
204 G4FragmentingString *currentSt << 168 //---------------------------------------
205 toCmsI = currentString->Transf << 169 Baryon[0][0][0][0]=1114; // Delta-
206 toObserverFrameI = toCmsI.inve << 170 BaryonWeight[0][0][0][0]=1.;
207 171
208 G4LorentzRotation toCms, toObserverFrame; << 172 //---------------------------------------
>> 173 Baryon[0][0][1][0]=2112; // neutron
>> 174 BaryonWeight[0][0][1][0]=pspin_barion;
209 175
210 //G4cout<<"Main loop start whilecounter "< << 176 Baryon[0][0][1][1]=2114; // Delta0
>> 177 BaryonWeight[0][0][1][1]=(1.-pspin_barion);
211 178
212 // Cleaning up the previously produced had << 179 //---------------------------------------
213 std::for_each(LeftVector->begin(), LeftVec << 180 Baryon[0][0][2][0]=3112; // Sigma-
214 LeftVector->clear(); << 181 BaryonWeight[0][0][2][0]=pspin_barion;
215 std::for_each(RightVector->begin(), RightV <<
216 RightVector->clear(); <<
217 182
218 // Main fragmentation loop until the strin << 183 Baryon[0][0][2][1]=3114; // Sigma*-
219 inner_success=true; // set false on failu << 184 BaryonWeight[0][0][2][1]=(1.-pspin_barion);
220 const G4int maxNumberOfLoops = <<
221 G4int loopCounter = -1; <<
222 <<
223 while ( (! StopFragmenting(currentString)) <<
224 { // Split current string into hadro <<
225 #ifdef debug_LUNDfragm <<
226 G4cout<<"The string wi <<
227 //G4cout<<"1 "<<curren <<
228 #endif <<
229 G4FragmentingString *newString=0; // us <<
230 <<
231 toCms=currentString->TransformToAlignedC <<
232 toObserverFrame= toCms <<
233 <<
234 #ifdef debug_LUNDfragm <<
235 //G4cout<<"CMS Left m <<
236 //G4cout<<"CMS Right m <<
237 //G4cout<<"CMS String <<
238 #endif <<
239 185
240 G4KineticTrack * Hadron=Splitup(currentS << 186 //---------------------------------------
>> 187 Baryon[0][1][0][0]=2112; // neutron
>> 188 BaryonWeight[0][1][0][0]=pspin_barion;
241 189
242 if ( Hadron != 0 ) // Store the hadron << 190 Baryon[0][1][0][1]=2114; // Delta0
243 { << 191 BaryonWeight[0][1][0][1]=(1.-pspin_barion);
244 #ifdef debug_L <<
245 G4cout<<"Hadro <<
246 //G4cout<<"2 " <<
247 #endif <<
248 <<
249 Hadron->Set4Momentum(toObserverFrame*H <<
250 <<
251 G4double TimeOftheStringCreation=theSt <<
252 G4ThreeVector PositionOftheStringCreat <<
253 <<
254 G4LorentzVector Coordinate(Hadron->Get <<
255 G4LorentzVector Momentum = toObserverF <<
256 Hadron->SetFormationTime(TimeOftheStri <<
257 G4ThreeVector aPosition(Momentum.vect( <<
258 Hadron->SetPosition(PositionOftheStrin <<
259 <<
260 // Open to pro <<
261 if ( currentString->GetDecayDirection( <<
262 { <<
263 LeftVector->push_back(Hadron); <<
264 } else <<
265 { <<
266 RightVector->push_back(Hadron); <<
267 } <<
268 delete currentString; <<
269 currentString=newString; <<
270 } else { <<
271 if ( newString ) del <<
272 } <<
273 <<
274 currentString->Lorentz <<
275 }; <<
276 <<
277 if ( loopCounter >= maxNumberO <<
278 inner_success=false; <<
279 } <<
280 <<
281 // Split remaining string into 2 final had <<
282 #ifdef debug_LUNDfragmentation <<
283 if (inner_success) G4cout<<"Sp <<
284 #endif <<
285 192
286 if ( inner_success && SplitLast(currentStr << 193 //---------------------------------------
287 { << 194 Baryon[0][1][1][0]=2212; // proton
288 final_success = true; << 195 BaryonWeight[0][1][1][0]=pspin_barion;
289 } <<
290 196
291 delete currentString; << 197 Baryon[0][1][1][1]=2214; // Delta+
292 } // End of the loop where we try to fragme << 198 BaryonWeight[0][1][1][1]=(1.-pspin_barion);
293 199
294 G4int sign = +1; << 200 //---------------------------------------
295 if ( theString.GetDirection() < 0 ) si << 201 Baryon[0][1][2][0]=3122; // Lambda
296 for ( unsigned int hadronI = 0; hadron << 202 BaryonWeight[0][1][2][0]=pspin_barion*0.5;
297 G4LorentzVector Tmp = LeftVector->o <<
298 Tmp.setZ(sign*Tmp.getZ()); <<
299 Tmp *= toObserverFrameI; <<
300 LeftVector->operator[](hadronI)->Se <<
301 } <<
302 for ( unsigned int hadronI = 0; hadron <<
303 G4LorentzVector Tmp = RightVector-> <<
304 Tmp.setZ(sign*Tmp.getZ()); <<
305 Tmp *= toObserverFrameI; <<
306 RightVector->operator[](hadronI)->S <<
307 } <<
308 203
309 return final_success; << 204 Baryon[0][1][2][1]=3212; // Sigma0
310 } << 205 BaryonWeight[0][1][2][2]=pspin_barion*0.5;
311 206
312 //-------------------------------------------- << 207 Baryon[0][1][2][2]=3214; // Sigma*0
>> 208 BaryonWeight[0][1][2][2]=(1.-pspin_barion);
313 209
314 G4bool G4LundStringFragmentation::StopFragment << 210 //---------------------------------------
315 { << 211 Baryon[0][2][0][0]=3112; // Sigma-
316 SetMinimalStringMass(string); << 212 BaryonWeight[0][2][0][0]=pspin_barion;
317 213
318 if ( MinimalStringMass < 0.) return true; << 214 Baryon[0][2][0][1]=3114; // Sigma*-
>> 215 BaryonWeight[0][2][0][1]=(1.-pspin_barion);
319 216
320 if (string->IsAFourQuarkString()) << 217 //---------------------------------------
321 { << 218 Baryon[0][2][1][0]=3122; // Lambda
322 return G4UniformRand() < G4Exp(-0.0005*(st << 219 BaryonWeight[0][2][1][0]=pspin_barion*0.5;
323 } else { <<
324 <<
325 if (MinimalStringMass < 0.0 ) <<
326 <<
327 G4bool Result = G4UniformRand() < <<
328 G4Exp(-0.66e-6*(string->Mass()*string- <<
329 // G4bool Result = string->Mas <<
330 <<
331 #ifdef debug_LUNDfragmentation <<
332 G4cout<<"StopFragmenting Minim <<
333 <<" "<<string->Mass()<<G <<
334 G4cout<<"StopFragmenting - Yes <<
335 #endif <<
336 return Result; <<
337 } <<
338 } <<
339 220
340 //-------------------------------------------- << 221 Baryon[0][2][1][1]=3212; // Sigma0
>> 222 BaryonWeight[0][2][1][1]=pspin_barion*0.5;
341 223
342 G4KineticTrack * G4LundStringFragmentation::Sp << 224 Baryon[0][2][1][2]=3214; // Sigma*0
343 G4Fragmentin << 225 BaryonWeight[0][2][1][2]=(1.-pspin_barion);
344 { <<
345 #ifdef debug_LUNDfragmentation <<
346 G4cout<<G4endl; <<
347 G4cout<<"Start SplitUP ================ <<
348 G4cout<<"String partons: " <<string->Ge <<
349 <<string->Ge <<
350 <<"Direction " <<string->Ge <<
351 #endif <<
352 <<
353 //... random choice of string end to us <<
354 G4int SideOfDecay = (G4UniformRand() < <<
355 if (SideOfDecay < 0) <<
356 { <<
357 string->SetLeftPartonStable(); <<
358 } else <<
359 { <<
360 string->SetRightPartonStable(); <<
361 } <<
362 <<
363 G4ParticleDefinition *newStringEnd; <<
364 G4ParticleDefinition * HadronDefinition <<
365 <<
366 G4double StringMass=string->Mass(); <<
367 <<
368 G4double ProbDqADq = GetDiquarkSuppress <<
369 G4double ProbSaS = 1.0 - 2.0 * GetStr <<
370 <<
371 #ifdef debug_LUNDfragmentation <<
372 G4cout<<"StrMass DiquarkSuppression <<
373 #endif <<
374 <<
375 G4int NumberOfpossibleBaryons = 2; <<
376 <<
377 if (string->GetLeftParton()->GetParticl <<
378 if (string->GetRightParton()->GetPartic <<
379 <<
380 G4double ActualProb = ProbDqADq ; <<
381 ActualProb *= (1.0-G4Pow::GetInstance() <<
382 if(ActualProb <0.0) ActualProb = 0.; <<
383 <<
384 SetDiquarkSuppression(ActualProb); <<
385 <<
386 G4double Mth = 1250.0; <<
387 if ( NumberOfpossibleBaryons == 3 ){Mth <<
388 else if ( NumberOfpossibleBaryons == 4 <<
389 else {} <<
390 <<
391 ActualProb = ProbSaS; <<
392 ActualProb *= (1.0 - G4Pow::GetInstance <<
393 if ( ActualProb < 0.0 ) ActualProb = 0. <<
394 SetStrangenessSuppression((1.0-ActualPr <<
395 <<
396 #ifdef debug_LUNDfragmentation <<
397 G4cout<<"StrMass DiquarkSuppression cor <<
398 #endif <<
399 <<
400 if (string->DecayIsQuark()) <<
401 { <<
402 HadronDefinition= QuarkSplitup(strin <<
403 } else { <<
404 HadronDefinition= DiQuarkSplitup(str <<
405 } <<
406 <<
407 SetDiquarkSuppression(ProbDqADq); <<
408 SetStrangenessSuppression((1.0-ProbSaS) <<
409 <<
410 if ( HadronDefinition == NULL ) { G4Kin <<
411 <<
412 #ifdef debug_LUNDfragmentation <<
413 G4cout<<"The parton "<<string->GetDecay <<
414 <<" produces hadron "<<HadronDefi <<
415 <<" and is transformed to "<<newS <<
416 G4cout<<"The side of the string decay L <<
417 #endif <<
418 // create new String from old, ie. keep <<
419 <<
420 if ( newString ) delete newString; <<
421 <<
422 newString=new G4FragmentingString(*stri <<
423 <<
424 #ifdef debug_LUNDfragmentation <<
425 G4cout<<"An attempt to determine its en <<
426 #endif <<
427 G4LorentzVector* HadronMomentum=SplitEa <<
428 226
429 delete newString; newString=0; << 227 //---------------------------------------
430 << 228 Baryon[0][2][2][0]=3312; // Theta-
431 G4KineticTrack * Hadron =0; << 229 BaryonWeight[0][2][2][0]=pspin_barion;
432 if ( HadronMomentum != 0 ) { <<
433 230
434 #ifdef debug_LUNDfragmentation << 231 Baryon[0][2][2][1]=3314; // Theta*-
435 G4cout<<"The attempt was successful << 232 BaryonWeight[0][2][2][1]=(1.-pspin_barion);
436 #endif <<
437 G4ThreeVector Pos; <<
438 Hadron = new G4KineticTrack(HadronDefinit <<
439 <<
440 if ( newString ) delete newString; <<
441 <<
442 newString=new G4FragmentingString(*string <<
443 HadronMomentum); <<
444 delete HadronMomentum; <<
445 } <<
446 else <<
447 { <<
448 #ifdef debug_LUNDfragmentation <<
449 G4cout<<"The attempt was not succes <<
450 #endif <<
451 } <<
452 <<
453 #ifdef debug_LUNDfragmentation <<
454 G4cout<<"End SplitUP (G4VLongitudinalSt <<
455 #endif <<
456 233
457 return Hadron; << 234 //---------------------------------------
458 } << 235 //---------------------------------------
>> 236 Baryon[1][0][0][0]=2112; // neutron
>> 237 BaryonWeight[1][0][0][0]=pspin_barion;
459 238
460 //-------------------------------------------- << 239 Baryon[1][0][0][1]=2114; // Delta0
>> 240 BaryonWeight[1][0][0][1]=(1.-pspin_barion);
461 241
462 G4ParticleDefinition * G4LundStringFragmentati << 242 //---------------------------------------
463 << 243 Baryon[1][0][1][0]=2212; // proton
464 { << 244 BaryonWeight[1][0][1][0]=pspin_barion;
465 G4double StrSup=GetStrangeSuppress(); <<
466 G4double ProbQQbar = (1.0 - 2.0*StrSup)*1.2 <<
467 245
468 //... can Diquark break or not? << 246 Baryon[1][0][1][1]=2214; // Delta+
469 if (G4UniformRand() < DiquarkBreakProb ){ << 247 BaryonWeight[1][0][1][1]=(1.-pspin_barion);
470 248
471 //... Diquark break << 249 //---------------------------------------
472 G4int stableQuarkEncoding = decay->GetPD << 250 Baryon[1][0][2][0]=3122; // Lambda
473 G4int decayQuarkEncoding = (decay->GetPD << 251 BaryonWeight[1][0][2][0]=pspin_barion*0.5;
474 if (G4UniformRand() < 0.5) <<
475 { <<
476 G4int Swap = stableQuarkEncoding; <<
477 stableQuarkEncoding = decayQuarkEncod <<
478 decayQuarkEncoding = Swap; <<
479 } <<
480 252
481 G4int IsParticle=(decayQuarkEncoding>0) << 253 Baryon[1][0][2][1]=3212; // Sigma0
>> 254 BaryonWeight[1][0][2][1]=pspin_barion*0.5;
482 255
483 SetStrangenessSuppression((1.0-ProbQQbar << 256 Baryon[1][0][2][2]=3214; // Sigma*0
484 pDefPair QuarkPair = CreatePartonPair(Is << 257 BaryonWeight[1][0][2][2]=(1.-pspin_barion);
485 SetStrangenessSuppression((1.0-StrSup)/2 <<
486 258
487 //... Build new Diquark << 259 //---------------------------------------
488 G4int QuarkEncoding=QuarkPair.second->Ge << 260 Baryon[1][1][0][0]=2212; // proton
489 G4int i10 = std::max(std::abs(QuarkEnco << 261 BaryonWeight[1][1][0][0]=pspin_barion;
490 G4int i20 = std::min(std::abs(QuarkEnco <<
491 G4int spin = (i10 != i20 && G4UniformRan <<
492 G4int NewDecayEncoding = -1*IsParticle*( <<
493 created = FindParticle(NewDecayEncoding) <<
494 G4ParticleDefinition * decayQuark=FindPa <<
495 G4ParticleDefinition * had=hadronizer->B <<
496 StrangeSuppress=StrSup; <<
497 262
498 return had; << 263 Baryon[1][1][0][1]=2214; // Delta+
>> 264 BaryonWeight[1][1][0][1]=(1.-pspin_barion);
499 265
500 } else { << 266 //---------------------------------------
501 //... Diquark does not break << 267 Baryon[1][1][1][0]=2224; // Delta++
>> 268 BaryonWeight[1][1][1][0]=1.;
502 269
503 G4int IsParticle=(decay->GetPDGEncoding( << 270 //---------------------------------------
>> 271 Baryon[1][1][2][0]=3222; // Sigma+
>> 272 BaryonWeight[1][1][2][0]=pspin_barion;
504 273
505 StrangeSuppress=(1.0 - ProbQQbar)/2.0; << 274 Baryon[1][1][2][1]=3224; // Sigma*+
506 pDefPair QuarkPair = CreatePartonPair(Is << 275 BaryonWeight[1][1][2][1]=(1.-pspin_barion);
507 276
508 created = QuarkPair.second; << 277 //---------------------------------------
>> 278 Baryon[1][2][0][0]=3122; // Lambda
>> 279 BaryonWeight[1][2][0][0]=pspin_barion*0.5;
509 280
510 G4ParticleDefinition * had=hadronizer->B << 281 Baryon[1][2][0][1]=3212; // Sigma0
511 StrangeSuppress=StrSup; << 282 BaryonWeight[1][2][0][1]=pspin_barion*0.5;
512 283
513 return had; << 284 Baryon[1][2][0][2]=3214; // Sigma*0
514 } << 285 BaryonWeight[1][2][0][2]=(1.-pspin_barion);
515 } <<
516 286
517 //-------------------------------------------- << 287 //---------------------------------------
>> 288 Baryon[1][2][1][0]=3222; // Sigma+
>> 289 BaryonWeight[1][2][1][0]=pspin_barion;
518 290
519 G4LorentzVector * G4LundStringFragmentation::S << 291 Baryon[1][2][1][1]=3224; // Sigma*+
520 << 292 BaryonWeight[1][2][1][1]=(1.-pspin_barion);
521 <<
522 { <<
523 G4LorentzVector String4Momentum=string->Get4 <<
524 G4double StringMT2=string->MassT2(); <<
525 G4double StringMT =std::sqrt(StringMT2); <<
526 <<
527 G4double HadronMass = pHadron->GetPDGMass(); <<
528 SetMinimalStringMass(newString); <<
529 <<
530 if ( MinimalStringMass < 0.0 ) return n <<
531 <<
532 #ifdef debug_LUNDfragmentation <<
533 G4cout<<G4endl<<"Start LUND SplitEandP <<
534 G4cout<<"String 4 mom, String M and Mt <<
535 <<" "<<std::sqrt(StringMT2)<<G4e <<
536 G4cout<<"Hadron "<<pHadron->GetParticl <<
537 G4cout<<"HadM MinimalStringMassLeft St <<
538 <<String4Momentum.mag()<<" "<<Ha <<
539 #endif <<
540 293
541 if ((HadronMass + MinimalStringMass > string << 294 //---------------------------------------
542 { << 295 Baryon[1][2][2][0]=3322; // Theta0
543 #ifdef debug_LUNDfragmentation << 296 BaryonWeight[1][2][2][0]=pspin_barion;
544 G4cout<<"Mass of the string is not s <<
545 #endif <<
546 return 0; <<
547 } // have to start all over! <<
548 <<
549 String4Momentum.setPz(0.); <<
550 G4ThreeVector StringPt=String4Momentum.vect( <<
551 StringPt.setZ(0.); <<
552 <<
553 // calculate and assign hadron transverse mo <<
554 G4ThreeVector HadronPt , RemSysPt; <<
555 G4double HadronMassT2, ResidualMassT2; <<
556 G4double HadronMt, Pt, Pt2, phi; <<
557 <<
558 G4double TmtCur = Tmt; <<
559 <<
560 if ( (string->GetDecayParton()->GetPar <<
561 (pHadron->GetBaryonNumber() != 0) <<
562 TmtCur = Tmt*0.37; // q <<
563 } else if ( (string->GetDecayParton()- <<
564 (pHadron->GetBaryonNumber( <<
565 //TmtCur = Tmt; <<
566 } else if ( (string->GetDecayParton()->GetPa <<
567 (pHadron->GetBaryonNumber( <<
568 //TmtCur = Tmt*0.89; <<
569 } else if ( (string->GetDecayParton()- <<
570 (pHadron->GetBaryonNumber( <<
571 TmtCur = Tmt*1.35; <<
572 } <<
573 297
574 //... sample Pt of the hadron << 298 Baryon[1][2][2][1]=3324; // Theta*0
575 G4int attempt=0; << 299 BaryonWeight[1][2][2][1]=(1.-pspin_barion);
576 do <<
577 { <<
578 attempt++; if (attempt > StringLoopI <<
579 300
580 HadronMt = HadronMass - TmtCur*G4Log << 301 //---------------------------------------
581 Pt2 = sqr(HadronMt)-sqr(HadronMass); Pt=st << 302 //---------------------------------------
582 phi = 2.*pi*G4UniformRand(); << 303 Baryon[2][0][0][0]=3112; // Sigma-
583 HadronPt = G4ThreeVector( Pt*std::co << 304 BaryonWeight[2][0][0][0]=pspin_barion;
584 RemSysPt = StringPt - HadronPt; <<
585 HadronMassT2 = sqr(HadronMass) + Had <<
586 ResidualMassT2=sqr(MinimalStringMass <<
587 305
588 } while (std::sqrt(HadronMassT2) + std << 306 Baryon[2][0][0][1]=3114; // Sigma*-
>> 307 BaryonWeight[2][0][0][1]=(1.-pspin_barion);
589 308
590 //... sample z to define hadron longitudina << 309 //---------------------------------------
591 //... but first check the available phase sp << 310 Baryon[2][0][1][0]=3122; // Lambda
>> 311 BaryonWeight[2][0][1][0]=pspin_barion*0.5;
592 312
593 G4double Pz2 = (sqr(StringMT2 - HadronMassT2 << 313 Baryon[2][0][1][1]=3212; // Sigma0
594 4*HadronMassT2 * ResidualMassT2)/4./Stri << 314 BaryonWeight[2][0][1][1]=pspin_barion*0.5;
595 315
596 if (Pz2 < 0 ) {return 0;} // have t << 316 Baryon[2][0][1][2]=3214; // Sigma*0
>> 317 BaryonWeight[2][0][1][2]=(1.-pspin_barion);
597 318
598 //... then compute allowed z region z_min < << 319 //---------------------------------------
>> 320 Baryon[2][0][2][0]=3312; // Sigma-
>> 321 BaryonWeight[2][0][2][0]=pspin_barion;
599 322
600 G4double Pz = std::sqrt(Pz2); << 323 Baryon[2][0][2][1]=3314; // Sigma*-
601 G4double zMin = (std::sqrt(HadronMassT2+Pz2) << 324 BaryonWeight[2][0][2][1]=(1.-pspin_barion);
602 // G4double zMin = (std::sqrt(HadronMa <<
603 G4double zMax = (std::sqrt(HadronMassT2+Pz2) <<
604 325
605 if (zMin >= zMax) return 0; // have to sta << 326 //---------------------------------------
>> 327 Baryon[2][1][0][0]=3122; // Lambda
>> 328 BaryonWeight[2][1][0][0]=pspin_barion*0.5;
606 329
607 G4double z = GetLightConeZ(zMin, zMax, << 330 Baryon[2][1][0][1]=3212; // Sigma0
608 string->GetDecayParton()->Get << 331 BaryonWeight[2][1][0][1]=pspin_barion*0.5;
609 HadronPt.x(), HadronPt.y()); <<
610 332
611 //... now compute hadron longitudinal moment << 333 Baryon[2][1][0][2]=3214; // Sigma*0
612 // longitudinal hadron momentum component Ha << 334 BaryonWeight[2][1][0][2]=(1.-pspin_barion);
613 335
614 HadronPt.setZ(0.5* string->GetDecayDirection << 336 //---------------------------------------
615 (z * string->LightConeDecay() - << 337 Baryon[2][1][1][0]=3222; // Sigma+
616 HadronMassT2/(z * string->LightConeD << 338 BaryonWeight[2][1][1][0]=pspin_barion;
617 G4double HadronE = 0.5* (z * string->LightC <<
618 HadronMassT2/(z * string->LightConeD <<
619 339
620 G4LorentzVector * a4Momentum= new G4LorentzV << 340 Baryon[2][1][1][1]=3224; // Sigma*+
>> 341 BaryonWeight[2][1][1][1]=(1.-pspin_barion);
621 342
622 #ifdef debug_LUNDfragmentation << 343 //---------------------------------------
623 G4cout<<G4endl<<" string->GetDecayDire << 344 Baryon[2][1][2][0]=3322; // Theta0
624 G4cout<<"string->LightConeDecay() "<<s << 345 BaryonWeight[2][1][2][0]=pspin_barion;
625 G4cout<<"HadronPt,HadronE "<<HadronPt< <<
626 G4cout<<"String4Momentum "<<String4Mom <<
627 G4cout<<"Out of LUND SplitEandP "<<G4e <<
628 #endif <<
629 346
630 return a4Momentum; << 347 Baryon[2][1][2][1]=3324; // Theta*0
631 } << 348 BaryonWeight[2][1][2][2]=(1.-pspin_barion);
632 349
633 //-------------------------------------------- << 350 //---------------------------------------
>> 351 Baryon[2][2][0][0]=3312; // Theta-
>> 352 BaryonWeight[2][2][0][0]=pspin_barion;
634 353
635 G4double G4LundStringFragmentation::GetLightCo << 354 Baryon[2][2][0][1]=3314; // Theta*-
636 G4int PD << 355 BaryonWeight[2][2][0][1]=(1.-pspin_barion);
637 G4Partic <<
638 G4double <<
639 { <<
640 G4double Mass = pHadron->GetPDGMass(); <<
641 G4int HadronEncoding=std::abs(pHadron- <<
642 356
643 G4double Mt2 = Px*Px + Py*Py + Mass*Mass; << 357 //---------------------------------------
>> 358 Baryon[2][2][1][0]=3322; // Theta0
>> 359 BaryonWeight[2][2][1][0]=pspin_barion;
644 360
645 G4double Alund, Blund; << 361 Baryon[2][2][1][1]=3324; // Theta*0
646 G4double zOfMaxyf(0.), maxYf(1.), z(0.), yf( << 362 BaryonWeight[2][2][1][1]=(1.-pspin_barion);
647 363
648 if (!((std::abs(PDGEncodingOfDecayParton) > << 364 //---------------------------------------
649 { // ---------------- Quark fragmentation << 365 Baryon[2][2][2][0]=3334; // Omega
650 Alund=1.; << 366 BaryonWeight[2][2][2][0]=1.;
651 Blund=0.7/GeV/GeV; <<
652 <<
653 G4double BMt2 = Blund*Mt2; <<
654 if (Alund == 1.0) { <<
655 zOfMaxyf=BMt2/(Blund*Mt2 + 1.);} <<
656 else { <<
657 zOfMaxyf = ((1.0+BMt2) - std::sqrt(sq <<
658 } <<
659 <<
660 if (zOfMaxyf < zmin) {zOfMaxyf=zmin;} <<
661 if (zOfMaxyf > zmax) {zOfMaxyf=zmax;} <<
662 maxYf=(1-zOfMaxyf)/zOfMaxyf * G4Exp(-Blun <<
663 <<
664 const G4int maxNumberOfLoops = 1000 <<
665 G4int loopCounter = 0; <<
666 do <<
667 { <<
668 z = zmin + G4UniformRand()*(zmax-zmin); <<
669 //yf = (1-z)/z * G4Exp(-Blund* <<
670 yf = G4Pow::GetInstance()->powA(1.0-z,Alun <<
671 } <<
672 while ( (G4UniformRand()*maxYf > yf) && + <<
673 if ( loopCounter >= maxNumberOfLoop <<
674 z = 0.5*(zmin + zmax); // Just a <<
675 } <<
676 return z; <<
677 } <<
678 367
679 if (std::abs(PDGEncodingOfDecayParton) > 100 << 368 //---------------------------------------
680 { << 369 /*
681 G4double an = 2.5; << 370 for(G4int i=0; i<3; i++)
682 an +=(sqr(Px)+sqr(Py))/sqr(GeV << 371 { for(G4int j=0; j<3; j++)
683 z=zmin + (zmax-zmin)*G4Pow::Ge << 372 { for(G4int k=0; k<3; k++)
684 if( PDGEncodingOfDecayParton > << 373 { for(G4int l=0; l<4; l++)
685 } << 374 { G4cout<<i<<" "<<j<<" "<<k<<" "<<l<<" "<<Baryon[i][j][k][l]<<G4endl;}
>> 375 }
>> 376 }
>> 377 }
>> 378 G4int Uzhi;
>> 379 G4cin>>Uzhi;
>> 380 */
>> 381 //StrangeSuppress=0.38;
>> 382 Prob_QQbar[0]=StrangeSuppress; // Probability of ddbar production
>> 383 Prob_QQbar[1]=StrangeSuppress; // Probability of uubar production
>> 384 Prob_QQbar[2]=StrangeSuppress/(2.+StrangeSuppress);//(1.-2.*StrangeSuppress); // Probability of ssbar production
>> 385 }
686 386
687 return z; << 387 // --------------------------------------------------------------
688 } << 388 G4LundStringFragmentation::G4LundStringFragmentation(const G4LundStringFragmentation &) : G4VLongitudinalStringDecay()
>> 389 {
>> 390 }
689 391
690 //-------------------------------------------- << 392 G4LundStringFragmentation::~G4LundStringFragmentation()
>> 393 {
>> 394 }
691 395
692 G4bool G4LundStringFragmentation::SplitLast(G4 << 396 //*************************************************************************************
693 G4 << 397
694 G4 << 398 const G4LundStringFragmentation & G4LundStringFragmentation::operator=(const G4LundStringFragmentation &)
>> 399 {
>> 400 throw G4HadronicException(__FILE__, __LINE__, "G4LundStringFragmentation::operator= meant to not be accessable");
>> 401 return *this;
>> 402 }
>> 403
>> 404 int G4LundStringFragmentation::operator==(const G4LundStringFragmentation &right) const
>> 405 {
>> 406 return !memcmp(this, &right, sizeof(G4LundStringFragmentation));
>> 407 }
>> 408
>> 409 int G4LundStringFragmentation::operator!=(const G4LundStringFragmentation &right) const
>> 410 {
>> 411 return memcmp(this, &right, sizeof(G4LundStringFragmentation));
>> 412 }
>> 413
>> 414 //--------------------------------------------------------------------------------------
>> 415 void G4LundStringFragmentation::SetMinimalStringMass(const G4FragmentingString * const string)
695 { 416 {
696 //... perform last cluster decay << 417 G4double EstimatedMass=0.;
697 SetMinimalStringMass( string); << 418 G4int Number_of_quarks=0;
698 if ( MinimalStringMass < 0.) return fa <<
699 #ifdef debug_LUNDfragmentation <<
700 G4cout<<G4endl<<"Split last----------- <<
701 G4cout<<"MinimalStringMass "<<MinimalS <<
702 G4cout<<"Left "<<string->GetLeftParto <<
703 G4cout<<"Right "<<string->GetRightPart <<
704 G4cout<<"String4mom "<<string->GetPstr <<
705 #endif <<
706 <<
707 G4LorentzVector Str4Mom=string->Get4Mo <<
708 G4LorentzRotation toCms(-1*Str4Mom.boo <<
709 G4LorentzVector Pleft = toCms * string <<
710 toCms.rotateZ(-1*Pleft.phi()); <<
711 toCms.rotateY(-1*Pleft.theta()); <<
712 <<
713 G4LorentzRotation toObserverFrame= toC <<
714 419
715 G4double StringMass=string->Mass(); << 420 G4int Qleft =std::abs(string->GetLeftParton()->GetPDGEncoding());
716 421
717 G4ParticleDefinition * LeftHadron(0), * Righ << 422 if( Qleft > 1000)
>> 423 {
>> 424 Number_of_quarks+=2;
>> 425 G4int q1=Qleft/1000;
>> 426 if( q1 < 3) {EstimatedMass +=Mass_of_light_quark;}
>> 427 if( q1 > 2) {EstimatedMass +=Mass_of_heavy_quark;}
>> 428
>> 429 G4int q2=(Qleft/100)%10;
>> 430 if( q2 < 3) {EstimatedMass +=Mass_of_light_quark;}
>> 431 if( q2 > 2) {EstimatedMass +=Mass_of_heavy_quark;}
>> 432 EstimatedMass +=Mass_of_string_junction;
>> 433 }
>> 434 else
>> 435 {
>> 436 Number_of_quarks++;
>> 437 if( Qleft < 3) {EstimatedMass +=Mass_of_light_quark;}
>> 438 if( Qleft > 2) {EstimatedMass +=Mass_of_heavy_quark;}
>> 439 }
718 440
719 NumberOf_FS=0; << 441 G4int Qright=std::abs(string->GetRightParton()->GetPDGEncoding());
720 for (G4int i=0; i<350; i++) {FS_Weight[i]=0. <<
721 442
722 G4int sampledState = 0; << 443 if( Qright > 1000)
>> 444 {
>> 445 Number_of_quarks+=2;
>> 446 G4int q1=Qright/1000;
>> 447 if( q1 < 3) {EstimatedMass +=Mass_of_light_quark;}
>> 448 if( q1 > 2) {EstimatedMass +=Mass_of_heavy_quark;}
>> 449
>> 450 G4int q2=(Qright/100)%10;
>> 451 if( q2 < 3) {EstimatedMass +=Mass_of_light_quark;}
>> 452 if( q2 > 2) {EstimatedMass +=Mass_of_heavy_quark;}
>> 453 EstimatedMass +=Mass_of_string_junction;
>> 454 }
>> 455 else
>> 456 {
>> 457 Number_of_quarks++;
>> 458 if( Qright < 3) {EstimatedMass +=Mass_of_light_quark;}
>> 459 if( Qright > 2) {EstimatedMass +=Mass_of_heavy_quark;}
>> 460 }
723 461
724 #ifdef debug_LUNDfragmentation << 462 if(Number_of_quarks==2){EstimatedMass +=100.*MeV;}
725 G4cout<<"StrMass "<<StringMass<<" q's << 463 if(Number_of_quarks==3){EstimatedMass += 20.*MeV;}
726 <<string->GetLeftParton()->GetPa << 464 if(Number_of_quarks==4){EstimatedMass -=2.*Mass_of_string_junction;
727 <<string->GetRightParton()->GetP << 465 if(EstimatedMass <= 1600.*MeV){EstimatedMass-=200.*MeV;}
728 #endif << 466 else {EstimatedMass+=100.*MeV;}
>> 467 }
>> 468 MinimalStringMass=EstimatedMass;
>> 469 SetMinimalStringMass2(EstimatedMass);
>> 470 }
729 471
730 string->SetLeftPartonStable(); // to query q << 472 //--------------------------------------------------------------------------------------
>> 473 void G4LundStringFragmentation::SetMinimalStringMass2(
>> 474 const G4double aValue)
>> 475 {
>> 476 MinimalStringMass2=aValue * aValue;
>> 477 }
731 478
732 if (string->IsAFourQuarkString() ) << 479 //--------------------------------------------------------------------------------------
733 { << 480 G4KineticTrackVector* G4LundStringFragmentation::FragmentString(
734 G4int IDleft =std::abs(string->GetLe << 481 const G4ExcitedString& theString)
735 G4int IDright=std::abs(string->GetRi << 482 {
>> 483 // Can no longer modify Parameters for Fragmentation.
>> 484 PastInitPhase=true;
>> 485 //SetVectorMesonProbability(1.);
>> 486 //SetSpinThreeHalfBarionProbability(1.);
736 487
737 if ( (IDleft > 3000) || (IDright > 3 << 488 // check if string has enough mass to fragment...
738 if ( ! Diquark_AntiDiquark_belowTh << 489
739 { << 490 SetMassCut(160.*MeV); // For LightFragmentationTest it is required
740 return false; << 491 // that no one pi-meson can be produced
741 } << 492 //
742 } else { << 493 //G4cout<<G4endl<<"G4LundStringFragmentation::"<<G4endl;
743 // The string is qq-qqbar type. Diquarks a << 494 //G4cout<<"FragmentString Position"<<theString.GetPosition()/fermi<<" "<<theString.GetTimeOfCreation()/fermi<<G4endl;
744 if (StringMass-MinimalStringMass < 0 << 495 //G4cout<<"FragmentString Momentum"<<theString.Get4Momentum()<<theString.Get4Momentum().mag()<<G4endl;
745 { << 496 //
746 if (! Diquark_AntiDiquark_belowThreshold << 497 G4FragmentingString aString(theString);
747 { << 498 SetMinimalStringMass(&aString);
748 return false; <<
749 } <<
750 } else <<
751 { <<
752 Diquark_AntiDiquark_aboveThreshold_lastS <<
753 499
754 if (NumberOf_FS == 0) return false; << 500 /*
>> 501 G4cout<<"St Frag "<<MinimalStringMass<<" "<<WminLUND<<" "<<(1.-SmoothParam)<<" "<< theString.Get4Momentum().mag()<<G4endl;
>> 502 G4cout<<(MinimalStringMass+WminLUND)*(1.-SmoothParam)<<" "<<theString.Get4Momentum().mag()<<G4endl;
>> 503
>> 504 if((MinimalStringMass+WminLUND)*(1.-SmoothParam) > theString.Get4Momentum().mag())
>> 505 {SetMassCut(1000.*MeV);}
>> 506 else {SetMassCut(160.*MeV);}
>> 507 */
>> 508 // V.U. 20.06.10 in order to put in correspondence LightFragTest and MinStrMass
>> 509
>> 510 G4KineticTrackVector * LeftVector(0);
>> 511 // G4KineticTrackVector * LeftVector=LightFragmentationTest(&theString);
>> 512 //G4cout<<"Min Str Mass "<<MinimalStringMass<<G4endl;
>> 513 if(!IsFragmentable(&aString)) // produce 1 hadron
>> 514 {
>> 515 //G4cout<<"Non fragmentable"<<G4endl;
>> 516 SetMassCut(1000.*MeV);
>> 517 LeftVector=LightFragmentationTest(&theString);
>> 518 SetMassCut(160.*MeV);
>> 519
>> 520
>> 521 //G4cout<<"Prod hadron "<<LeftVector->operator[](0)->GetDefinition()->GetParticleName()<<G4endl;
>> 522 /*
>> 523 if(LeftVector->size() == 1)
>> 524 {
>> 525 if(! (std::abs(LeftVector->operator[](0)->GetDefinition()->GetPDGMass()-
>> 526 theString.Get4Momentum().mag()) < 100*MeV))
>> 527 { // produce a particle with arbitrary isospin
>> 528 G4cout<<" produce a particle with arbitrary isospin"<<G4endl;
>> 529
>> 530 pDefPair hadrons((G4ParticleDefinition *)0,(G4ParticleDefinition *)0);
>> 531 Pcreate build=&G4HadronBuilder::Build;
>> 532 FragmentationMass(&aString,build,&hadrons); // 0->1
>> 533 G4cout<<"Hadron PDG "<<hadrons.first->GetPDGEncoding()<<G4endl;
>> 534 if ( hadrons.second ==0 )
>> 535 {// Substitute string by light hadron, Note that Energy is not conserved here!
>> 536 // Cleaning up the previously produced hadrons ------------------------------
>> 537 std::for_each(LeftVector->begin() , LeftVector->end() , DeleteKineticTrack());
>> 538 LeftVector->clear();
>> 539
>> 540 G4ThreeVector Mom3 = aString.Get4Momentum().vect();
>> 541 G4LorentzVector Mom(Mom3,std::sqrt(Mom3.mag2() + sqr(hadrons.first->GetPDGMass())));
>> 542 LeftVector->push_back(new G4KineticTrack(hadrons.first, 0,
>> 543 theString.GetPosition(),
>> 544 Mom));
>> 545 } // end of if ( hadrons.second ==0 )
>> 546 } // end of produce a particle with arbitrary isospin
>> 547
>> 548 } // end of if(LeftVector->size() == 1)
>> 549 */
>> 550 } // end of if(!IsFragmentable(&aString))
>> 551
>> 552 if ( LeftVector != 0 ) {
>> 553 // Uzhi insert 6.05.08 start
>> 554 if(LeftVector->size() == 1){
>> 555 // One hadron is saved in the interaction
>> 556 LeftVector->operator[](0)->SetFormationTime(theString.GetTimeOfCreation());
>> 557 LeftVector->operator[](0)->SetPosition(theString.GetPosition());
>> 558
>> 559 /* // To set large formation time open *
>> 560 LeftVector->operator[](0)->SetFormationTime(theString.GetTimeOfCreation()+100.*fermi);
>> 561 LeftVector->operator[](0)->SetPosition(theString.GetPosition());
>> 562 G4ThreeVector aPosition(theString.GetPosition().x(),
>> 563 theString.GetPosition().y(),
>> 564 theString.GetPosition().z()+100.*fermi);
>> 565 LeftVector->operator[](0)->SetPosition(aPosition);
>> 566 */
>> 567 } else { // 2 hadrons created from qq-qqbar are stored
>> 568
>> 569 LeftVector->operator[](0)->SetFormationTime(theString.GetTimeOfCreation());
>> 570 LeftVector->operator[](0)->SetPosition(theString.GetPosition());
>> 571 LeftVector->operator[](1)->SetFormationTime(theString.GetTimeOfCreation());
>> 572 LeftVector->operator[](1)->SetPosition(theString.GetPosition());
>> 573 }
>> 574 return LeftVector;
>> 575 }
755 576
756 sampledState = SampleS << 577 //--------------------- The string can fragment -------------------------------
757 if (string->GetLeftParton()->GetPDGEncod << 578 //--------------- At least two particles can be produced ----------------------
758 { << 579 //G4cout<<"Fragmentable"<<G4endl;
759 LeftHadron =FS_LeftHadron[sampledState << 580 LeftVector =new G4KineticTrackVector;
760 RightHadron=FS_RightHadron[sampledStat << 581 G4KineticTrackVector * RightVector=new G4KineticTrackVector;
761 } else <<
762 { <<
763 LeftHadron =FS_RightHadron[sampledStat <<
764 RightHadron=FS_LeftHadron[sampledState <<
765 } <<
766 } <<
767 } // ID > 3300 <<
768 } else { <<
769 if (string->DecayIsQuark() && string->Stab <<
770 { //... there are quarks on cluster <<
771 #ifdef debug_LUNDfragm <<
772 G4cout<<"Q Q string La <<
773 #endif <<
774 582
775 Quark_AntiQuark_lastSplitting(string, Le << 583 G4ExcitedString *theStringInCMS=CPExcited(theString);
>> 584 G4LorentzRotation toCms=theStringInCMS->TransformToAlignedCms();
>> 585
>> 586 G4bool success=false, inner_sucess=true;
>> 587 G4int attempt=0;
>> 588 while ( !success && attempt++ < StringLoopInterrupt )
>> 589 { // If the string fragmentation do not be happend, repeat the fragmentation---
>> 590 G4FragmentingString *currentString=new G4FragmentingString(*theStringInCMS);
>> 591 //G4cout<<"Main loop start whilecounter "<<attempt<<G4endl;
>> 592 // Cleaning up the previously produced hadrons ------------------------------
>> 593 std::for_each(LeftVector->begin() , LeftVector->end() , DeleteKineticTrack());
>> 594 LeftVector->clear();
>> 595
>> 596 std::for_each(RightVector->begin(), RightVector->end(), DeleteKineticTrack());
>> 597 RightVector->clear();
776 598
777 if (NumberOf_FS == 0) return false; << 599 // Main fragmentation loop until the string will not be able to fragment ----
778 sampledState = SampleState() << 600 inner_sucess=true; // set false on failure..
779 if (string->GetLeftParton()->GetPDGEncod <<
780 { <<
781 LeftHadron =FS_RightHadron[sampledStat <<
782 RightHadron=FS_LeftHadron[sampledState <<
783 } else <<
784 { <<
785 LeftHadron =FS_LeftHadron[sampledState <<
786 RightHadron=FS_RightHadron[sampledStat <<
787 } <<
788 } else <<
789 { //... there is a Diquark on one of <<
790 #ifdef debug_LUNDfragm <<
791 G4cout<<"DiQ Q string <<
792 #endif <<
793 601
794 Quark_Diquark_lastSplitting(string, Left << 602 while (! StopFragmenting(currentString) )
>> 603 { // Split current string into hadron + new string
795 604
796 if (NumberOf_FS == 0) return false; << 605 G4FragmentingString *newString=0; // used as output from SplitUp...
797 sampledState = SampleState() <<
798 606
799 if (string->GetLeftParton()->GetParticle << 607 G4KineticTrack * Hadron=Splitup(currentString,newString);
800 { << 608 //G4cout<<"SplitUp------"<<Hadron<<G4endl;
801 LeftHadron =FS_LeftHadron[sampledState << 609
802 RightHadron=FS_RightHadron[sampledStat << 610 if ( Hadron != 0 ) // Store the hadron
803 } else <<
804 { 611 {
805 LeftHadron =FS_RightHadron[sampledStat << 612 //G4cout<<"Hadron prod at fragm. "<<Hadron->GetDefinition()->GetParticleName()<<G4endl;
806 RightHadron=FS_LeftHadron[sampledState << 613 if ( currentString->GetDecayDirection() > 0 )
807 } << 614 LeftVector->push_back(Hadron);
>> 615 else
>> 616 RightVector->push_back(Hadron);
>> 617 delete currentString;
>> 618 currentString=newString;
>> 619 }
>> 620 };
>> 621 // Split remaining string into 2 final Hadrons ------------------------
>> 622 //G4cout<<"Split Last"<<G4endl;
>> 623 if ( inner_sucess &&
>> 624 SplitLast(currentString,LeftVector, RightVector) )
>> 625 {
>> 626 success=true;
808 } 627 }
809 << 628 delete currentString;
>> 629 } // End of the loop in attemps to fragment the string
>> 630
>> 631 delete theStringInCMS;
>> 632
>> 633 if ( ! success )
>> 634 {
>> 635 std::for_each(LeftVector->begin(), LeftVector->end(), DeleteKineticTrack());
>> 636 LeftVector->clear();
>> 637 std::for_each(RightVector->begin(), RightVector->end(), DeleteKineticTrack());
>> 638 delete RightVector;
>> 639 return LeftVector;
810 } 640 }
811 << 641
812 #ifdef debug_LUNDfragmentation << 642 // Join Left- and RightVector into LeftVector in correct order.
813 G4cout<<"Sampled hadrons: "<<LeftHadro << 643 while(!RightVector->empty())
814 #endif << 644 {
815 << 645 LeftVector->push_back(RightVector->back());
816 G4LorentzVector P_left =string->GetPleft(), << 646 RightVector->erase(RightVector->end()-1);
817 <<
818 G4LorentzVector LeftMom, RightMom; <<
819 G4ThreeVector Pos; <<
820 <<
821 Sample4Momentum(&LeftMom, LeftHadron->GetPD <<
822 &RightMom, RightHadron->GetPDGMass(), <<
823 StringMass); <<
824 <<
825 // Sample4Momentum ascribes LeftMom.pz <<
826 // It must be negative in case when th <<
827 <<
828 if (!(string->DecayIsQuark() && string->Stab <<
829 { // Only for qq - q, q - qq, and qq - qqbar <<
830 <<
831 if ( G4UniformRand() <= 0.5 ) <<
832 { <<
833 if (P_left.z() <= 0.) {G4LorentzVector t <<
834 } <<
835 else <<
836 { <<
837 if (P_right.z() >= 0.) {G4LorentzVector <<
838 } <<
839 } 647 }
>> 648 delete RightVector;
840 649
841 LeftMom *=toObserverFrame; << 650 CalculateHadronTimePosition(theString.Get4Momentum().mag(), LeftVector);
842 RightMom*=toObserverFrame; <<
843 651
844 LeftVector->push_back(new G4KineticTrack(Lef << 652 G4LorentzRotation toObserverFrame(toCms.inverse());
845 RightVector->push_back(new G4KineticTrack(Ri <<
846 653
847 string->LorentzRotate(toObserverFrame); << 654 // LeftVector->operator[](0)->SetFormationTime(theString.GetTimeOfCreation());
848 return true; << 655 // LeftVector->operator[](0)->SetPosition(theString.GetPosition());
849 } <<
850 656
851 //-------------------------------------------- << 657 G4double TimeOftheStringCreation=theString.GetTimeOfCreation();
852 << 658 G4ThreeVector PositionOftheStringCreation(theString.GetPosition());
853 G4bool G4LundStringFragmentation:: <<
854 Diquark_AntiDiquark_belowThreshold_lastSplitti <<
855 <<
856 <<
857 { <<
858 G4double StringMass = string->Mass(); <<
859 659
860 G4int cClusterInterrupt = 0; << 660 //G4cout<<"# prod hadrons "<<LeftVector->size()<<G4endl;
861 G4bool isOK = false; << 661 for(size_t C1 = 0; C1 < LeftVector->size(); C1++)
862 do <<
863 { 662 {
864 G4int LeftQuark1= string->GetLeftParton()- << 663 G4KineticTrack* Hadron = LeftVector->operator[](C1);
865 G4int LeftQuark2=(string->GetLeftParton()- << 664 G4LorentzVector Momentum = Hadron->Get4Momentum();
>> 665 //G4cout<<"Hadron "<<Hadron->GetDefinition()->GetParticleName()<<" "<<Momentum<<G4endl;
>> 666 Momentum = toObserverFrame*Momentum;
>> 667 Hadron->Set4Momentum(Momentum);
>> 668
>> 669 G4LorentzVector Coordinate(Hadron->GetPosition(), Hadron->GetFormationTime());
>> 670 Momentum = toObserverFrame*Coordinate;
>> 671 Hadron->SetFormationTime(TimeOftheStringCreation+Momentum.e()
>> 672 -fermi/c_light);
>> 673 G4ThreeVector aPosition(Momentum.vect());
>> 674 // Hadron->SetPosition(theString.GetPosition()+aPosition);
>> 675 Hadron->SetPosition(PositionOftheStringCreation+aPosition);
>> 676 };
866 677
867 G4int RightQuark1= string->GetRightParton( << 678 return LeftVector;
868 G4int RightQuark2=(string->GetRightParton( << 679
>> 680 }
869 681
870 if (G4UniformRand()<0.5) << 682 //----------------------------------------------------------------------------------
871 { << 683 G4bool G4LundStringFragmentation::IsFragmentable(const G4FragmentingString * const string)
872 LeftHadron =hadronizer->Build(FindPartic << 684 {
873 FindParticle(RightQuark1)); << 685 SetMinimalStringMass(string);
874 RightHadron= (LeftHadron == nullptr) ? n << 686 // return sqr(MinimalStringMass + WminLUND) < string->Get4Momentum().mag2();
875 << 687 return MinimalStringMass < string->Get4Momentum().mag(); // 21.07.2010
876 FindParticle(RightQuark2)); << 688 }
877 } else <<
878 { <<
879 LeftHadron =hadronizer->Build(FindPartic <<
880 FindParticle(RightQuark2)); <<
881 RightHadron=(LeftHadron == nullptr) ? nu <<
882 hadronizer <<
883 FindParticle(RightQuark1)); <<
884 } <<
885 689
886 isOK = (LeftHadron != nullptr) && (RightHa << 690 //----------------------------------------------------------------------------------------
>> 691 G4bool G4LundStringFragmentation::StopFragmenting(const G4FragmentingString * const string)
>> 692 {
>> 693 SetMinimalStringMass(string);
887 694
888 if(isOK) { isOK = (StringMass > LeftHadron << 695 /*
889 ++cClusterInterrupt; << 696 G4cout<<"StopFragm MinMass "<<MinimalStringMass<<" String Mass "<<string->Get4Momentum().mag()<<G4endl;
890 //... repeat procedure, if mass of cluster << 697 G4cout<<"WminLUND "<<WminLUND<<" SmoothParam "<<SmoothParam<<" "<<string->Mass()<<G4endl;
891 //... ClusterMassCut = 0.15*GeV model para << 698 //G4cout<<std::exp(-(string->Mass()*string->Mass()-MinimalStringMass2)/WminLUND/WminLUND)<<G4endl;
892 } << 699 //G4int Uzhi; G4cin>>Uzhi;
893 while (isOK == false && cClusterInterrupt < << 700 */
894 /* Loop checking, 07.08.2015, A.Ribon */ << 701 //
895 return isOK; << 702 return (MinimalStringMass + WminLUND)*
>> 703 (1 + SmoothParam * (1.-2*G4UniformRand())) >
>> 704 string->Mass();
>> 705 //
>> 706 // return G4UniformRand() < std::exp(-(string->Mass()*string->Mass()-MinimalStringMass2)/WminLUND/WminLUND);
896 } 707 }
897 708
898 //-------------------------------------------- 709 //----------------------------------------------------------------------------------------
899 << 710 G4bool G4LundStringFragmentation::SplitLast(G4FragmentingString * string,
900 G4bool G4LundStringFragmentation:: << 711 G4KineticTrackVector * LeftVector,
901 Diquark_AntiDiquark_aboveThreshold_lastSplitti << 712 G4KineticTrackVector * RightVector)
902 <<
903 <<
904 { 713 {
905 // StringMass-MinimalStringMass > 0. Creatio << 714 //... perform last cluster decay
>> 715 //G4cout<<"Split last-----------------------------------------"<<G4endl;
>> 716 G4LorentzVector Str4Mom=string->Get4Momentum();
906 717
907 G4double StringMass = string->Mass(); << 718 G4ThreeVector ClusterVel =string->Get4Momentum().boostVector();
908 G4double StringMassSqr= sqr(StringMass); <<
909 G4ParticleDefinition * Di_Quark; <<
910 G4ParticleDefinition * Anti_Di_Quark; <<
911 719
912 if (string->GetLeftParton()->GetPDGEncoding( << 720 G4double StringMass = string->Mass();
913 { << 721 G4double StringMassSqr= sqr(StringMass);
914 Anti_Di_Quark =string->GetLeftParton(); <<
915 Di_Quark=string->GetRightParton(); <<
916 } else <<
917 { <<
918 Anti_Di_Quark =string->GetRightParton(); <<
919 Di_Quark=string->GetLeftParton(); <<
920 } <<
921 722
922 G4int IDAnti_di_quark =Anti_Di_Quark->Get << 723 G4ParticleDefinition * LeftHadron(0), * RightHadron(0);
923 G4int AbsIDAnti_di_quark =std::abs(IDAnti_di << 724 G4double LeftHadronMass(0.), RightHadronMass(0.);
924 G4int IDdi_quark =Di_Quark->GetPDGEn <<
925 G4int AbsIDdi_quark =std::abs(IDdi_quar <<
926 725
927 G4int ADi_q1=AbsIDAnti_di_quark/1000; << 726 G4ParticleDefinition * FS_LeftHadron[25], * FS_RightHadron[25];
928 G4int ADi_q2=(AbsIDAnti_di_quark-ADi_q1*1000 << 727 G4double FS_Weight[25];
>> 728 G4int NumberOf_FS=0;
929 729
930 G4int Di_q1=AbsIDdi_quark/1000; << 730 for(G4int i=0; i<25; i++) {FS_Weight[i]=0.;}
931 G4int Di_q2=(AbsIDdi_quark-Di_q1*1000)/100; << 731 //***********************************************
>> 732 //G4cout<<"StrMass "<<StringMass<<" q "<<string->GetLeftParton()->GetParticleName()<<" "<<string->GetRightParton()->GetParticleName()<<" StringMassSqr "<<StringMassSqr<<G4endl;
932 733
933 NumberOf_FS=0; <<
934 for (G4int ProdQ=1; ProdQ < 6; ProdQ++) <<
935 { <<
936 G4int StateADiQ=0; <<
937 const G4int maxNumberOfLoops = <<
938 G4int loopCounter = 0; <<
939 do // while(Meson[AbsIDquark-1][ProdQ-1][ <<
940 { <<
941 LeftHadron=G4ParticleTable::GetParticleT <<
942 -Baryon[ADi_q1-1][ADi_q2-1][Prod <<
943 734
944 if (LeftHadron == NULL) continue; << 735 string->SetLeftPartonStable(); // to query quark contents..
945 G4double LeftHadronMass=LeftHadron->GetP <<
946 736
947 G4int StateDiQ=0; << 737 if (string->FourQuarkString() )
948 const G4int maxNumberO << 738 {
949 G4int internalLoopCoun << 739 // The string is qq-qqbar type. Diquarks are on the string ends
950 do // while(Baryon[Di_q1-1][Di_q2-1][Pro << 740 G4int cClusterInterrupt = 0;
951 { << 741 do
952 RightHadron=G4ParticleTable::GetPartic << 742 {
953 +Baryon[Di_q1-1][Di_q2-1][Prod << 743 //G4cout<<"cClusterInterrupt "<<cClusterInterrupt<<G4endl;
>> 744 if (cClusterInterrupt++ >= ClusterLoopInterrupt)
>> 745 {
>> 746 return false;
>> 747 }
954 748
955 if (RightHadron == NULL) continue; << 749 G4int LeftQuark1= string->GetLeftParton()->GetPDGEncoding()/1000;
956 G4double RightHadronMass=RightHadron-> << 750 G4int LeftQuark2=(string->GetLeftParton()->GetPDGEncoding()/100)%10;
957 751
958 if (StringMass > LeftHadronMass + Righ << 752 G4int RightQuark1= string->GetRightParton()->GetPDGEncoding()/1000;
959 { << 753 G4int RightQuark2=(string->GetRightParton()->GetPDGEncoding()/100)%10;
960 if ( N <<
961 G4Ex <<
962 ed < <<
963 G4Ex <<
964 <<
965 Numb <<
966 } <<
967 <<
968 G4double FS_Psqr=lambda(StringMassSq <<
969 sqr(RightHadronMass)); <<
970 //FS_Psqr=1.; <<
971 FS_Weight[NumberOf_FS]=std::sqrt(FS_ <<
972 BaryonWeight[ADi_q1-1][AD <<
973 BaryonWeight[Di_q1-1][Di_ <<
974 Prob_QQbar[ProdQ-1]; <<
975 <<
976 FS_LeftHadron[NumberOf_FS] = LeftHad <<
977 FS_RightHadron[NumberOf_FS]= RightHa <<
978 NumberOf_FS++; <<
979 } // End of if (StringMass > LeftHadro <<
980 <<
981 StateDiQ++; <<
982 <<
983 } while( (Baryon[Di_q1-1][Di_q2-1][ProdQ <<
984 ++internalLoo <<
985 if ( internalLoopCount <<
986 return false; <<
987 } <<
988 <<
989 StateADiQ++; <<
990 } while( (Baryon[ADi_q1-1][ADi_q2-1][ProdQ <<
991 ++loopCounter < maxNu <<
992 if ( loopCounter >= maxNumberO <<
993 return false; <<
994 } <<
995 } // End of for (G4int ProdQ=1; ProdQ < 4; P <<
996 754
997 return true; << 755 if(G4UniformRand()<0.5)
998 } << 756 {
>> 757 LeftHadron =hadronizer->Build(FindParticle( LeftQuark1),
>> 758 FindParticle(RightQuark1));
>> 759 RightHadron=hadronizer->Build(FindParticle( LeftQuark2),
>> 760 FindParticle(RightQuark2));
>> 761 } else
>> 762 {
>> 763 LeftHadron =hadronizer->Build(FindParticle( LeftQuark1),
>> 764 FindParticle(RightQuark2));
>> 765 RightHadron=hadronizer->Build(FindParticle( LeftQuark2),
>> 766 FindParticle(RightQuark1));
>> 767 }
>> 768
>> 769 //... repeat procedure, if mass of cluster is too low to produce hadrons
>> 770 //... ClusterMassCut = 0.15*GeV model parameter
>> 771 }
>> 772 while ((StringMass <= LeftHadron->GetPDGMass() + RightHadron->GetPDGMass()));
>> 773
>> 774 } // End of if (string->FourQuarkString() )
>> 775
>> 776 //***********************************************
>> 777
>> 778 if (!string->FourQuarkString())
>> 779 {
>> 780 if (string->DecayIsQuark() && string->StableIsQuark() )
>> 781 {//... there are quarks on cluster ends
>> 782 //G4cout<<"Q Q string"<<G4endl;
>> 783 G4ParticleDefinition * Quark;
>> 784 G4ParticleDefinition * Anti_Quark;
>> 785
>> 786 if(string->GetLeftParton()->GetPDGEncoding()>0)
>> 787 {
>> 788 Quark =string->GetLeftParton();
>> 789 Anti_Quark=string->GetRightParton();
>> 790 } else
>> 791 {
>> 792 Quark =string->GetRightParton();
>> 793 Anti_Quark=string->GetLeftParton();
>> 794 }
999 795
1000 //------------------------------------------- << 796 G4int IDquark =Quark->GetPDGEncoding();
>> 797 G4int AbsIDquark =std::abs(IDquark);
>> 798 G4int IDanti_quark =Anti_Quark->GetPDGEncoding();
>> 799 G4int AbsIDanti_quark=std::abs(IDanti_quark);
1001 800
1002 G4bool G4LundStringFragmentation::Quark_Diqua << 801 NumberOf_FS=0;
1003 << 802 for(G4int ProdQ=1; ProdQ < 4; ProdQ++)
1004 << 803 {
1005 { << 804 G4int SignQ=-1;
1006 G4double StringMass = string->Mass(); << 805 if(IDquark == 2) SignQ= 1;
1007 G4double StringMassSqr= sqr(StringMass); << 806 if((IDquark == 1) && (ProdQ == 3)) SignQ= 1; // K0
>> 807 if((IDquark == 3) && (ProdQ == 1)) SignQ=-1; // K0bar
>> 808 if(IDquark == ProdQ) SignQ= 1;
>> 809
>> 810 G4int SignAQ= 1;
>> 811 if(IDanti_quark == -2) SignAQ=-1;
>> 812 if((IDanti_quark ==-1) && (ProdQ == 3)) SignAQ=-1; // K0bar
>> 813 if((IDanti_quark ==-3) && (ProdQ == 1)) SignAQ= 1; // K0
>> 814 if(AbsIDanti_quark == ProdQ) SignAQ= 1;
>> 815
>> 816 G4int StateQ=0;
>> 817 do // while(Meson[AbsIDquark-1][ProdQ-1][StateQ]<>0);
>> 818 {
>> 819 LeftHadron=G4ParticleTable::GetParticleTable()->FindParticle(SignQ*
>> 820 Meson[AbsIDquark-1][ProdQ-1][StateQ]);
>> 821 LeftHadronMass=LeftHadron->GetPDGMass();
>> 822 StateQ++;
>> 823
>> 824 G4int StateAQ=0;
>> 825 do // while(Meson[AbsIDanti_quark-1][ProdQ-1][StateAQ]<>0);
>> 826 {
>> 827 RightHadron=G4ParticleTable::GetParticleTable()->FindParticle(SignAQ*
>> 828 Meson[AbsIDanti_quark-1][ProdQ-1][StateAQ]);
>> 829 RightHadronMass=RightHadron->GetPDGMass();
>> 830 StateAQ++;
>> 831
>> 832 if(StringMass >= LeftHadronMass + RightHadronMass)
>> 833 {
>> 834 G4double FS_Psqr=lambda(StringMassSqr,sqr(LeftHadronMass),
>> 835 sqr(RightHadronMass));
>> 836 FS_Weight[NumberOf_FS]=std::sqrt(FS_Psqr)*
>> 837 MesonWeight[AbsIDquark-1][ProdQ-1][StateQ]*
>> 838 MesonWeight[AbsIDanti_quark-1][ProdQ-1][StateAQ]*
>> 839 Prob_QQbar[ProdQ-1];
>> 840
>> 841 FS_LeftHadron[NumberOf_FS] = LeftHadron;
>> 842 FS_RightHadron[NumberOf_FS]= RightHadron;
>> 843 NumberOf_FS++;
>> 844 if(NumberOf_FS > 24)
>> 845 {G4int Uzhi; G4cout<<"QQbar string #_FS "<<NumberOf_FS<<G4endl; G4cin>>Uzhi;}
>> 846 } // End of if(StringMass >= LeftHadronMass + RightHadronMass)
>> 847 } while(Meson[AbsIDanti_quark-1][ProdQ-1][StateAQ]!=0);
>> 848 } while(Meson[AbsIDquark-1][ProdQ-1][StateQ]!=0);
>> 849 } // End of for(G4int ProdQ=1; ProdQ < 4; ProdQ++)
>> 850 //
>> 851 } else //---------------------------------------------
>> 852 { //... there is a Diquark on one of the cluster ends
>> 853 //G4cout<<"DiQ Q string"<<G4endl;
>> 854 G4ParticleDefinition * Di_Quark;
>> 855 G4ParticleDefinition * Quark;
>> 856
>> 857 if(string->GetLeftParton()->GetParticleSubType()== "quark")
>> 858 {
>> 859 Quark =string->GetLeftParton();
>> 860 Di_Quark=string->GetRightParton();
>> 861 } else
>> 862 {
>> 863 Quark =string->GetRightParton();
>> 864 Di_Quark=string->GetLeftParton();
>> 865 }
1008 866
1009 G4ParticleDefinition * Di_Quark; << 867 G4int IDquark =Quark->GetPDGEncoding();
1010 G4ParticleDefinition * Quark; << 868 G4int AbsIDquark =std::abs(IDquark);
>> 869 G4int IDdi_quark =Di_Quark->GetPDGEncoding();
>> 870 G4int AbsIDdi_quark=std::abs(IDdi_quark);
>> 871 G4int Di_q1=AbsIDdi_quark/1000;
>> 872 G4int Di_q2=(AbsIDdi_quark-Di_q1*1000)/100;
>> 873 //G4cout<<"IDs "<<IDdi_quark<<" "<<IDquark<<G4endl;
1011 874
1012 if (string->GetLeftParton()->GetParticleSub << 875 G4int SignDiQ= 1;
1013 { << 876 if(IDdi_quark < 0) SignDiQ=-1;
1014 Quark =string->GetLeftParton(); <<
1015 Di_Quark=string->GetRightParton(); <<
1016 } else <<
1017 { <<
1018 Quark =string->GetRightParton(); <<
1019 Di_Quark=string->GetLeftParton(); <<
1020 } <<
1021 877
1022 G4int IDquark =Quark->GetPDGEncoding << 878 NumberOf_FS=0;
1023 G4int AbsIDquark =std::abs(IDquark); << 879 for(G4int ProdQ=1; ProdQ < 4; ProdQ++)
1024 G4int IDdi_quark =Di_Quark->GetPDGEncodin << 880 {
1025 G4int AbsIDdi_quark=std::abs(IDdi_quark); << 881 G4int SignQ;
1026 G4int Di_q1=AbsIDdi_quark/1000; << 882 if(IDquark > 0)
1027 G4int Di_q2=(AbsIDdi_quark-Di_q1*1000)/100; << 883 { SignQ=-1;
1028 G4int SignDiQ= 1; << 884 if(IDquark == 2) SignQ= 1;
1029 if (IDdi_quark < 0) SignDiQ=-1; << 885 if((IDquark == 1) && (ProdQ == 3)) SignQ= 1; // K0
1030 << 886 if((IDquark == 3) && (ProdQ == 1)) SignQ=-1; // K0bar
1031 NumberOf_FS=0; << 887 } else
1032 for (G4int ProdQ=1; ProdQ < 4; ProdQ++) // << 888 {
1033 { // << 889 SignQ= 1;
1034 G4int SignQ; << 890 if(IDquark == -2) SignQ=-1;
1035 if (IDquark > 0) << 891 if((IDquark ==-1) && (ProdQ == 3)) SignQ=-1; // K0bar
1036 { << 892 if((IDquark ==-3) && (ProdQ == 1)) SignQ= 1; // K0
1037 SignQ=-1; << 893 }
1038 if (IDquark == 2) << 894
1039 if ((IDquark == 1) && (ProdQ == 3)) Sig << 895 if(AbsIDquark == ProdQ) SignQ= 1;
1040 if ((IDquark == 3) && (ProdQ == 1)) Sig << 896
1041 if (IDquark == 4) << 897 //G4cout<<G4endl;
1042 if (IDquark == 5) Sig << 898 //G4cout<<"-------------------------------------------"<<G4endl;
1043 } else << 899 //G4cout<<"Insert QQbar "<<ProdQ<<" Sign "<<SignQ<<G4endl;
1044 { << 900
1045 SignQ= 1; << 901 G4int StateQ=0;
1046 if (IDquark == -2) Sig << 902 do // while(Meson[AbsIDquark-1][ProdQ-1][StateQ]<>0);
1047 if ((IDquark ==-1) && (ProdQ == 3)) Sig << 903 {
1048 if ((IDquark ==-3) && (ProdQ == 1)) Sig << 904 //G4cout<<G4endl<<"AbsIDquark ProdQ StateQ "<<MesonWeight[AbsIDquark-1][ProdQ-1][StateQ]<<" "<<AbsIDquark<<" "<<ProdQ<<" "<<StateQ<<G4endl;
1049 if (IDquark == -4) Sig << 905 //G4cout<<G4endl<<"AbsIDquark ProdQ StateQ "<<SignQ*Meson[AbsIDquark-1][ProdQ-1][StateQ]<<" "<<AbsIDquark<<" "<<ProdQ<<" "<<StateQ<<G4endl;
1050 if (IDquark == -5) Sig << 906 LeftHadron=G4ParticleTable::GetParticleTable()->FindParticle(SignQ*
1051 } << 907 Meson[AbsIDquark-1][ProdQ-1][StateQ]);
>> 908 LeftHadronMass=LeftHadron->GetPDGMass();
>> 909
>> 910 //G4cout<<"Meson "<<LeftHadron->GetParticleName()<<G4endl;
>> 911
>> 912 G4int StateDiQ=0;
>> 913 do // while(Baryon[Di_q1-1][Di_q2-1][ProdQ-1][StateDiQ]<>0);
>> 914 {
>> 915 //G4cout<<G4endl<<"Di_q1 Di_q2 ProdQ StateDiQ "<<BaryonWeight[Di_q1-1][Di_q2-1][ProdQ-1][StateDiQ]<<" "<<Di_q1-1<<" "<<Di_q2-1<<" "<<ProdQ-1<<" "<<StateDiQ<<G4endl;
>> 916 RightHadron=G4ParticleTable::GetParticleTable()->FindParticle(SignDiQ*
>> 917 Baryon[Di_q1-1][Di_q2-1][ProdQ-1][StateDiQ]);
>> 918 RightHadronMass=RightHadron->GetPDGMass();
>> 919
>> 920 //G4cout<<"Baryon "<<RightHadron->GetParticleName()<<G4endl;
>> 921
>> 922 //G4cout<<"StringMass LeftHadronMass RightHadronMass "<<StringMass<<" "<<LeftHadronMass<<" "<< RightHadronMass<<G4endl;
>> 923
>> 924 if(StringMass >= LeftHadronMass + RightHadronMass)
>> 925 {
>> 926 G4double FS_Psqr=lambda(StringMassSqr,sqr(LeftHadronMass),
>> 927 sqr(RightHadronMass));
>> 928 FS_Weight[NumberOf_FS]=std::sqrt(FS_Psqr)*
>> 929 MesonWeight[AbsIDquark-1][ProdQ-1][StateQ]*
>> 930 BaryonWeight[Di_q1-1][Di_q2-1][ProdQ-1][StateDiQ]*
>> 931 Prob_QQbar[ProdQ-1];
>> 932
>> 933 FS_LeftHadron[NumberOf_FS] = LeftHadron;
>> 934 FS_RightHadron[NumberOf_FS]= RightHadron;
>> 935
>> 936 //G4cout<<"State "<<NumberOf_FS<<" "<<std::sqrt(FS_Psqr/4./StringMassSqr)<<" "<<std::sqrt(FS_Psqr)<<" "<<FS_Weight[NumberOf_FS]<<G4endl;
>> 937 //G4cout<<"++++++++++++++++++++++++++++++++"<<G4endl;
>> 938 NumberOf_FS++;
>> 939
>> 940 if(NumberOf_FS > 24)
>> 941 {G4int Uzhi; G4cout<<"QQbar string #_FS "<<NumberOf_FS<<G4endl; G4cin>>Uzhi;}
>> 942 } // End of if(StringMass >= LeftHadronMass + RightHadronMass)
>> 943
>> 944 StateDiQ++;
>> 945 //G4cout<<Baryon[Di_q1-1][Di_q2-1][ProdQ-1][StateDiQ]<<" "<<Di_q1-1<<" "<<Di_q2-1<<" "<<ProdQ-1<<" "<<StateDiQ<<G4endl;
>> 946 } while(Baryon[Di_q1-1][Di_q2-1][ProdQ-1][StateDiQ]!=0);
>> 947
>> 948 StateQ++;
>> 949 } while(Meson[AbsIDquark-1][ProdQ-1][StateQ]!=0);
>> 950 } // End of for(G4int ProdQ=1; ProdQ < 4; ProdQ++)
>> 951 }
>> 952 //====================================
>> 953
>> 954 if(NumberOf_FS == 0) return false;
>> 955 //G4cout<<"NumberOf_FS "<<NumberOf_FS<<G4endl;
>> 956 G4double SumWeights=0.;
>> 957 for(G4int i=0; i<NumberOf_FS; i++) {SumWeights+=FS_Weight[i];}
>> 958
>> 959 G4double ksi=G4UniformRand();
>> 960 G4double Sum=0.;
>> 961 G4int SampledState(0);
>> 962
>> 963 for(G4int i=0; i<NumberOf_FS; i++)
>> 964 {
>> 965 Sum+=(FS_Weight[i]/SumWeights);
>> 966 SampledState=i;
>> 967 if(Sum >= ksi) break;
>> 968 }
>> 969
>> 970 LeftHadron =FS_LeftHadron[SampledState];
>> 971 RightHadron=FS_RightHadron[SampledState];
>> 972
>> 973 //G4cout<<"Selected "<<SampledState<<" "<<LeftHadron->GetParticleName()<<" "<<RightHadron->GetParticleName()<<G4endl;
>> 974
>> 975 } // End of if(!string->FourQuarkString())
>> 976
>> 977 G4LorentzVector LeftMom, RightMom;
>> 978 G4ThreeVector Pos;
>> 979
>> 980 Sample4Momentum(&LeftMom, LeftHadron->GetPDGMass(),
>> 981 &RightMom, RightHadron->GetPDGMass(),
>> 982 StringMass);
1052 983
1053 if (AbsIDquark == ProdQ) SignQ << 984 LeftMom.boost(ClusterVel);
>> 985 RightMom.boost(ClusterVel);
1054 986
1055 G4int StateQ=0; << 987 LeftVector->push_back(new G4KineticTrack(LeftHadron, 0, Pos, LeftMom));
1056 const G4int maxNumberOfLoops << 988 RightVector->push_back(new G4KineticTrack(RightHadron, 0, Pos, RightMom));
1057 G4int loopCounter = 0; << 989
1058 do // while(Meson[AbsIDquark-1][ProdQ-1] << 990 return true;
1059 { <<
1060 LeftHadron=G4ParticleTable::GetParticle <<
1061 Meson[AbsIDquark-1][ProdQ-1][St <<
1062 if (LeftHadron == NULL) continue; <<
1063 G4double LeftHadronMass=LeftHadron->Get <<
1064 <<
1065 G4int StateDiQ=0; <<
1066 const G4int maxNumber <<
1067 G4int internalLoopCou <<
1068 do // while(Baryon[Di_q1-1][Di_q2-1][Pr <<
1069 { <<
1070 RightHadron=G4ParticleTable::GetParti <<
1071 Baryon[Di_q1-1][Di_q2-1][Prod <<
1072 if (RightHadron == NULL) continue; <<
1073 G4double RightHadronMass=RightHadron- <<
1074 <<
1075 if (StringMass > LeftHadronMass + Rig <<
1076 { <<
1077 if ( <<
1078 G4E <<
1079 ed <<
1080 G4E <<
1081 <<
1082 Num <<
1083 } <<
1084 <<
1085 G4double FS_Psqr=lambda(StringMassS <<
1086 sqr(RightHadronMass)); <<
1087 FS_Weight[NumberOf_FS]=std::sqrt(FS <<
1088 MesonWeight[AbsIDquark-1 <<
1089 BaryonWeight[Di_q1-1][Di <<
1090 Prob_QQbar[ProdQ-1]; <<
1091 <<
1092 FS_LeftHadron[NumberOf_FS] = LeftHa <<
1093 FS_RightHadron[NumberOf_FS]= RightH <<
1094 <<
1095 NumberOf_FS++; <<
1096 } // End of if (StringMass > LeftHadr <<
1097 <<
1098 StateDiQ++; <<
1099 <<
1100 } while( (Baryon[Di_q1-1][Di_q2-1][Prod <<
1101 ++internalLo <<
1102 if ( internalLoopCoun <<
1103 return false; <<
1104 } <<
1105 <<
1106 StateQ++; <<
1107 } while( (Meson[AbsIDquark-1][ProdQ-1][St <<
1108 ++loopCounter < maxN <<
1109 <<
1110 if ( loopCounter >= maxNumb <<
1111 return false; <<
1112 } <<
1113 } <<
1114 991
1115 return true; <<
1116 } 992 }
1117 993
1118 //------------------------------------------- << 994 //----------------------------------------------------------------------------------------------------------
1119 995
1120 G4bool G4LundStringFragmentation::Quark_AntiQ << 996 void G4LundStringFragmentation::Sample4Momentum(G4LorentzVector* Mom, G4double Mass, G4LorentzVector* AntiMom, G4double AntiMass, G4double InitialMass)
1121 << 997 {
1122 << 998 // ------ Sampling of momenta of 2 last produced hadrons --------------------
1123 { << 999 G4ThreeVector Pt;
1124 G4double StringMass = string->Mass(); << 1000 G4double MassMt2, AntiMassMt2;
1125 G4double StringMassSqr= sqr(StringMass); << 1001 G4double AvailablePz, AvailablePz2;
>> 1002
>> 1003 //G4cout<<"Masses "<<InitialMass<<" "<<Mass<<" "<<AntiMass<<G4endl;
>> 1004 //
>> 1005
>> 1006 if((Mass > 930. || AntiMass > 930.)) //If there is a baryon
>> 1007 {
>> 1008 // ----------------- Isotropic decay ------------------------------------
>> 1009 G4double r_val = sqr(InitialMass*InitialMass - Mass*Mass - AntiMass*AntiMass) -
>> 1010 sqr(2.*Mass*AntiMass);
>> 1011 G4double Pabs = (r_val > 0.)? std::sqrt(r_val)/(2.*InitialMass) : 0;
>> 1012 //G4cout<<"P for isotr decay "<<Pabs<<G4endl;
>> 1013
>> 1014 //... sample unit vector
>> 1015 G4double pz =1. - 2.*G4UniformRand();
>> 1016 G4double st = std::sqrt(1. - pz * pz)*Pabs;
>> 1017 G4double phi = 2.*pi*G4UniformRand();
>> 1018 G4double px = st*std::cos(phi);
>> 1019 G4double py = st*std::sin(phi);
>> 1020 pz *= Pabs;
>> 1021
>> 1022 Mom->setPx(px); Mom->setPy(py); Mom->setPz(pz);
>> 1023 Mom->setE(std::sqrt(Pabs*Pabs + Mass*Mass));
1126 1024
1127 G4ParticleDefinition * Quark; << 1025 AntiMom->setPx(-px); AntiMom->setPy(-py); AntiMom->setPz(-pz);
1128 G4ParticleDefinition * Anti_Quark; << 1026 AntiMom->setE (std::sqrt(Pabs*Pabs + AntiMass*AntiMass));
>> 1027 //G4int Uzhi; G4cin>>Uzhi;
>> 1028 }
>> 1029 else
>> 1030 //
>> 1031 {
>> 1032 do
>> 1033 {
>> 1034 // GF 22-May-09, limit sampled pt to allowed range
>> 1035
>> 1036 G4double termD = InitialMass*InitialMass -Mass*Mass - AntiMass*AntiMass;
>> 1037 G4double termab = 4*sqr(Mass*AntiMass);
>> 1038 G4double termN = 2*termD + 4*Mass*Mass + 4*AntiMass*AntiMass;
>> 1039 G4double pt2max=(termD*termD - termab )/ termN ;
>> 1040 //G4cout<<"Anis "<<pt2max<<" "<<(termD*termD-termab)/(4.*InitialMass*InitialMass)<<G4endl;
>> 1041
>> 1042 Pt=SampleQuarkPt(std::sqrt(pt2max)); Pt.setZ(0); G4double Pt2=Pt.mag2();
>> 1043 //G4cout<<"Sampl pt2 "<<Pt2<<G4endl;
>> 1044 MassMt2 = Mass * Mass + Pt2;
>> 1045 AntiMassMt2= AntiMass * AntiMass + Pt2;
>> 1046
>> 1047 AvailablePz2= sqr(InitialMass*InitialMass - MassMt2 - AntiMassMt2) -
>> 1048 4.*MassMt2*AntiMassMt2;
>> 1049 }
>> 1050 while(AvailablePz2 < 0.); // GF will occur only for numerical precision problem with limit in sampled pt
>> 1051
>> 1052 AvailablePz2 /=(4.*InitialMass*InitialMass);
>> 1053 AvailablePz = std::sqrt(AvailablePz2);
>> 1054
>> 1055 G4double Px=Pt.getX();
>> 1056 G4double Py=Pt.getY();
1129 1057
1130 if (string->GetLeftParton()->GetPDGEncoding << 1058 Mom->setPx(Px); Mom->setPy(Py); Mom->setPz(AvailablePz);
1131 { << 1059 Mom->setE(std::sqrt(MassMt2+AvailablePz2));
1132 Quark =string->GetLeftParton(); <<
1133 Anti_Quark=string->GetRightParton(); <<
1134 } else <<
1135 { <<
1136 Quark =string->GetRightParton(); <<
1137 Anti_Quark=string->GetLeftParton(); <<
1138 } <<
1139 1060
1140 G4int IDquark =Quark->GetPDGEncodin << 1061 AntiMom->setPx(-Px); AntiMom->setPy(-Py); AntiMom->setPz(-AvailablePz);
1141 G4int AbsIDquark =std::abs(IDquark); << 1062 AntiMom->setE (std::sqrt(AntiMassMt2+AvailablePz2));
1142 G4int QuarkCharge =Qcharge[IDquar << 1063 }
1143 << 1064 }
1144 G4int IDanti_quark =Anti_Quark->GetPDGEn <<
1145 G4int AbsIDanti_quark =std::abs(IDanti_quar <<
1146 G4int AntiQuarkCharge =-Qcharge[AbsID <<
1147 <<
1148 G4int LeftHadronCharge(0), RightHadro <<
1149 <<
1150 //G4cout<<"Q Qbar "<<IDquark<<" "<<ID <<
1151 <<
1152 NumberOf_FS=0; <<
1153 for (G4int ProdQ=1; ProdQ < 4; ProdQ++) // <<
1154 { // <<
1155 //G4cout<<"NumberOf_FS ProdQ <<
1156 LeftHadronCharge = QuarkCharge - Qcharge[ <<
1157 G4int SignQ = LeftHadronCharge/3; if (Sig <<
1158 <<
1159 if ((IDquark == 1) && (ProdQ == 3)) SignQ <<
1160 if ((IDquark == 3) && (ProdQ == 1)) SignQ <<
1161 if ((IDquark == 4) && (ProdQ <<
1162 if ((IDquark == 5) && (ProdQ <<
1163 if ((IDquark == 5) && (ProdQ <<
1164 <<
1165 RightHadronCharge = AntiQuark <<
1166 G4int SignAQ = RightHadronCharge/3; if (S <<
1167 1065
1168 if ((IDanti_quark ==-1) && (ProdQ == 3)) << 1066 //-----------------------------------------------------------------------------
1169 if ((IDanti_quark ==-3) && (ProdQ == 1)) <<
1170 if ((IDanti_quark ==-4) && (ProdQ == 2)) <<
1171 if ((IDanti_quark ==-5) && (ProdQ == 1)) <<
1172 if ((IDanti_quark ==-5) && (ProdQ == 3)) <<
1173 <<
1174 //G4cout<<"ProQ signs "<<Prod <<
1175 <<
1176 G4int StateQ=0; <<
1177 const G4int maxNumberOfLoops <<
1178 G4int loopCounter = 0; <<
1179 do <<
1180 { <<
1181 //G4cout<<"[AbsIDquar <<
1182 //<<ProdQ-1<<" "<<Sta <<
1183 LeftHadron=G4ParticleTable::GetParticle <<
1184 Meson[AbsIDquark-1][ProdQ- <<
1185 //G4cout<<"LeftHadron <<
1186 if (LeftHadron == NULL) { StateQ++; con <<
1187 //G4cout<<"LeftHadron <<
1188 G4double LeftHadronMass=LeftHadron->Get <<
1189 <<
1190 G4int StateAQ=0; <<
1191 const G4int maxNumber <<
1192 G4int internalLoopCou <<
1193 do <<
1194 { <<
1195 //G4cout<<" <<
1196 // <<Pro <<
1197 RightHadron=G4ParticleTable::GetParti <<
1198 Meson[AbsIDanti_quark-1][Prod <<
1199 //G4cout<<"Ri <<
1200 if(RightHadron == NULL) { StateAQ++; <<
1201 //G4cout<<"Ri <<
1202 G4double RightHadronMass=RightHadron- <<
1203 <<
1204 if (StringMass > LeftHadronMass + Rig <<
1205 { <<
1206 if ( <<
1207 G4E <<
1208 ed <<
1209 G4E <<
1210 <<
1211 Num <<
1212 } <<
1213 <<
1214 G4dou <<
1215 sqr(RightHadronMass)); <<
1216 //FS_Psqr=1.; <<
1217 FS_Weight[NumberOf_FS]=std::sqrt(FS <<
1218 MesonWeight[AbsIDquark-1 <<
1219 MesonWeight[AbsIDanti_qu <<
1220 Prob_QQbar[ProdQ-1]; <<
1221 if (string->GetLeftParton()->GetPDG <<
1222 { <<
1223 FS_LeftHadron[NumberOf_FS] = Righ <<
1224 FS_RightHadron[NumberOf_FS]= Left <<
1225 } else <<
1226 { <<
1227 FS_LeftHadron[NumberOf_FS] = Left <<
1228 FS_RightHadron[NumberOf_FS]= Righ <<
1229 } <<
1230 <<
1231 NumberOf_FS++; <<
1232 <<
1233 } <<
1234 <<
1235 StateAQ++; <<
1236 //G4cout<<" <<
1237 // <<Mes <<
1238 } while ( (Meson[AbsIDanti_quark-1][Pro <<
1239 ++internalL <<
1240 if ( internalLoopCo <<
1241 return false; <<
1242 } <<
1243 <<
1244 StateQ++; <<
1245 //G4cout<<"StateQ Mes <<
1246 // <<Meson[AbsID <<
1247 <<
1248 } while ( (Meson[AbsIDquark-1][ProdQ-1][S <<
1249 ++loopCounter < maxN <<
1250 if ( loopCounter >= maxNumb <<
1251 return false; <<
1252 } <<
1253 } // End of for (G4int ProdQ=1; ProdQ < 4; <<
1254 1067
1255 return true; << 1068 G4LorentzVector * G4LundStringFragmentation::SplitEandP(G4ParticleDefinition * pHadron,
1256 } << 1069 G4FragmentingString * string, G4FragmentingString * newString)
>> 1070 {
>> 1071 //G4cout<<"Start SplitEandP "<<G4endl;
>> 1072 G4LorentzVector String4Momentum=string->Get4Momentum();
>> 1073 G4double StringMT2=string->Get4Momentum().mt2();
1257 1074
1258 //------------------------------------------- << 1075 G4double HadronMass = pHadron->GetPDGMass();
1259 1076
1260 G4int G4LundStringFragmentation::SampleState( << 1077 SetMinimalStringMass(newString);
1261 { << 1078 //G4cout<<"HadM MinimalStringMassLeft StringM "<<HadronMass<<" "<<MinimalStringMass<<" "<<String4Momentum.mag()<<G4endl;
1262 if ( NumberOf_FS > 349 ) { << 1079
1263 G4ExceptionDescription ed; << 1080 if(HadronMass + MinimalStringMass > String4Momentum.mag()) {return 0;}// have to start all over!
1264 ed << " NumberOf_FS exceeds its lim << 1081 String4Momentum.setPz(0.);
1265 G4Exception( "G4LundStringFragmenta << 1082 G4ThreeVector StringPt=String4Momentum.vect();
1266 NumberOf_FS = 349; <<
1267 } <<
1268 1083
1269 G4double SumWeights=0.; << 1084 // calculate and assign hadron transverse momentum component HadronPx and HadronPy
1270 for (G4int i=0; i<NumberOf_FS; i++) {SumWei << 1085 G4ThreeVector thePt;
>> 1086 thePt=SampleQuarkPt();
1271 1087
1272 G4double ksi=G4UniformRand(); << 1088 G4ThreeVector HadronPt = thePt +string->DecayPt();
1273 G4double Sum=0.; << 1089 HadronPt.setZ(0);
1274 G4int indexPosition = 0; <<
1275 1090
1276 for (G4int i=0; i<NumberOf_FS; i++) << 1091 G4ThreeVector RemSysPt = StringPt - HadronPt;
1277 { <<
1278 Sum+=(FS_Weight[i]/SumWeights); <<
1279 indexPosition=i; <<
1280 if (Sum >= ksi) break; <<
1281 } <<
1282 return indexPosition; <<
1283 } <<
1284 1092
1285 //------------------------------------------- << 1093 //... sample z to define hadron longitudinal momentum and energy
>> 1094 //... but first check the available phase space
1286 1095
1287 void G4LundStringFragmentation::Sample4Moment << 1096 G4double HadronMassT2 = sqr(HadronMass) + HadronPt.mag2();
1288 << 1097 G4double ResidualMassT2=sqr(MinimalStringMass) + RemSysPt.mag2();
1289 <<
1290 { <<
1291 // ------ Sampling of momenta of 2 last pro <<
1292 G4ThreeVector Pt; <<
1293 G4double MassMt, AntiMassMt; <<
1294 G4double AvailablePz, AvailablePz2; <<
1295 #ifdef debug_LUNDfragmentation <<
1296 G4cout<<"Sampling of momenta of 2 las <<
1297 G4cout<<"Init Mass "<<InitialMass<<" <<
1298 #endif <<
1299 <<
1300 G4double r_val = sqr(InitialMass*InitialMas <<
1301 sqr(2.*Mass*AntiMass); <<
1302 G4double Pabs = (r_val > 0.)? std::sqrt(r_v <<
1303 <<
1304 const G4int maxNumberOfLoops = 1000; <<
1305 G4double SigmaQTw = SigmaQT; <<
1306 if ( Mass > 930. || AntiMass > 930. ) { <<
1307 SigmaQT *= ( 1.0 - 0.55*sqr( (Mass+AntiMa <<
1308 } <<
1309 if ( Mass < 930. && AntiMass < 930. ) <<
1310 if ( ( Mass < 930. && AntiMass > 930. <<
1311 ( Mass > 930. && AntiMass < 930. ) ) { <<
1312 //SigmaQT = -1.; // isotropical de <<
1313 } <<
1314 if ( Mass > 930. && AntiMass > 930. ) <<
1315 SigmaQT *= ( 1.0 - 0.55*sqr( (Mass+ <<
1316 } <<
1317 1098
1318 G4int loopCounter = 0; << 1099 G4double Pz2 = (sqr(StringMT2 - HadronMassT2 - ResidualMassT2) -
1319 do << 1100 4*HadronMassT2 * ResidualMassT2)/4./StringMT2;
1320 { << 1101 //G4cout<<"Pz2 "<<Pz2<<G4endl;
1321 Pt=SampleQuarkPt(Pabs); Pt.setZ(0); G4dou << 1102 if(Pz2 < 0 ) {return 0;} // have to start all over!
1322 MassMt = std::sqrt( Mass * Mass <<
1323 AntiMassMt= std::sqrt(AntiMass * AntiMass <<
1324 } <<
1325 while ( (InitialMass < MassMt + AntiMassMt) <<
1326 1103
1327 SigmaQT = SigmaQTw; << 1104 //... then compute allowed z region z_min <= z <= z_max
>> 1105
>> 1106 G4double Pz = std::sqrt(Pz2);
>> 1107 G4double zMin = (std::sqrt(HadronMassT2+Pz2) - Pz)/std::sqrt(StringMT2);
>> 1108 G4double zMax = (std::sqrt(HadronMassT2+Pz2) + Pz)/std::sqrt(StringMT2);
1328 1109
1329 AvailablePz2= sqr(InitialMass*InitialMass - << 1110 //G4cout<<"if (zMin >= zMax) return 0 "<<zMin<<" "<<zMax<<G4endl;
1330 4.*sqr(MassMt*AntiMassMt); << 1111 if (zMin >= zMax) return 0; // have to start all over!
>> 1112
>> 1113 G4double z = GetLightConeZ(zMin, zMax,
>> 1114 string->GetDecayParton()->GetPDGEncoding(), pHadron,
>> 1115 HadronPt.x(), HadronPt.y());
>> 1116 //G4cout<<"z "<<z<<G4endl;
>> 1117 //... now compute hadron longitudinal momentum and energy
>> 1118 // longitudinal hadron momentum component HadronPz
>> 1119
>> 1120 HadronPt.setZ(0.5* string->GetDecayDirection() *
>> 1121 (z * string->LightConeDecay() -
>> 1122 HadronMassT2/(z * string->LightConeDecay())));
>> 1123
>> 1124 G4double HadronE = 0.5* (z * string->LightConeDecay() +
>> 1125 HadronMassT2/(z * string->LightConeDecay()));
>> 1126
>> 1127 G4LorentzVector * a4Momentum= new G4LorentzVector(HadronPt,HadronE);
>> 1128 //G4cout<<"Out SplitEandP "<<G4endl;
>> 1129 return a4Momentum;
>> 1130 }
1331 1131
1332 AvailablePz2 /=(4.*InitialMass*InitialMass) << 1132 //-----------------------------------------------------------------------------------------
1333 AvailablePz = std::sqrt(AvailablePz2); << 1133 G4double G4LundStringFragmentation::GetLightConeZ(G4double zmin, G4double zmax,
>> 1134 G4int PDGEncodingOfDecayParton,
>> 1135 G4ParticleDefinition* pHadron,
>> 1136 G4double Px, G4double Py)
>> 1137 {
>> 1138 G4double alund;
1334 1139
1335 G4double Px=Pt.getX(); << 1140 // If blund get restored, you MUST adapt the calculation of zOfMaxyf.
1336 G4double Py=Pt.getY(); << 1141 // const G4double blund = 1;
1337 1142
1338 Mom->setPx(Px); Mom->setPy(Py); Mom->setPz( << 1143 G4double z, yf;
1339 Mom->setE(std::sqrt(sqr(MassMt)+AvailablePz << 1144 G4double Mass = pHadron->GetPDGMass();
>> 1145 // G4int HadronEncoding=pHadron->GetPDGEncoding();
>> 1146
>> 1147 G4double Mt2 = Px*Px + Py*Py + Mass*Mass;
1340 1148
1341 AntiMom->setPx(-Px); AntiMom->setPy(-Py); A << 1149 if(std::abs(PDGEncodingOfDecayParton) < 1000)
1342 AntiMom->setE (std::sqrt(sqr(AntiMassMt)+Av << 1150 {
>> 1151 // ---------------- Quark fragmentation ----------------------
>> 1152 alund=0.35/GeV/GeV; // Instead of 0.7 because kinks are not considered
>> 1153 G4double zOfMaxyf=alund*Mt2/(alund*Mt2 + 1.);
>> 1154 G4double maxYf=(1-zOfMaxyf)/zOfMaxyf * std::exp(-alund*Mt2/zOfMaxyf);
>> 1155
>> 1156 do
>> 1157 {
>> 1158 z = zmin + G4UniformRand()*(zmax-zmin);
>> 1159 // yf = std::pow(1. - z, blund)/z*std::exp(-alund*Mt2/z);
>> 1160 yf = (1-z)/z * std::exp(-alund*Mt2/z);
>> 1161 }
>> 1162 while (G4UniformRand()*maxYf > yf);
>> 1163 }
>> 1164 else
>> 1165 {
>> 1166 // ---------------- Di-quark fragmentation ----------------------
>> 1167 alund=0.7/GeV/GeV; // 0.7 2.0
>> 1168 G4double zOfMaxyf=alund*Mt2/(alund*Mt2 + 1.);
>> 1169 G4double maxYf=(1-zOfMaxyf)/zOfMaxyf * std::exp(-alund*Mt2/zOfMaxyf);
>> 1170 do
>> 1171 {
>> 1172 z = zmin + G4UniformRand()*(zmax-zmin);
>> 1173 // yf = std::pow(1. - z, blund)/z*std::exp(-alund*Mt2/z);
>> 1174 yf = (1-z)/z * std::exp(-alund*Mt2/z);
>> 1175 }
>> 1176 while (G4UniformRand()*maxYf > yf);
>> 1177 };
1343 1178
1344 #ifdef debug_LUNDfragmentation << 1179 return z;
1345 G4cout<<"Fmass Mom "<<Mom->getX()<<" << 1180 }
1346 G4cout<<"Smass Mom "<<AntiMom->getX() <<
1347 <<" "<<AntiMom->getT()<<G4endl; <<
1348 #endif <<
1349 } <<
1350 1181
1351 //------------------------------------------- 1182 //------------------------------------------------------------------------
1352 << 1183 G4double G4LundStringFragmentation::lambda(G4double s, G4double m1_Sqr, G4double m2_Sqr)
1353 G4double G4LundStringFragmentation::lambda(G4 <<
1354 { 1184 {
1355 G4double lam = sqr(S - m1_Sqr - m2_Sqr) - 4 << 1185 G4double lam = sqr(s - m1_Sqr - m2_Sqr) - 4.*m1_Sqr*m2_Sqr;
1356 return lam; << 1186 return lam;
1357 } 1187 }
1358 <<
1359 // ------------------------------------------ <<
1360 G4LundStringFragmentation::~G4LundStringFragm <<
1361 {} <<
1362 1188
1363 1189