| 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 106967 2017-10-31 08:41:49Z gcosmo $
>> 28 // GEANT4 tag $Name: $ 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" 36 #include "G4PhysicalConstants.hh"
35 #include "G4SystemOfUnits.hh" 37 #include "G4SystemOfUnits.hh"
36 #include "Randomize.hh" 38 #include "Randomize.hh"
37 #include "G4FragmentingString.hh" 39 #include "G4FragmentingString.hh"
38 #include "G4DiQuarks.hh" 40 #include "G4DiQuarks.hh"
39 #include "G4Quarks.hh" 41 #include "G4Quarks.hh"
40 #include "G4HadronicParameters.hh" << 42
41 #include "G4Exp.hh" 43 #include "G4Exp.hh"
42 #include "G4Pow.hh" 44 #include "G4Pow.hh"
43 45
44 //#define debug_LUNDfragmentation 46 //#define debug_LUNDfragmentation
45 47
46 // Class G4LundStringFragmentation 48 // Class G4LundStringFragmentation
47 //******************************************** 49 //*************************************************************************************
48 50
49 G4LundStringFragmentation::G4LundStringFragmen 51 G4LundStringFragmentation::G4LundStringFragmentation()
50 : G4VLongitudinalStringDecay("LundStringFrag <<
51 { 52 {
52 SetMassCut(210.*MeV); // Mpi + Delta << 53 // ------ For estimation of a minimal string mass ---------------
53 // For ProduceOneH << 54 Mass_of_light_quark =140.*MeV;
54 // that no one pi- << 55 Mass_of_heavy_quark =500.*MeV;
55 SigmaQT = 0.435 * GeV; << 56 Mass_of_string_junction=720.*MeV;
56 Tmt = 190.0 * MeV; << 57 // ------ An estimated minimal string mass ----------------------
57 << 58 MinimalStringMass = 0.;
58 SetStringTensionParameter(1.*GeV/fermi); << 59 MinimalStringMass2 = 0.;
59 SetDiquarkBreakProbability(0.3); << 60 // ------ Minimal invariant mass used at a string fragmentation -
60 << 61 WminLUND = 0.45*GeV; //0.23*GeV; // Uzhi 0.7 -> 0.23 3.8.10 //0.8 1.5
61 SetStrangenessSuppression((1.0 - 0.12)/2.0 << 62 // ------ Smooth parameter used at a string fragmentation for ---
62 SetDiquarkSuppression(0.07); << 63 // ------ smearing sharp mass cut-off ---------------------------
63 << 64 SmoothParam = 0.2;
64 // Check if charmed and bottom hadrons are << 65
65 // set the non-zero probabilities for c-cb << 66 SetStringTensionParameter(1.);
66 // else set them to 0.0. If these probabil << 67 SetDiquarkBreakProbability(0.05);
67 // hadrons can't/can be created during the << 68 SetStrangenessSuppression(0.46); //(0.447);
68 // (i.e. not heavy) projectile hadron nucl << 69 SetDiquarkSuppression(0.05);
69 if ( G4HadronicParameters::Instance()->Ena << 70
70 SetProbCCbar(0.0002); // According to O << 71 // For treating of small string decays
71 SetProbBBbar(5.0e-5); // According to O << 72 for(G4int i=0; i<3; i++)
72 } else { << 73 { for(G4int j=0; j<3; j++)
73 SetProbCCbar(0.0); << 74 { for(G4int k=0; k<6; k++)
74 SetProbBBbar(0.0); << 75 {
>> 76 Meson[i][j][k]=0; MesonWeight[i][j][k]=0.;
>> 77 }
75 } 78 }
76 << 79 }
77 SetMinMasses(); // For treating of small << 80 //--------------------------
78 } << 81 Meson[0][0][0]=111; // dbar-d Pi0
>> 82 MesonWeight[0][0][0]=(1.-pspin_meson)*(1.-scalarMesonMix[0]);
79 83
80 //-------------------------------------------- << 84 Meson[0][0][1]=221; // dbar-d Eta
>> 85 MesonWeight[0][0][1]=(1.-pspin_meson)*(scalarMesonMix[0]-scalarMesonMix[1]);
81 86
82 G4KineticTrackVector* G4LundStringFragmentatio << 87 Meson[0][0][2]=331; // dbar-d EtaPrime
83 { << 88 MesonWeight[0][0][2]=(1.-pspin_meson)*(scalarMesonMix[1]);
84 // Can no longer modify Parameters for Fragm <<
85 89
86 PastInitPhase=true; << 90 Meson[0][0][3]=113; // dbar-d Rho0
>> 91 MesonWeight[0][0][3]=pspin_meson*(1.-vectorMesonMix[0]);
87 92
88 G4FragmentingString aString(theString); << 93 Meson[0][0][4]=223; // dbar-d Omega
89 SetMinimalStringMass(&aString); << 94 MesonWeight[0][0][4]=pspin_meson*(vectorMesonMix[0]);
>> 95 //--------------------------
90 96
91 #ifdef debug_LUNDfragmentation << 97 Meson[0][1][0]=211; // dbar-u Pi+
92 G4cout<<G4endl<<"LUND StringFragmentat << 98 MesonWeight[0][1][0]=(1.-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 99
105 G4KineticTrackVector * LeftVector(0); << 100 Meson[0][1][1]=213; // dbar-u Rho+
>> 101 MesonWeight[0][1][1]=pspin_meson;
>> 102 //--------------------------
106 103
107 if (!aString.IsAFourQuarkString() && !IsItFr << 104 Meson[0][2][0]=311; // dbar-s K0bar
108 { << 105 MesonWeight[0][2][0]=(1.-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 <<
120 if ( LeftVector ) <<
121 { <<
122 if ( LeftVector->size() > 0) <<
123 { <<
124 LeftVector->operator[](0)->SetForm <<
125 LeftVector->operator[](0)->SetPosi <<
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[0][2][1]=313; // dbar-s K*0bar
138 G4cout<<"The string will be fragmented << 108 MesonWeight[0][2][1]=pspin_meson;
139 #endif << 109 //--------------------------
140 110
141 // The string can fragment. At least two par << 111 Meson[1][0][0]=211; // ubar-d Pi-
142 LeftVector =new G4KineticTrackVec << 112 MesonWeight[1][0][0]=(1.-pspin_meson);
143 G4KineticTrackVector * RightVector=new G4Kin <<
144 113
145 G4bool success = Loop_toFragmentString << 114 Meson[1][0][1]=213; // ubar-d Rho-
>> 115 MesonWeight[1][0][1]=pspin_meson;
>> 116 //--------------------------
146 117
147 if ( ! success ) << 118 Meson[1][1][0]=111; // ubar-u Pi0
148 { << 119 MesonWeight[1][1][0]=(1.-pspin_meson)*(1.-scalarMesonMix[0]);
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 120
156 // Join Left- and RightVector into LeftVecto << 121 Meson[1][1][1]=221; // ubar-u Eta
157 while (!RightVector->empty()) << 122 MesonWeight[1][1][1]=(1.-pspin_meson)*(scalarMesonMix[0]-scalarMesonMix[1]);
158 { <<
159 LeftVector->push_back(RightVector->back()) <<
160 RightVector->erase(RightVector->end()-1); <<
161 } <<
162 delete RightVector; <<
163 123
164 return LeftVector; << 124 Meson[1][1][2]=331; // ubar-u EtaPrime
165 } << 125 MesonWeight[1][1][2]=(1.-pspin_meson)*(scalarMesonMix[1]);
166 126
167 //-------------------------------------------- << 127 Meson[1][1][3]=113; // ubar-u Rho0
>> 128 MesonWeight[1][1][3]=pspin_meson*(1.-vectorMesonMix[0]);
168 129
169 G4bool G4LundStringFragmentation::IsItFragment << 130 Meson[1][1][4]=223; // ubar-u Omega
170 { << 131 //MesonWeight[1][1][4]=pspin_meson*(scalarMesonMix[0]);
171 SetMinimalStringMass(string); << 132 MesonWeight[1][1][4]=pspin_meson*(vectorMesonMix[0]); // Uzhi 2015 scalar -> vector
172 //G4cout<<"MinM StrM "<<MinimalStringM << 133 //--------------------------
173 134
174 return std::abs(MinimalStringMass) < string- << 135 Meson[1][2][0]=321; // ubar-s K-
>> 136 MesonWeight[1][2][0]=(1.-pspin_meson);
175 137
176 //MinimalStringMass is negative and la << 138 Meson[1][2][1]=323; // ubar-s K*-bar -
177 } << 139 MesonWeight[1][2][1]=pspin_meson;
>> 140 //--------------------------
178 141
179 //-------------------------------------------- << 142 Meson[2][0][0]=311; // sbar-d K0
>> 143 MesonWeight[2][0][0]=(1.-pspin_meson);
180 144
181 G4bool G4LundStringFragmentation::Loop_toFragm << 145 Meson[2][0][1]=313; // sbar-d K*0
182 << 146 MesonWeight[2][0][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][1][0]=321; // sbar-u K+
>> 150 MesonWeight[2][1][0]=(1.-pspin_meson);
194 151
195 G4bool final_success=false; << 152 Meson[2][1][1]=323; // sbar-u K*+
196 G4bool inner_success=true; << 153 MesonWeight[2][1][1]=pspin_meson;
>> 154 //--------------------------
197 155
198 G4int attempt=0; << 156 Meson[2][2][0]=221; // sbar-s Eta
>> 157 MesonWeight[2][2][0]=(1.-pspin_meson)*(1.-scalarMesonMix[5]);
199 158
200 while ( ! final_success && attempt++ < Strin << 159 Meson[2][2][1]=331; // sbar-s EtaPrime
201 { // If the string fragmentation does << 160 MesonWeight[2][2][1]=(1.-pspin_meson)*(1.-scalarMesonMix[5]);
202 // repeat the fragmentation. <<
203 161
204 G4FragmentingString *currentSt << 162 Meson[2][2][3]=333; // sbar-s EtaPrime
205 toCmsI = currentString->Transf << 163 MesonWeight[2][2][3]=pspin_meson*(vectorMesonMix[5]);
206 toObserverFrameI = toCmsI.inve << 164 //--------------------------
207 165
208 G4LorentzRotation toCms, toObserverFrame; << 166 for(G4int i=0; i<3; i++)
>> 167 { for(G4int j=0; j<3; j++)
>> 168 { for(G4int k=0; k<3; k++)
>> 169 { for(G4int l=0; l<4; l++)
>> 170 { Baryon[i][j][k][l]=0; BaryonWeight[i][j][k][l]=0.;}
>> 171 }
>> 172 }
>> 173 }
209 174
210 //G4cout<<"Main loop start whilecounter "< << 175 G4double pspin_barion_in=pspin_barion;
>> 176 //pspin_barion=0.75;
>> 177 //---------------------------------------
211 178
212 // Cleaning up the previously produced had << 179 Baryon[0][0][0][0]=1114; // Delta-
213 std::for_each(LeftVector->begin(), LeftVec << 180 BaryonWeight[0][0][0][0]=1.;
214 LeftVector->clear(); << 181 //---------------------------------------
215 std::for_each(RightVector->begin(), RightV <<
216 RightVector->clear(); <<
217 182
218 // Main fragmentation loop until the strin << 183 Baryon[0][0][1][0]=2112; // neutron
219 inner_success=true; // set false on failu << 184 BaryonWeight[0][0][1][0]=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 <<
240 G4KineticTrack * Hadron=Splitup(currentS <<
241 <<
242 if ( Hadron != 0 ) // Store the hadron <<
243 { <<
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 <<
286 if ( inner_success && SplitLast(currentStr <<
287 { <<
288 final_success = true; <<
289 } <<
290 <<
291 delete currentString; <<
292 } // End of the loop where we try to fragme <<
293 <<
294 G4int sign = +1; <<
295 if ( theString.GetDirection() < 0 ) si <<
296 for ( unsigned int hadronI = 0; hadron <<
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 185
309 return final_success; << 186 Baryon[0][0][1][1]=2114; // Delta0
310 } << 187 BaryonWeight[0][0][1][1]=pspin_barion;
>> 188 //---------------------------------------
311 189
312 //-------------------------------------------- << 190 Baryon[0][0][2][0]=3112; // Sigma-
>> 191 BaryonWeight[0][0][2][0]=1.-pspin_barion;
313 192
314 G4bool G4LundStringFragmentation::StopFragment << 193 Baryon[0][0][2][1]=3114; // Sigma*-
315 { << 194 BaryonWeight[0][0][2][1]=pspin_barion;
316 SetMinimalStringMass(string); << 195 //---------------------------------------
317 196
318 if ( MinimalStringMass < 0.) return true; << 197 Baryon[0][1][0][0]=2112; // neutron
>> 198 BaryonWeight[0][1][0][0]=1.-pspin_barion;
319 199
320 if (string->IsAFourQuarkString()) << 200 Baryon[0][1][0][1]=2114; // Delta0
321 { << 201 BaryonWeight[0][1][0][1]=pspin_barion;
322 return G4UniformRand() < G4Exp(-0.0005*(st << 202 //---------------------------------------
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 203
340 //-------------------------------------------- << 204 Baryon[0][1][1][0]=2212; // proton
>> 205 BaryonWeight[0][1][1][0]=1.-pspin_barion;
341 206
342 G4KineticTrack * G4LundStringFragmentation::Sp << 207 Baryon[0][1][1][1]=2214; // Delta+
343 G4Fragmentin << 208 BaryonWeight[0][1][1][1]=pspin_barion;
344 { << 209 //---------------------------------------
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 210
353 //... random choice of string end to us << 211 Baryon[0][1][2][0]=3122; // Lambda
354 G4int SideOfDecay = (G4UniformRand() < << 212 BaryonWeight[0][1][2][0]=(1.-pspin_barion)*0.5;
355 if (SideOfDecay < 0) <<
356 { <<
357 string->SetLeftPartonStable(); <<
358 } else <<
359 { <<
360 string->SetRightPartonStable(); <<
361 } <<
362 213
363 G4ParticleDefinition *newStringEnd; << 214 Baryon[0][1][2][1]=3212; // Sigma0
364 G4ParticleDefinition * HadronDefinition << 215 BaryonWeight[0][1][2][1]=(1.-pspin_barion)*0.5;
365 216
366 G4double StringMass=string->Mass(); << 217 Baryon[0][1][2][2]=3214; // Sigma*0
>> 218 BaryonWeight[0][1][2][2]=pspin_barion;
>> 219 //---------------------------------------
367 220
368 G4double ProbDqADq = GetDiquarkSuppress << 221 Baryon[0][2][0][0]=3112; // Sigma-
369 G4double ProbSaS = 1.0 - 2.0 * GetStr << 222 BaryonWeight[0][2][0][0]=1.-pspin_barion;
370 223
371 #ifdef debug_LUNDfragmentation << 224 Baryon[0][2][0][1]=3114; // Sigma*-
372 G4cout<<"StrMass DiquarkSuppression << 225 BaryonWeight[0][2][0][1]=pspin_barion;
373 #endif << 226 //---------------------------------------
374 227
375 G4int NumberOfpossibleBaryons = 2; << 228 Baryon[0][2][1][0]=3122; // Lambda
>> 229 BaryonWeight[0][2][1][0]=(1.-pspin_barion)*0.5;
376 230
377 if (string->GetLeftParton()->GetParticl << 231 Baryon[0][2][1][1]=3212; // Sigma0
378 if (string->GetRightParton()->GetPartic << 232 BaryonWeight[0][2][1][1]=(1.-pspin_barion)*0.5;
379 233
380 G4double ActualProb = ProbDqADq ; << 234 Baryon[0][2][1][2]=3214; // Sigma*0
381 ActualProb *= (1.0-G4Pow::GetInstance() << 235 BaryonWeight[0][2][1][2]=pspin_barion;
382 if(ActualProb <0.0) ActualProb = 0.; << 236 //---------------------------------------
383 237
384 SetDiquarkSuppression(ActualProb); << 238 Baryon[0][2][2][0]=3312; // Theta-
>> 239 BaryonWeight[0][2][2][0]=1.-pspin_barion;
385 240
386 G4double Mth = 1250.0; << 241 Baryon[0][2][2][1]=3314; // Theta*-
387 if ( NumberOfpossibleBaryons == 3 ){Mth << 242 BaryonWeight[0][2][2][1]=pspin_barion;
388 else if ( NumberOfpossibleBaryons == 4 << 243 //---------------------------------------
389 else {} <<
390 244
391 ActualProb = ProbSaS; << 245 Baryon[1][0][0][0]=2112; // neutron
392 ActualProb *= (1.0 - G4Pow::GetInstance << 246 BaryonWeight[1][0][0][0]=1.-pspin_barion;
393 if ( ActualProb < 0.0 ) ActualProb = 0. <<
394 SetStrangenessSuppression((1.0-ActualPr <<
395 247
396 #ifdef debug_LUNDfragmentation << 248 Baryon[1][0][0][1]=2114; // Delta0
397 G4cout<<"StrMass DiquarkSuppression cor << 249 BaryonWeight[1][0][0][1]=pspin_barion;
398 #endif << 250 //---------------------------------------
399 251
400 if (string->DecayIsQuark()) << 252 Baryon[1][0][1][0]=2212; // proton
401 { << 253 BaryonWeight[1][0][1][0]=1.-pspin_barion;
402 HadronDefinition= QuarkSplitup(strin <<
403 } else { <<
404 HadronDefinition= DiQuarkSplitup(str <<
405 } <<
406 254
407 SetDiquarkSuppression(ProbDqADq); << 255 Baryon[1][0][1][1]=2214; // Delta+
408 SetStrangenessSuppression((1.0-ProbSaS) << 256 BaryonWeight[1][0][1][1]=pspin_barion;
>> 257 //---------------------------------------
409 258
410 if ( HadronDefinition == NULL ) { G4Kin << 259 Baryon[1][0][2][0]=3122; // Lambda
>> 260 BaryonWeight[1][0][2][0]=(1.-pspin_barion)*0.5;
411 261
412 #ifdef debug_LUNDfragmentation << 262 Baryon[1][0][2][1]=3212; // Sigma0
413 G4cout<<"The parton "<<string->GetDecay << 263 BaryonWeight[1][0][2][1]=(1.-pspin_barion)*0.5;
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 264
420 if ( newString ) delete newString; << 265 Baryon[1][0][2][2]=3214; // Sigma*0
>> 266 BaryonWeight[1][0][2][2]=pspin_barion;
>> 267 //---------------------------------------
421 268
422 newString=new G4FragmentingString(*stri << 269 Baryon[1][1][0][0]=2212; // proton
>> 270 BaryonWeight[1][1][0][0]=1.-pspin_barion;
423 271
424 #ifdef debug_LUNDfragmentation << 272 Baryon[1][1][0][1]=2214; // Delta+
425 G4cout<<"An attempt to determine its en << 273 BaryonWeight[1][1][0][1]=pspin_barion;
426 #endif << 274 //---------------------------------------
427 G4LorentzVector* HadronMomentum=SplitEa <<
428 275
429 delete newString; newString=0; << 276 Baryon[1][1][1][0]=2224; // Delta++
430 << 277 BaryonWeight[1][1][1][0]=1.;
431 G4KineticTrack * Hadron =0; << 278 //---------------------------------------
432 if ( HadronMomentum != 0 ) { <<
433 <<
434 #ifdef debug_LUNDfragmentation <<
435 G4cout<<"The attempt was successful <<
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 279
453 #ifdef debug_LUNDfragmentation << 280 Baryon[1][1][2][0]=3222; // Sigma+
454 G4cout<<"End SplitUP (G4VLongitudinalSt << 281 BaryonWeight[1][1][2][0]=1.-pspin_barion;
455 #endif <<
456 282
457 return Hadron; << 283 Baryon[1][1][2][1]=3224; // Sigma*+
458 } << 284 BaryonWeight[1][1][2][1]=pspin_barion;
>> 285 //---------------------------------------
459 286
460 //-------------------------------------------- << 287 Baryon[1][2][0][0]=3122; // Lambda
>> 288 BaryonWeight[1][2][0][0]=(1.-pspin_barion)*0.5;
461 289
462 G4ParticleDefinition * G4LundStringFragmentati << 290 Baryon[1][2][0][1]=3212; // Sigma0
463 << 291 BaryonWeight[1][2][0][1]=(1.-pspin_barion)*0.5;
464 { <<
465 G4double StrSup=GetStrangeSuppress(); <<
466 G4double ProbQQbar = (1.0 - 2.0*StrSup)*1.2 <<
467 292
468 //... can Diquark break or not? << 293 Baryon[1][2][0][2]=3214; // Sigma*0
469 if (G4UniformRand() < DiquarkBreakProb ){ << 294 BaryonWeight[1][2][0][2]=pspin_barion;
>> 295 //---------------------------------------
470 296
471 //... Diquark break << 297 Baryon[1][2][1][0]=3222; // Sigma+
472 G4int stableQuarkEncoding = decay->GetPD << 298 BaryonWeight[1][2][1][0]=1.-pspin_barion;
473 G4int decayQuarkEncoding = (decay->GetPD <<
474 if (G4UniformRand() < 0.5) <<
475 { <<
476 G4int Swap = stableQuarkEncoding; <<
477 stableQuarkEncoding = decayQuarkEncod <<
478 decayQuarkEncoding = Swap; <<
479 } <<
480 299
481 G4int IsParticle=(decayQuarkEncoding>0) << 300 Baryon[1][2][1][1]=3224; // Sigma*+
>> 301 BaryonWeight[1][2][1][1]=pspin_barion;
>> 302 //---------------------------------------
482 303
483 SetStrangenessSuppression((1.0-ProbQQbar << 304 Baryon[1][2][2][0]=3322; // Theta0
484 pDefPair QuarkPair = CreatePartonPair(Is << 305 BaryonWeight[1][2][2][0]=1.-pspin_barion;
485 SetStrangenessSuppression((1.0-StrSup)/2 <<
486 306
487 //... Build new Diquark << 307 Baryon[1][2][2][1]=3324; // Theta*0
488 G4int QuarkEncoding=QuarkPair.second->Ge << 308 BaryonWeight[1][2][2][1]=pspin_barion;
489 G4int i10 = std::max(std::abs(QuarkEnco << 309 //---------------------------------------
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 310
498 return had; << 311 Baryon[2][0][0][0]=3112; // Sigma-
>> 312 BaryonWeight[2][0][0][0]=1.-pspin_barion;
499 313
500 } else { << 314 Baryon[2][0][0][1]=3114; // Sigma*-
501 //... Diquark does not break << 315 BaryonWeight[2][0][0][1]=pspin_barion;
>> 316 //---------------------------------------
502 317
503 G4int IsParticle=(decay->GetPDGEncoding( << 318 Baryon[2][0][1][0]=3122; // Lambda
>> 319 BaryonWeight[2][0][1][0]=(1.-pspin_barion)*0.5;
504 320
505 StrangeSuppress=(1.0 - ProbQQbar)/2.0; << 321 Baryon[2][0][1][1]=3212; // Sigma0
506 pDefPair QuarkPair = CreatePartonPair(Is << 322 BaryonWeight[2][0][1][1]=(1.-pspin_barion)*0.5;
507 323
508 created = QuarkPair.second; << 324 Baryon[2][0][1][2]=3214; // Sigma*0
>> 325 BaryonWeight[2][0][1][2]=pspin_barion;
>> 326 //---------------------------------------
509 327
510 G4ParticleDefinition * had=hadronizer->B << 328 Baryon[2][0][2][0]=3312; // Sigma-
511 StrangeSuppress=StrSup; << 329 BaryonWeight[2][0][2][0]=1.-pspin_barion;
512 330
513 return had; << 331 Baryon[2][0][2][1]=3314; // Sigma*-
514 } << 332 BaryonWeight[2][0][2][1]=pspin_barion;
515 } << 333 //---------------------------------------
516 334
517 //-------------------------------------------- << 335 Baryon[2][1][0][0]=3122; // Lambda
>> 336 BaryonWeight[2][1][0][0]=(1.-pspin_barion)*0.5;
518 337
519 G4LorentzVector * G4LundStringFragmentation::S << 338 Baryon[2][1][0][1]=3212; // Sigma0
520 << 339 BaryonWeight[2][1][0][1]=(1.-pspin_barion)*0.5;
521 <<
522 { <<
523 G4LorentzVector String4Momentum=string->Get4 <<
524 G4double StringMT2=string->MassT2(); <<
525 G4double StringMT =std::sqrt(StringMT2); <<
526 340
527 G4double HadronMass = pHadron->GetPDGMass(); << 341 Baryon[2][1][0][2]=3214; // Sigma*0
528 SetMinimalStringMass(newString); << 342 BaryonWeight[2][1][0][2]=pspin_barion;
>> 343 //---------------------------------------
529 344
530 if ( MinimalStringMass < 0.0 ) return n << 345 Baryon[2][1][1][0]=3222; // Sigma+
>> 346 BaryonWeight[2][1][1][0]=1.-pspin_barion;
531 347
532 #ifdef debug_LUNDfragmentation << 348 Baryon[2][1][1][1]=3224; // Sigma*+
533 G4cout<<G4endl<<"Start LUND SplitEandP << 349 BaryonWeight[2][1][1][1]=pspin_barion;
534 G4cout<<"String 4 mom, String M and Mt << 350 //---------------------------------------
535 <<" "<<std::sqrt(StringMT2)<<G4e <<
536 G4cout<<"Hadron "<<pHadron->GetParticl <<
537 G4cout<<"HadM MinimalStringMassLeft St <<
538 <<String4Momentum.mag()<<" "<<Ha <<
539 #endif <<
540 351
541 if ((HadronMass + MinimalStringMass > string << 352 Baryon[2][1][2][0]=3322; // Theta0
542 { << 353 BaryonWeight[2][1][2][0]=1.-pspin_barion;
543 #ifdef debug_LUNDfragmentation <<
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 354
574 //... sample Pt of the hadron << 355 Baryon[2][1][2][1]=3324; // Theta*0
575 G4int attempt=0; << 356 BaryonWeight[2][1][2][1]=pspin_barion;
576 do << 357 //---------------------------------------
577 { <<
578 attempt++; if (attempt > StringLoopI <<
579 358
580 HadronMt = HadronMass - TmtCur*G4Log << 359 Baryon[2][2][0][0]=3312; // Theta-
581 Pt2 = sqr(HadronMt)-sqr(HadronMass); Pt=st << 360 BaryonWeight[2][2][0][0]=1.-pspin_barion;
582 phi = 2.*pi*G4UniformRand(); <<
583 HadronPt = G4ThreeVector( Pt*std::co <<
584 RemSysPt = StringPt - HadronPt; <<
585 HadronMassT2 = sqr(HadronMass) + Had <<
586 ResidualMassT2=sqr(MinimalStringMass <<
587 361
588 } while (std::sqrt(HadronMassT2) + std << 362 Baryon[2][2][0][1]=3314; // Theta*-
>> 363 BaryonWeight[2][2][0][1]=pspin_barion;
>> 364 //---------------------------------------
589 365
590 //... sample z to define hadron longitudina << 366 Baryon[2][2][1][0]=3322; // Theta0
591 //... but first check the available phase sp << 367 BaryonWeight[2][2][1][0]=1.-pspin_barion;
592 368
593 G4double Pz2 = (sqr(StringMT2 - HadronMassT2 << 369 Baryon[2][2][1][1]=3324; // Theta*0
594 4*HadronMassT2 * ResidualMassT2)/4./Stri << 370 BaryonWeight[2][2][1][1]=pspin_barion;
>> 371 //---------------------------------------
595 372
596 if (Pz2 < 0 ) {return 0;} // have t << 373 Baryon[2][2][2][0]=3334; // Omega
>> 374 BaryonWeight[2][2][2][0]=1.;
>> 375 //---------------------------------------
597 376
598 //... then compute allowed z region z_min < << 377 pspin_barion=pspin_barion_in;
>> 378 /*
>> 379 for(G4int i=0; i<3; i++)
>> 380 { for(G4int j=0; j<3; j++)
>> 381 { for(G4int k=0; k<3; k++)
>> 382 { for(G4int l=0; l<4; l++)
>> 383 { G4cout<<i<<" "<<j<<" "<<k<<" "<<l<<" "<<Baryon[i][j][k][l]<<G4endl;}
>> 384 }
>> 385 }
>> 386 }
>> 387 G4int Uzhi;
>> 388 G4cin>>Uzhi;
>> 389 */
>> 390
>> 391 SetStrangenessSuppression(0.375);
>> 392 Prob_QQbar[0]=StrangeSuppress; // Probability of ddbar production
>> 393 Prob_QQbar[1]=StrangeSuppress; // Probability of uubar production
>> 394 Prob_QQbar[2]=1.0-2.*StrangeSuppress; // Probability of ssbar production
>> 395 SetStrangenessSuppression(0.46); //(0.447);
>> 396
>> 397 //A.R. 25-Jul-2012 : Coverity fix.
>> 398 for ( G4int i=0 ; i<35 ; i++ ) {
>> 399 FS_LeftHadron[i] = 0;
>> 400 FS_RightHadron[i] = 0;
>> 401 FS_Weight[i] = 0.0;
>> 402 }
>> 403 NumberOf_FS = 0;
599 404
600 G4double Pz = std::sqrt(Pz2); << 405 }
601 G4double zMin = (std::sqrt(HadronMassT2+Pz2) <<
602 // G4double zMin = (std::sqrt(HadronMa <<
603 G4double zMax = (std::sqrt(HadronMassT2+Pz2) <<
604 406
605 if (zMin >= zMax) return 0; // have to sta << 407 // --------------------------------------------------------------
>> 408 G4LundStringFragmentation::~G4LundStringFragmentation()
>> 409 {}
606 410
607 G4double z = GetLightConeZ(zMin, zMax, <<
608 string->GetDecayParton()->Get <<
609 HadronPt.x(), HadronPt.y()); <<
610 411
611 //... now compute hadron longitudinal moment << 412 //--------------------------------------------------------------------------------------
612 // longitudinal hadron momentum component Ha << 413 void G4LundStringFragmentation::SetMinimalStringMass(const G4FragmentingString * const string)
>> 414 {
>> 415 G4double EstimatedMass=0.;
>> 416 G4int Number_of_quarks=0;
>> 417 G4int Number_of_squarks=0;
>> 418
>> 419 G4double StringM=string->Get4Momentum().mag();
>> 420
>> 421 G4int Qleft =std::abs(string->GetLeftParton()->GetPDGEncoding());
>> 422
>> 423 #ifdef debug_LUNDfragmentation
>> 424 //G4cout<<"MinStringMass// Input String mass "<<string->Get4Momentum().mag()<<" Qleft "<<Qleft<<G4endl;
>> 425 #endif
>> 426
>> 427 if( Qleft > 1000)
>> 428 {
>> 429 Number_of_quarks+=2;
>> 430 G4int q1=Qleft/1000;
>> 431 if( q1 < 3) {EstimatedMass +=Mass_of_light_quark;}
>> 432 if( q1 > 2) {EstimatedMass +=Mass_of_heavy_quark; Number_of_squarks++;}
>> 433
>> 434 G4int q2=(Qleft/100)%10;
>> 435 if( q2 < 3) {EstimatedMass +=Mass_of_light_quark;}
>> 436 if( q2 > 2) {EstimatedMass +=Mass_of_heavy_quark; Number_of_squarks++;}
>> 437 // EstimatedMass +=Mass_of_string_junction;
>> 438 //if((q1 > 2)||(q2 > 2)) EstimatedMass -= 120*MeV;
>> 439 }
>> 440 else
>> 441 {
>> 442 Number_of_quarks++;
>> 443 if( Qleft < 3) {EstimatedMass +=Mass_of_light_quark;}
>> 444 if( Qleft > 2) {EstimatedMass +=Mass_of_heavy_quark; Number_of_squarks++;}
>> 445 }
>> 446
>> 447 #ifdef debug_LUNDfragmentation
>> 448 //G4cout<<"Min mass with Qleft "<<Qleft<<" "<<EstimatedMass<<G4endl;
>> 449 #endif
>> 450 G4int Qright=std::abs(string->GetRightParton()->GetPDGEncoding());
>> 451 if( Qright > 1000)
>> 452 {
>> 453 Number_of_quarks+=2;
>> 454 G4int q1=Qright/1000;
>> 455 if( q1 < 3) {EstimatedMass +=Mass_of_light_quark;}
>> 456 if( q1 > 2) {EstimatedMass +=Mass_of_heavy_quark; Number_of_squarks++;}
>> 457
>> 458 G4int q2=(Qright/100)%10;
>> 459 if( q2 < 3) {EstimatedMass +=Mass_of_light_quark;}
>> 460 if( q2 > 2) {EstimatedMass +=Mass_of_heavy_quark; Number_of_squarks++;}
>> 461 //EstimatedMass +=Mass_of_string_junction;
>> 462 //if((q1 > 2)||(q2 > 2)) EstimatedMass -= 120*MeV;
>> 463 }
>> 464 else
>> 465 {
>> 466 Number_of_quarks++;
>> 467 if( Qright < 3) {EstimatedMass +=Mass_of_light_quark;}
>> 468 if( Qright > 2) {EstimatedMass +=Mass_of_heavy_quark; Number_of_squarks++;}
>> 469 }
>> 470
>> 471 #ifdef debug_LUNDfragmentation
>> 472 //G4cout<<"Min mass with Qleft and Qright "<<Qright<<" "<<EstimatedMass<<G4endl;
>> 473 //G4cout<<"Number_of_quarks "<<Number_of_quarks<<" Number_of_squarks "<<Number_of_squarks<<G4endl;
>> 474 #endif
>> 475
>> 476 if(Number_of_quarks==2){EstimatedMass += 70.*MeV;} //100.*MeV;}
>> 477 //if(Number_of_quarks==3){EstimatedMass += 20.*MeV;}
>> 478 if(Number_of_quarks==3)
>> 479 {
>> 480 if(Number_of_squarks==0) {EstimatedMass += 740.*MeV;} // 700
>> 481 if(Number_of_squarks==1) {EstimatedMass += 740.*MeV;} // 740
>> 482 if(Number_of_squarks==2) {EstimatedMass += 400.*MeV;}
>> 483 if(Number_of_squarks==3) {EstimatedMass += 382.*MeV;}
>> 484 }
>> 485 if(Number_of_quarks==4)
>> 486 {
>> 487 if((StringM > 1880.) && ( EstimatedMass < 2100)) {EstimatedMass = 2020.;}//1880.;}
>> 488 //if((StringM > 1880.) && ( EstimatedMass < 2100)) {EstimatedMass = 2051.;}
>> 489 else if((StringM > 2232.) && ( EstimatedMass < 2730)){EstimatedMass = 2570.;}
>> 490 else if((StringM > 5130.) && ( EstimatedMass < 3450)){EstimatedMass = 5130.;}
>> 491 else {
>> 492 // EstimatedMass -=2.*Mass_of_string_junction;
>> 493 if(EstimatedMass <= 1600.*MeV){EstimatedMass-=200.*MeV;}
>> 494 else {EstimatedMass+=100.*MeV;}
>> 495 }
>> 496 }
613 497
614 HadronPt.setZ(0.5* string->GetDecayDirection << 498 #ifdef debug_LUNDfragmentation
615 (z * string->LightConeDecay() - << 499 //G4cout<<"EstimatedMass -------------------- "<<EstimatedMass <<G4endl;
616 HadronMassT2/(z * string->LightConeD << 500 #endif
617 G4double HadronE = 0.5* (z * string->LightC << 501 MinimalStringMass=EstimatedMass;
618 HadronMassT2/(z * string->LightConeD << 502 SetMinimalStringMass2(EstimatedMass);
>> 503 }
619 504
620 G4LorentzVector * a4Momentum= new G4LorentzV << 505 //--------------------------------------------------------------------------------------
>> 506 void G4LundStringFragmentation::SetMinimalStringMass2(const G4double aValue)
>> 507 {
>> 508 MinimalStringMass2=aValue * aValue;
>> 509 }
621 510
622 #ifdef debug_LUNDfragmentation << 511 //--------------------------------------------------------------------------------------
623 G4cout<<G4endl<<" string->GetDecayDire << 512 G4KineticTrackVector* G4LundStringFragmentation::FragmentString(const G4ExcitedString& theString)
624 G4cout<<"string->LightConeDecay() "<<s << 513 {
625 G4cout<<"HadronPt,HadronE "<<HadronPt< << 514 // Can no longer modify Parameters for Fragmentation.
626 G4cout<<"String4Momentum "<<String4Mom << 515 PastInitPhase=true;
627 G4cout<<"Out of LUND SplitEandP "<<G4e <<
628 #endif <<
629 516
630 return a4Momentum; << 517 SetMassCut(160.*MeV); // For LightFragmentationTest it is required that no one pi-meson can be produced
631 } <<
632 518
633 //-------------------------------------------- << 519 G4FragmentingString aString(theString);
>> 520 SetMinimalStringMass(&aString);
634 521
635 G4double G4LundStringFragmentation::GetLightCo << 522 #ifdef debug_LUNDfragmentation
636 G4int PD << 523 G4cout<<G4endl<<"LUND StringFragm: String Mass " <<theString.Get4Momentum().mag()<<" Pz "
637 G4Partic << 524 <<theString.Get4Momentum().pz()
638 G4double << 525 <<"------------------------------------"<<G4endl;
639 { << 526 G4cout<<"String ends Direct "<<theString.GetLeftParton()->GetPDGcode()<<" "
640 G4double Mass = pHadron->GetPDGMass(); << 527 <<theString.GetRightParton()->GetPDGcode()<<" "
641 G4int HadronEncoding=std::abs(pHadron- << 528 <<theString.GetDirection()<< G4endl;
>> 529 G4cout<<"Left mom "<<theString.GetLeftParton()->Get4Momentum()<<G4endl;
>> 530 G4cout<<"Right mom "<<theString.GetRightParton()->Get4Momentum()<<G4endl;
>> 531 G4cout<<"Check for Fragmentation "<<G4endl;
>> 532 #endif
>> 533
>> 534 G4KineticTrackVector * LeftVector(0);
>> 535
>> 536 //if(!IsFragmentable(&aString)) // produce 1 hadron True ===============
>> 537 if(!aString.FourQuarkString() && !IsFragmentable(&aString))
>> 538 {
>> 539 #ifdef debug_LUNDfragmentation
>> 540 G4cout<<"Non fragmentable - the string is converted to one hadron "<<G4endl;
>> 541 #endif
>> 542
>> 543 SetMassCut(10000.*MeV);
>> 544 LeftVector=LightFragmentationTest(&theString);
>> 545 SetMassCut(160.*MeV);
>> 546
>> 547 LeftVector->operator[](0)->SetFormationTime(theString.GetTimeOfCreation());
>> 548 LeftVector->operator[](0)->SetPosition(theString.GetPosition());
>> 549
>> 550 if(LeftVector->size() > 1)
>> 551 {
>> 552 // 2 hadrons created from qq-qqbar are stored
>> 553 LeftVector->operator[](1)->SetFormationTime(theString.GetTimeOfCreation());
>> 554 LeftVector->operator[](1)->SetPosition(theString.GetPosition());
>> 555 }
>> 556 return LeftVector;
>> 557 }
642 558
643 G4double Mt2 = Px*Px + Py*Py + Mass*Mass; << 559 #ifdef debug_LUNDfragmentation
>> 560 G4cout<<"The string will be fragmented. "<<G4endl;
>> 561 #endif
>> 562
>> 563 // The string can fragment. At least two particles can be produced.
>> 564 LeftVector =new G4KineticTrackVector;
>> 565 G4KineticTrackVector * RightVector=new G4KineticTrackVector;
>> 566
>> 567 G4ExcitedString *theStringInCMS=CPExcited(theString);
>> 568
>> 569 #ifdef debug_LUNDfragmentation
>> 570 G4cout<<"CMS Left mom "<<theStringInCMS->GetLeftParton()->Get4Momentum()<<G4endl;
>> 571 G4cout<<"CMS Right mom "<<theStringInCMS->GetRightParton()->Get4Momentum()<<G4endl;
>> 572 #endif
>> 573
>> 574 G4LorentzRotation toCms=theStringInCMS->TransformToAlignedCms();
>> 575
>> 576 #ifdef debug_LUNDfragmentation
>> 577 G4cout<<"aligCMS Left mom "<<theStringInCMS->GetLeftParton()->Get4Momentum()<<G4endl;
>> 578 G4cout<<"aligCMS Right mom "<<theStringInCMS->GetRightParton()->Get4Momentum()<<G4endl;
>> 579 G4cout<<G4endl<<"LUND StringFragm: String Mass " <<theStringInCMS->Get4Momentum().mag()<<" Pz Lab "
>> 580 <<theStringInCMS->Get4Momentum().pz()
>> 581 <<"------------------------------------"<<G4endl;
>> 582 G4cout<<"String ends and Direction "<<theStringInCMS->GetLeftParton()->GetPDGcode()<<" "
>> 583 <<theStringInCMS->GetRightParton()->GetPDGcode()<<" "
>> 584 <<theStringInCMS->GetDirection()<< G4endl;
>> 585 #endif
>> 586
>> 587 G4bool success = Loop_toFragmentString(theStringInCMS, LeftVector, RightVector);
>> 588
>> 589 delete theStringInCMS;
>> 590
>> 591 if ( ! success )
>> 592 {
>> 593 std::for_each(LeftVector->begin(), LeftVector->end(), DeleteKineticTrack());
>> 594 LeftVector->clear();
>> 595 std::for_each(RightVector->begin(), RightVector->end(), DeleteKineticTrack());
>> 596 delete RightVector;
>> 597 return LeftVector;
>> 598 }
>> 599
>> 600 // Join Left- and RightVector into LeftVector in correct order.
>> 601 while(!RightVector->empty()) /* Loop checking, 07.08.2015, A.Ribon */
>> 602 {
>> 603 LeftVector->push_back(RightVector->back());
>> 604 RightVector->erase(RightVector->end()-1);
>> 605 }
>> 606 delete RightVector;
>> 607
>> 608 CalculateHadronTimePosition(theString.Get4Momentum().mag(), LeftVector);
>> 609
>> 610 G4LorentzRotation toObserverFrame(toCms.inverse());
>> 611
>> 612 G4double TimeOftheStringCreation=theString.GetTimeOfCreation();
>> 613 G4ThreeVector PositionOftheStringCreation(theString.GetPosition());
>> 614
>> 615 for(size_t C1 = 0; C1 < LeftVector->size(); C1++)
>> 616 {
>> 617 G4KineticTrack* Hadron = LeftVector->operator[](C1);
>> 618 G4LorentzVector Momentum = Hadron->Get4Momentum();
>> 619 //G4cout<<"Hadron "<<Hadron->GetDefinition()->GetParticleName()<<" "<<Momentum<<G4endl;
>> 620 Momentum = toObserverFrame*Momentum;
>> 621 Hadron->Set4Momentum(Momentum);
>> 622
>> 623 G4LorentzVector Coordinate(Hadron->GetPosition(), Hadron->GetFormationTime());
>> 624 Momentum = toObserverFrame*Coordinate;
>> 625 Hadron->SetFormationTime(TimeOftheStringCreation + Momentum.e() - fermi/c_light);
>> 626 G4ThreeVector aPosition(Momentum.vect());
>> 627 Hadron->SetPosition(PositionOftheStringCreation+aPosition);
>> 628 }
644 629
645 G4double Alund, Blund; << 630 return LeftVector;
646 G4double zOfMaxyf(0.), maxYf(1.), z(0.), yf( << 631 }
647 632
648 if (!((std::abs(PDGEncodingOfDecayParton) > << 633 //----------------------------------------------------------------------------------
649 { // ---------------- Quark fragmentation << 634 G4bool G4LundStringFragmentation::IsFragmentable(const G4FragmentingString * const string)
650 Alund=1.; << 635 {
651 Blund=0.7/GeV/GeV; << 636 SetMinimalStringMass(string);
652 << 637 //return sqr(MinimalStringMass + WminLUND) < string->Get4Momentum().mag2();
653 G4double BMt2 = Blund*Mt2; << 638 //G4cout<<"MinM StrM "<<MinimalStringMass<<" "<< string->Get4Momentum().mag()<<G4endl;
654 if (Alund == 1.0) { << 639 return MinimalStringMass < string->Get4Momentum().mag();
655 zOfMaxyf=BMt2/(Blund*Mt2 + 1.);} << 640 }
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 641
679 if (std::abs(PDGEncodingOfDecayParton) > 100 << 642 //----------------------------------------------------------------------------------------
680 { << 643 G4bool G4LundStringFragmentation::StopFragmenting(const G4FragmentingString * const string)
681 G4double an = 2.5; << 644 {
682 an +=(sqr(Px)+sqr(Py))/sqr(GeV << 645 SetMinimalStringMass(string);
683 z=zmin + (zmax-zmin)*G4Pow::Ge <<
684 if( PDGEncodingOfDecayParton > <<
685 } <<
686 646
687 return z; << 647 if (string->FourQuarkString())
>> 648 {
>> 649 return G4UniformRand() < G4Exp(-0.0005*(string->Mass() - MinimalStringMass));
>> 650 } else {
>> 651 // Uzhi 23 Jan. 2015 0.88 -> 0.66
>> 652 return G4UniformRand() < G4Exp(-0.66e-6*(string->Mass()*string->Mass() -
>> 653 MinimalStringMass*MinimalStringMass));
>> 654 }
688 } 655 }
689 656
690 //-------------------------------------------- 657 //----------------------------------------------------------------------------------------------------------
691 <<
692 G4bool G4LundStringFragmentation::SplitLast(G4 658 G4bool G4LundStringFragmentation::SplitLast(G4FragmentingString * string,
693 G4 659 G4KineticTrackVector * LeftVector,
694 G4 660 G4KineticTrackVector * RightVector)
695 { 661 {
696 //... perform last cluster decay << 662 //... perform last cluster decay
697 SetMinimalStringMass( string); <<
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 663
707 G4LorentzVector Str4Mom=string->Get4Mo << 664 #ifdef debug_LUNDfragmentation
708 G4LorentzRotation toCms(-1*Str4Mom.boo << 665 G4cout<<"Split last-----------------------------------------"<<G4endl;
709 G4LorentzVector Pleft = toCms * string << 666 #endif
710 toCms.rotateZ(-1*Pleft.phi()); <<
711 toCms.rotateY(-1*Pleft.theta()); <<
712 <<
713 G4LorentzRotation toObserverFrame= toC <<
714 667
715 G4double StringMass=string->Mass(); << 668 G4LorentzVector Str4Mom=string->Get4Momentum();
>> 669 G4ThreeVector ClusterVel=string->Get4Momentum().boostVector();
>> 670 G4double StringMass=string->Mass();
716 671
717 G4ParticleDefinition * LeftHadron(0), * Righ << 672 G4ParticleDefinition * LeftHadron(0), * RightHadron(0);
718 673
719 NumberOf_FS=0; << 674 NumberOf_FS=0;
720 for (G4int i=0; i<350; i++) {FS_Weight[i]=0. << 675 for(G4int i=0; i<35; i++) {FS_Weight[i]=0.;}
721 676
722 G4int sampledState = 0; << 677 #ifdef debug_LUNDfragmentation
>> 678 G4cout<<"StrMass "<<StringMass<<" q "
>> 679 <<string->GetLeftParton()->GetParticleName()<<" "
>> 680 <<string->GetRightParton()->GetParticleName()<<G4endl;
>> 681 #endif
723 682
724 #ifdef debug_LUNDfragmentation << 683 string->SetLeftPartonStable(); // to query quark contents..
725 G4cout<<"StrMass "<<StringMass<<" q's <<
726 <<string->GetLeftParton()->GetPa <<
727 <<string->GetRightParton()->GetP <<
728 #endif <<
729 684
730 string->SetLeftPartonStable(); // to query q << 685 if (string->FourQuarkString() )
>> 686 {
>> 687 // The string is qq-qqbar type. Diquarks are on the string ends
>> 688 //G4cout<<"The string is qq-qqbar type. Diquarks are on the string ends"<<G4endl;
731 689
732 if (string->IsAFourQuarkString() ) << 690 if(StringMass-MinimalStringMass < 0.)
733 { << 691 {
734 G4int IDleft =std::abs(string->GetLe << 692 if (! Diquark_AntiDiquark_belowThreshold_lastSplitting(string, LeftHadron, RightHadron) )
735 G4int IDright=std::abs(string->GetRi << 693 {
>> 694 return false;
>> 695 }
>> 696 } else {
>> 697 Diquark_AntiDiquark_aboveThreshold_lastSplitting(string, LeftHadron, RightHadron);
736 698
737 if ( (IDleft > 3000) || (IDright > 3 << 699 if(NumberOf_FS == 0) return false;
738 if ( ! Diquark_AntiDiquark_belowTh << 700
739 { << 701 G4int sampledState = SampleState();
740 return false; << 702 if(string->GetLeftParton()->GetPDGEncoding() < 0)
741 } << 703 {
742 } else { << 704 LeftHadron =FS_LeftHadron[sampledState];
743 // The string is qq-qqbar type. Diquarks a << 705 RightHadron=FS_RightHadron[sampledState];
744 if (StringMass-MinimalStringMass < 0 << 706 } else {
745 { << 707 LeftHadron =FS_RightHadron[sampledState];
746 if (! Diquark_AntiDiquark_belowThreshold << 708 RightHadron=FS_LeftHadron[sampledState];
747 { << 709 }
748 return false; << 710 //G4cout<<"Selected "<<SampledState<<" "<<LeftHadron->GetParticleName()<<" "<<RightHadron->GetParticleName()<<G4endl;
749 } << 711 }
750 } else << 712 } else {
751 { << 713 if (string->DecayIsQuark() && string->StableIsQuark() ) { //... there are quarks on cluster ends
752 Diquark_AntiDiquark_aboveThreshold_lastS << 714
753 << 715 #ifdef debug_LUNDfragmentation
754 if (NumberOf_FS == 0) return false; << 716 G4cout<<"Q Q string LastSplit"<<G4endl;
755 << 717 #endif
756 sampledState = SampleS << 718
757 if (string->GetLeftParton()->GetPDGEncod << 719 Quark_AntiQuark_lastSplitting(string, LeftHadron, RightHadron);
758 { << 720 } else { //... there is a Diquark on one of the cluster ends
759 LeftHadron =FS_LeftHadron[sampledState <<
760 RightHadron=FS_RightHadron[sampledStat <<
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 721
775 Quark_AntiQuark_lastSplitting(string, Le << 722 #ifdef debug_LUNDfragmentation
>> 723 G4cout<<"DiQ Q string Last Split"<<G4endl;
>> 724 #endif
>> 725
>> 726 Quark_Diquark_lastSplitting(string, LeftHadron, RightHadron);
>> 727 }
776 728
777 if (NumberOf_FS == 0) return false; << 729 if(NumberOf_FS == 0) return false;
778 sampledState = SampleState() << 730 G4int sampledState = SampleState();
779 if (string->GetLeftParton()->GetPDGEncod << 731 LeftHadron =FS_LeftHadron[sampledState];
780 { << 732 RightHadron=FS_RightHadron[sampledState];
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 <<
794 Quark_Diquark_lastSplitting(string, Left <<
795 <<
796 if (NumberOf_FS == 0) return false; <<
797 sampledState = SampleState() <<
798 <<
799 if (string->GetLeftParton()->GetParticle <<
800 { <<
801 LeftHadron =FS_LeftHadron[sampledState <<
802 RightHadron=FS_RightHadron[sampledStat <<
803 } else <<
804 { <<
805 LeftHadron =FS_RightHadron[sampledStat <<
806 RightHadron=FS_LeftHadron[sampledState <<
807 } <<
808 } <<
809 733
810 } << 734 #ifdef debug_LUNDfragmentation
>> 735 G4cout<<"Selected LeftHad RightHad "<<sampledState<<" "
>> 736 <<LeftHadron->GetParticleName()<<" "<<RightHadron->GetParticleName()<<G4endl;
>> 737 #endif
811 738
812 #ifdef debug_LUNDfragmentation << 739 }
813 G4cout<<"Sampled hadrons: "<<LeftHadro <<
814 #endif <<
815 740
816 G4LorentzVector P_left =string->GetPleft(), << 741 G4LorentzVector LeftMom, RightMom;
>> 742 G4ThreeVector Pos;
817 743
818 G4LorentzVector LeftMom, RightMom; << 744 Sample4Momentum(&LeftMom, LeftHadron->GetPDGMass(), &RightMom, RightHadron->GetPDGMass(), StringMass);
819 G4ThreeVector Pos; <<
820 745
821 Sample4Momentum(&LeftMom, LeftHadron->GetPD << 746 LeftMom.boost(ClusterVel);
822 &RightMom, RightHadron->GetPDGMass(), << 747 RightMom.boost(ClusterVel);
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 } <<
840 748
841 LeftMom *=toObserverFrame; << 749 LeftVector->push_back(new G4KineticTrack(LeftHadron, 0, Pos, LeftMom));
842 RightMom*=toObserverFrame; << 750 RightVector->push_back(new G4KineticTrack(RightHadron, 0, Pos, RightMom));
843 751
844 LeftVector->push_back(new G4KineticTrack(Lef << 752 return true;
845 RightVector->push_back(new G4KineticTrack(Ri << 753 }
846 754
847 string->LorentzRotate(toObserverFrame); << 755 //----------------------------------------------------------------------------------------------------------
848 return true; << 756 void G4LundStringFragmentation::Sample4Momentum(G4LorentzVector* Mom, G4double Mass,
>> 757 G4LorentzVector* AntiMom, G4double AntiMass,
>> 758 G4double InitialMass)
>> 759 {
>> 760 // ------ Sampling of momenta of 2 last produced hadrons --------------------
>> 761 G4ThreeVector Pt;
>> 762 G4double MassMt2, AntiMassMt2;
>> 763 G4double AvailablePz, AvailablePz2;
>> 764 //AR-Oct2017 G4double ProbIsotropy = sqr(sqr(938.0/InitialMass));
>> 765
>> 766 #ifdef debug_LUNDfragmentation
>> 767 G4cout<<"Sampling of momenta of 2 last produced hadrons ----------------"<<G4endl;
>> 768 G4cout<<"Masses "<<InitialMass<<" "<<Mass<<" "<<AntiMass<<G4endl;
>> 769 #endif
>> 770
>> 771 //AR-Oct2017 if((Mass > 930. || AntiMass > 930.) && (G4UniformRand() < ProbIsotropy)) { //If there is a baryon
>> 772 //AR-Oct2017 // ----------------- Isotropic decay ------------------------------------
>> 773 //AR-Oct2017 G4double r_val = sqr(InitialMass*InitialMass - Mass*Mass - AntiMass*AntiMass) - sqr(2.*Mass*AntiMass);
>> 774 //AR-Oct2017 G4double Pabs = (r_val > 0.)? std::sqrt(r_val)/(2.*InitialMass) : 0;
>> 775 //AR-Oct2017 G4cout<<"P for isotr decay "<<Pabs<<G4endl;
>> 776 //AR-Oct2017 //... sample unit vector
>> 777 //AR-Oct2017 G4double pz =1. - 2.*G4UniformRand();
>> 778 //AR-Oct2017 G4double st = std::sqrt(1. - pz * pz)*Pabs;
>> 779 //AR-Oct2017 G4double phi = 2.*pi*G4UniformRand();
>> 780 //AR-Oct2017 G4double px = st*std::cos(phi);
>> 781 //AR-Oct2017 G4double py = st*std::sin(phi);
>> 782 //AR-Oct2017 pz *= Pabs;
>> 783 //AR-Oct2017 Mom->setPx(px); Mom->setPy(py); Mom->setPz(pz);
>> 784 //AR-Oct2017 Mom->setE(std::sqrt(Pabs*Pabs + Mass*Mass));
>> 785 //AR-Oct2017 AntiMom->setPx(-px); AntiMom->setPy(-py); AntiMom->setPz(-pz);
>> 786 //AR-Oct2017 AntiMom->setE (std::sqrt(Pabs*Pabs + AntiMass*AntiMass));
>> 787 //AR-Oct2017 //G4int Uzhi; G4cin>>Uzhi;
>> 788 //AR-Oct2017 } else {
>> 789
>> 790 const G4int maxNumberOfLoops = 1000;
>> 791
>> 792 G4double SigmaQTw = SigmaQT; //AR-Oct2017
>> 793 if ( Mass > 930. && AntiMass > 930.) SigmaQT *=(1.0-0.45*sqr((Mass+AntiMass)/InitialMass)); //AR-Oct2017
>> 794 // In the above condition, "&&" is used instead of "||" (as in development) to be conservative
>> 795 // i.e. to be applied (for the time being) only in annihilations.
>> 796
>> 797 G4int loopCounter = 0;
>> 798 do
>> 799 {
>> 800 // GF 22-May-09, limit sampled pt to allowed range
>> 801
>> 802 G4double termD = InitialMass*InitialMass -Mass*Mass - AntiMass*AntiMass;
>> 803 G4double termab = 4*sqr(Mass*AntiMass);
>> 804 G4double termN = 2*termD + 4*Mass*Mass + 4*AntiMass*AntiMass;
>> 805 G4double pt2max=(termD*termD - termab )/ termN ;
>> 806 //G4cout<<"Anis "<<pt2max<<" "<<(termD*termD-termab)/(4.*InitialMass*InitialMass)<<G4endl;
>> 807
>> 808 Pt=SampleQuarkPt(std::sqrt(pt2max)); Pt.setZ(0); G4double Pt2=Pt.mag2();
>> 809 //G4cout<<"Sampl pt2 "<<Pt2<<G4endl;
>> 810 MassMt2 = Mass * Mass + Pt2;
>> 811 AntiMassMt2= AntiMass * AntiMass + Pt2;
>> 812
>> 813 AvailablePz2= sqr(InitialMass*InitialMass - MassMt2 - AntiMassMt2) - 4.*MassMt2*AntiMassMt2;
>> 814
>> 815 } while( (AvailablePz2 < 0.) && // GF will occur only for numerical precision problem with limit in sampled pt
>> 816 ++loopCounter < maxNumberOfLoops ); /* Loop checking, 07.08.2015, A.Ribon */
>> 817
>> 818 if ( loopCounter >= maxNumberOfLoops ) {
>> 819 AvailablePz2 = 0.0;
>> 820 }
>> 821
>> 822 AvailablePz2 /=(4.*InitialMass*InitialMass);
>> 823 AvailablePz = std::sqrt(AvailablePz2);
>> 824
>> 825 G4double Px=Pt.getX();
>> 826 G4double Py=Pt.getY();
>> 827
>> 828 Mom->setPx(Px); Mom->setPy(Py); Mom->setPz(AvailablePz);
>> 829 Mom->setE(std::sqrt(MassMt2+AvailablePz2));
>> 830
>> 831 AntiMom->setPx(-Px); AntiMom->setPy(-Py); AntiMom->setPz(-AvailablePz);
>> 832 AntiMom->setE (std::sqrt(AntiMassMt2+AvailablePz2));
>> 833
>> 834 if ( Mass > 930. && AntiMass > 930.) SigmaQT = SigmaQTw; //AR-Oct2017
>> 835 // In the above condition, "&&" is used instead of "||" (as in development) to be conservative
>> 836 // i.e. to be applied (for the time being) only in annihilations.
>> 837
>> 838 //AR-Oct2017 }
849 } 839 }
850 840
>> 841 //-----------------------------------------------------------------------------
>> 842 G4LorentzVector * G4LundStringFragmentation::SplitEandP(G4ParticleDefinition * pHadron,
>> 843 G4FragmentingString * string,
>> 844 G4FragmentingString * newString)
>> 845 {
>> 846 G4LorentzVector String4Momentum=string->Get4Momentum();
>> 847 G4double StringMT2=string->MassT2();
>> 848 G4double StringMT =std::sqrt(StringMT2);
>> 849
>> 850 G4double HadronMass = pHadron->GetPDGMass();
>> 851
>> 852 SetMinimalStringMass(newString);
>> 853
>> 854 #ifdef debug_LUNDfragmentation
>> 855 G4cout<<G4endl<<"Start LUND SplitEandP "<<G4endl;
>> 856 G4cout<<"String 4 mom, String M and Mt "<<String4Momentum<<" "<<String4Momentum.mag()<<" "<<std::sqrt(StringMT2)<<G4endl;
>> 857 G4cout<<"Hadron "<<pHadron->GetParticleName()<<G4endl;
>> 858 G4cout<<"HadM MinimalStringMassLeft StringM hM+sM "<<HadronMass<<" "<<MinimalStringMass<<" "
>> 859 <<String4Momentum.mag()<<" "<<HadronMass+MinimalStringMass<<G4endl;
>> 860 #endif
>> 861
>> 862 if(HadronMass + MinimalStringMass > string->Mass())
>> 863 {
>> 864 #ifdef debug_LUNDfragmentation
>> 865 G4cout<<"Mass of the string is not sufficient to produce the hadron!"<<G4endl;
>> 866 #endif
>> 867
>> 868 return 0;
>> 869 } // have to start all over!
>> 870
>> 871 String4Momentum.setPz(0.);
>> 872 G4ThreeVector StringPt=String4Momentum.vect();
>> 873
>> 874 // calculate and assign hadron transverse momentum component HadronPx and HadronPy
>> 875 G4ThreeVector HadronPt , RemSysPt;
>> 876 G4double HadronMassT2, ResidualMassT2;
>> 877
>> 878 //... sample Pt of the hadron
>> 879 G4int attempt=0;
>> 880 do
>> 881 {
>> 882 attempt++; if(attempt > StringLoopInterrupt) return 0;
>> 883
>> 884 HadronPt =SampleQuarkPt() + string->DecayPt();
>> 885 HadronPt.setZ(0);
>> 886 RemSysPt = StringPt - HadronPt;
>> 887
>> 888 HadronMassT2 = sqr(HadronMass) + HadronPt.mag2();
>> 889 ResidualMassT2=sqr(MinimalStringMass) + RemSysPt.mag2();
>> 890
>> 891 } while(std::sqrt(HadronMassT2) + std::sqrt(ResidualMassT2) > StringMT); /* Loop checking, 07.08.2015, A.Ribon */
>> 892
>> 893 //... sample z to define hadron longitudinal momentum and energy
>> 894 //... but first check the available phase space
>> 895
>> 896 G4double Pz2 = (sqr(StringMT2 - HadronMassT2 - ResidualMassT2) -
>> 897 4*HadronMassT2 * ResidualMassT2)/4./StringMT2;
>> 898
>> 899 if(Pz2 < 0 ) {return 0;} // have to start all over!
>> 900
>> 901 //... then compute allowed z region z_min <= z <= z_max
>> 902
>> 903 G4double Pz = std::sqrt(Pz2);
>> 904 G4double zMin = (std::sqrt(HadronMassT2+Pz2) - Pz)/std::sqrt(StringMT2);
>> 905 //G4double zMin = (std::sqrt(HadronMassT2+Pz2) - 0.)/std::sqrt(StringMT2);
>> 906 G4double zMax = (std::sqrt(HadronMassT2+Pz2) + Pz)/std::sqrt(StringMT2);
>> 907
>> 908 if (zMin >= zMax) return 0; // have to start all over!
>> 909
>> 910 G4double z = GetLightConeZ(zMin, zMax, string->GetDecayParton()->GetPDGEncoding(),
>> 911 pHadron, HadronPt.x(), HadronPt.y());
>> 912
>> 913 //... now compute hadron longitudinal momentum and energy
>> 914 // longitudinal hadron momentum component HadronPz
>> 915
>> 916 HadronPt.setZ(0.5* string->GetDecayDirection() *
>> 917 (z * string->LightConeDecay() - HadronMassT2/(z * string->LightConeDecay())));
>> 918 G4double HadronE = 0.5* (z * string->LightConeDecay() +
>> 919 HadronMassT2/(z * string->LightConeDecay()));
>> 920
>> 921 G4LorentzVector * a4Momentum= new G4LorentzVector(HadronPt,HadronE);
>> 922
>> 923 #ifdef debug_LUNDfragmentation
>> 924 G4cout<<"string->LightConeDecay() "<<string->LightConeDecay()<<G4endl;
>> 925 G4cout<<"HadronPt,HadronE "<<HadronPt<<" "<<HadronE<<G4endl;
>> 926 G4cout<<"String4Momentum "<<String4Momentum<<G4endl;
>> 927 //G4int Uzhi; G4cin>>Uzhi;
>> 928 G4cout<<"Out of LUND SplitEandP "<<G4endl;
>> 929 #endif
>> 930
>> 931 return a4Momentum;
>> 932 }
>> 933
>> 934 //-----------------------------------------------------------------------------------------
>> 935 G4double G4LundStringFragmentation::
>> 936 GetLightConeZ(G4double zmin, G4double zmax, G4int PDGEncodingOfDecayParton,
>> 937 G4ParticleDefinition* pHadron, G4double Px, G4double Py)
>> 938 {
>> 939 G4double Mass = pHadron->GetPDGMass();
>> 940 //G4int HadronEncoding=std::abs(pHadron->GetPDGEncoding());
>> 941
>> 942 G4double Mt2 = Px*Px + Py*Py + Mass*Mass;
>> 943
>> 944 G4double alund;
>> 945 G4double zOfMaxyf(0.), maxYf(1.), z(0.), yf(1.);
>> 946 if(std::abs(PDGEncodingOfDecayParton) < 1000)
>> 947 { // ---------------- Quark fragmentation ----------------------
>> 948 alund=0.7/GeV/GeV;
>> 949 // If blund get restored, you MUST adapt the calculation of zOfMaxyf.
>> 950 // const G4double blund = 1;
>> 951
>> 952 zOfMaxyf=alund*Mt2/(alund*Mt2 + 1.);
>> 953 maxYf=(1-zOfMaxyf)/zOfMaxyf * G4Exp(-alund*Mt2/zOfMaxyf);
>> 954
>> 955 const G4int maxNumberOfLoops = 1000;
>> 956 G4int loopCounter = 0;
>> 957 do
>> 958 {
>> 959 z = zmin + G4UniformRand()*(zmax-zmin);
>> 960 yf = (1-z)/z * G4Exp(-alund*Mt2/z);
>> 961 //yf = G4Pow::GetInstance()->powA(1.0-z,blund)/z*G4Exp(-alund*Mt2/z
>> 962 } while ( (G4UniformRand()*maxYf > yf) && ++loopCounter < maxNumberOfLoops ); /* Loop checking, 07.08.2015, A.Ribon */
>> 963 if ( loopCounter >= maxNumberOfLoops ) {
>> 964 z = 0.5*(zmin + zmax); // Just a value between zmin and zmax, no physics considerations at all!
>> 965 }
>> 966 return z;
>> 967 }
>> 968
>> 969 if(std::abs(PDGEncodingOfDecayParton) > 1000)
>> 970 {
>> 971 /*
>> 972 if(HadronEncoding < 3000)
>> 973 {
>> 974 maxYf=(zmax-zmin);
>> 975 do
>> 976 {
>> 977 z = zmin + G4UniformRand()*(zmax-zmin);
>> 978 //yf=G4Exp(-sqr(z-Zc)/2/sqr(0.28)); // 0.42 0.632 0.28 a'la UrQMD
>> 979 yf =(z-zmin);
>> 980 } while (G4UniformRand()*maxYf > yf);
>> 981 } else { // Strange baryons
>> 982 z = zmin + G4UniformRand()*(zmax-zmin);
>> 983 }
>> 984 */
>> 985 z = zmin + G4UniformRand()*(zmax-zmin);
>> 986 }
>> 987 return z;
>> 988 }
>> 989
>> 990 //------------------------------------------------------------------------
>> 991 G4double G4LundStringFragmentation::lambda(G4double S, G4double m1_Sqr, G4double m2_Sqr)
>> 992 {
>> 993 G4double lam = sqr(S - m1_Sqr - m2_Sqr) - 4.*m1_Sqr*m2_Sqr;
>> 994 return lam;
>> 995 }
>> 996
>> 997
>> 998 //------------------------------------------------------------------------
>> 999 //------------------------------------------------------------------------
>> 1000 // Internal methods introduced to improve the code structure (AR Nov 2011)
>> 1001 //------------------------------------------------------------------------
>> 1002 //------------------------------------------------------------------------
>> 1003
851 //-------------------------------------------- 1004 //----------------------------------------------------------------------------------------
>> 1005 G4bool G4LundStringFragmentation::Loop_toFragmentString(G4ExcitedString * & theStringInCMS,
>> 1006 G4KineticTrackVector * & LeftVector,
>> 1007 G4KineticTrackVector * & RightVector)
>> 1008 {
>> 1009 #ifdef debug_LUNDfragmentation
>> 1010 G4cout<<"Loop_toFrag "<<theStringInCMS->GetLeftParton()->GetPDGcode()<<" "
>> 1011 <<theStringInCMS->GetRightParton()->GetPDGcode()<<" "
>> 1012 <<theStringInCMS->GetDirection()<< G4endl;
>> 1013 #endif
>> 1014
>> 1015 G4bool final_success=false;
>> 1016 G4bool inner_success=true;
>> 1017 G4int attempt=0;
>> 1018 while ( ! final_success && attempt++ < StringLoopInterrupt ) /* Loop checking, 07.08.2015, A.Ribon */
>> 1019 { // If the string fragmentation do not be happend, repeat the fragmentation.
>> 1020
>> 1021 G4FragmentingString *currentString=new G4FragmentingString(*theStringInCMS);
>> 1022 //G4cout<<"Main loop start whilecounter "<<attempt<<G4endl;
>> 1023 // Cleaning up the previously produced hadrons
>> 1024 std::for_each(LeftVector->begin(), LeftVector->end(), DeleteKineticTrack());
>> 1025 LeftVector->clear();
>> 1026 std::for_each(RightVector->begin(), RightVector->end(), DeleteKineticTrack());
>> 1027 RightVector->clear();
>> 1028
>> 1029 // Main fragmentation loop until the string will not be able to fragment
>> 1030 inner_success=true; // set false on failure.
>> 1031 const G4int maxNumberOfLoops = 1000;
>> 1032 G4int loopCounter = -1;
>> 1033 while ( (! StopFragmenting(currentString)) && ++loopCounter < maxNumberOfLoops ) /* Loop checking, 07.08.2015, A.Ribon */
>> 1034 { // Split current string into hadron + new string
>> 1035 #ifdef debug_LUNDfragmentation
>> 1036 G4cout<<"The string can fragment. "<<G4endl;;
>> 1037 //G4cout<<"1 "<<currentString->GetDecayDirection()<<G4endl;
>> 1038 #endif
>> 1039 G4FragmentingString *newString=0; // used as output from SplitUp.
>> 1040 G4KineticTrack * Hadron=Splitup(currentString,newString);
>> 1041 if ( Hadron != 0 ) // Store the hadron
>> 1042 {
>> 1043 #ifdef debug_LUNDfragmentation
>> 1044 G4cout<<"Hadron prod at fragm. "<<Hadron->GetDefinition()->GetParticleName()<<G4endl;
>> 1045 //G4cout<<"2 "<<currentString->GetDecayDirection()<<G4endl;
>> 1046 #endif
852 1047
>> 1048 if ( currentString->GetDecayDirection() > 0 ) {
>> 1049 LeftVector->push_back(Hadron);
>> 1050 } else {
>> 1051 RightVector->push_back(Hadron);
>> 1052 }
>> 1053 delete currentString;
>> 1054 currentString=newString;
>> 1055 }
>> 1056 };
>> 1057 if ( loopCounter >= maxNumberOfLoops ) {
>> 1058 inner_success=false;
>> 1059 }
>> 1060
>> 1061 // Split remaining string into 2 final hadrons.
>> 1062 #ifdef debug_LUNDfragmentation
>> 1063 G4cout<<"Split remaining string into 2 final hadrons."<<G4endl;
>> 1064 #endif
>> 1065
>> 1066 if ( inner_success && SplitLast(currentString, LeftVector, RightVector) )
>> 1067 {
>> 1068 final_success=true;
>> 1069 }
>> 1070
>> 1071 //final_success=true;
>> 1072 delete currentString;
>> 1073 } // End of the loop where we try to fragment the string.
>> 1074 return final_success;
>> 1075 }
>> 1076
>> 1077 //----------------------------------------------------------------------------------------
853 G4bool G4LundStringFragmentation:: 1078 G4bool G4LundStringFragmentation::
854 Diquark_AntiDiquark_belowThreshold_lastSplitti 1079 Diquark_AntiDiquark_belowThreshold_lastSplitting(G4FragmentingString * & string,
855 1080 G4ParticleDefinition * & LeftHadron,
856 1081 G4ParticleDefinition * & RightHadron)
857 { 1082 {
858 G4double StringMass = string->Mass(); << 1083 G4double StringMass = string->Mass();
859 << 1084 G4int cClusterInterrupt = 0;
860 G4int cClusterInterrupt = 0; << 1085 do
861 G4bool isOK = false; << 1086 {
862 do << 1087 //G4cout<<"cClusterInterrupt "<<cClusterInterrupt<<G4endl;
863 { << 1088 if (cClusterInterrupt++ >= ClusterLoopInterrupt) return false;
864 G4int LeftQuark1= string->GetLeftParton()- << 1089
865 G4int LeftQuark2=(string->GetLeftParton()- << 1090 G4int LeftQuark1= string->GetLeftParton()->GetPDGEncoding()/1000;
>> 1091 G4int LeftQuark2=(string->GetLeftParton()->GetPDGEncoding()/100)%10;
>> 1092
>> 1093 G4int RightQuark1= string->GetRightParton()->GetPDGEncoding()/1000;
>> 1094 G4int RightQuark2=(string->GetRightParton()->GetPDGEncoding()/100)%10;
>> 1095
>> 1096 if(G4UniformRand()<0.5) {
>> 1097 LeftHadron =hadronizer->Build(FindParticle( LeftQuark1), FindParticle(RightQuark1));
>> 1098 RightHadron=hadronizer->Build(FindParticle( LeftQuark2), FindParticle(RightQuark2));
>> 1099 } else {
>> 1100 LeftHadron =hadronizer->Build(FindParticle( LeftQuark1), FindParticle(RightQuark2));
>> 1101 RightHadron=hadronizer->Build(FindParticle( LeftQuark2), FindParticle(RightQuark1));
>> 1102 }
866 1103
867 G4int RightQuark1= string->GetRightParton( << 1104 //... repeat procedure, if mass of cluster is too low to produce hadrons
868 G4int RightQuark2=(string->GetRightParton( << 1105 //... ClusterMassCut = 0.15*GeV model parameter
>> 1106 }
>> 1107 while ((StringMass <= LeftHadron->GetPDGMass() + RightHadron->GetPDGMass())); /* Loop checking, 07.08.2015, A.Ribon */
869 1108
870 if (G4UniformRand()<0.5) << 1109 return true;
871 { <<
872 LeftHadron =hadronizer->Build(FindPartic <<
873 FindParticle(RightQuark1)); <<
874 RightHadron= (LeftHadron == nullptr) ? n <<
875 <<
876 FindParticle(RightQuark2)); <<
877 } else <<
878 { <<
879 LeftHadron =hadronizer->Build(FindPartic <<
880 FindParticle(RightQuark2)); <<
881 RightHadron=(LeftHadron == nullptr) ? nu <<
882 hadronizer <<
883 FindParticle(RightQuark1)); <<
884 } <<
885 <<
886 isOK = (LeftHadron != nullptr) && (RightHa <<
887 <<
888 if(isOK) { isOK = (StringMass > LeftHadron <<
889 ++cClusterInterrupt; <<
890 //... repeat procedure, if mass of cluster <<
891 //... ClusterMassCut = 0.15*GeV model para <<
892 } <<
893 while (isOK == false && cClusterInterrupt < <<
894 /* Loop checking, 07.08.2015, A.Ribon */ <<
895 return isOK; <<
896 } 1110 }
897 1111
898 //-------------------------------------------- 1112 //----------------------------------------------------------------------------------------
899 <<
900 G4bool G4LundStringFragmentation:: 1113 G4bool G4LundStringFragmentation::
901 Diquark_AntiDiquark_aboveThreshold_lastSplitti 1114 Diquark_AntiDiquark_aboveThreshold_lastSplitting(G4FragmentingString * & string,
902 1115 G4ParticleDefinition * & LeftHadron,
903 1116 G4ParticleDefinition * & RightHadron)
904 { 1117 {
905 // StringMass-MinimalStringMass > 0. Creatio << 1118 // StringMass-MinimalStringMass > 0. Creation of 2 baryons is possible ----
906 <<
907 G4double StringMass = string->Mass(); <<
908 G4double StringMassSqr= sqr(StringMass); <<
909 G4ParticleDefinition * Di_Quark; <<
910 G4ParticleDefinition * Anti_Di_Quark; <<
911 1119
912 if (string->GetLeftParton()->GetPDGEncoding( << 1120 G4double StringMass = string->Mass();
913 { << 1121 G4double StringMassSqr= sqr(StringMass);
914 Anti_Di_Quark =string->GetLeftParton(); << 1122 G4ParticleDefinition * Di_Quark;
915 Di_Quark=string->GetRightParton(); << 1123 G4ParticleDefinition * Anti_Di_Quark;
916 } else << 1124
917 { << 1125 if(string->GetLeftParton()->GetPDGEncoding() < 0) {
918 Anti_Di_Quark =string->GetRightParton(); << 1126 Anti_Di_Quark =string->GetLeftParton();
919 Di_Quark=string->GetLeftParton(); << 1127 Di_Quark=string->GetRightParton();
920 } << 1128 } else {
921 << 1129 Anti_Di_Quark =string->GetRightParton();
922 G4int IDAnti_di_quark =Anti_Di_Quark->Get << 1130 Di_Quark=string->GetLeftParton();
923 G4int AbsIDAnti_di_quark =std::abs(IDAnti_di << 1131 }
924 G4int IDdi_quark =Di_Quark->GetPDGEn << 1132
925 G4int AbsIDdi_quark =std::abs(IDdi_quar << 1133 G4int IDAnti_di_quark =Anti_Di_Quark->GetPDGEncoding();
926 << 1134 G4int AbsIDAnti_di_quark =std::abs(IDAnti_di_quark);
927 G4int ADi_q1=AbsIDAnti_di_quark/1000; << 1135 G4int IDdi_quark =Di_Quark->GetPDGEncoding();
928 G4int ADi_q2=(AbsIDAnti_di_quark-ADi_q1*1000 << 1136 G4int AbsIDdi_quark =std::abs(IDdi_quark);
929 << 1137
930 G4int Di_q1=AbsIDdi_quark/1000; << 1138 G4int ADi_q1=AbsIDAnti_di_quark/1000;
931 G4int Di_q2=(AbsIDdi_quark-Di_q1*1000)/100; << 1139 G4int ADi_q2=(AbsIDAnti_di_quark-ADi_q1*1000)/100;
>> 1140
>> 1141 G4int Di_q1=AbsIDdi_quark/1000;
>> 1142 G4int Di_q2=(AbsIDdi_quark-Di_q1*1000)/100;
>> 1143
>> 1144 NumberOf_FS=0;
>> 1145 for(G4int ProdQ=1; ProdQ < 4; ProdQ++)
>> 1146 {
>> 1147 G4int StateADiQ=0;
>> 1148 const G4int maxNumberOfLoops = 1000;
>> 1149 G4int loopCounter = 0;
>> 1150 do
>> 1151 {
>> 1152 LeftHadron=G4ParticleTable::GetParticleTable()->FindParticle(
>> 1153 -Baryon[ADi_q1-1][ADi_q2-1][ProdQ-1][StateADiQ]);
>> 1154 G4double LeftHadronMass=LeftHadron->GetPDGMass();
>> 1155
>> 1156 //G4cout<<"Anti Bar "<<LeftHadron->GetParticleName()<<G4endl;
>> 1157
>> 1158 G4int StateDiQ=0;
>> 1159 const G4int maxNumberOfInternalLoops = 1000;
>> 1160 G4int internalLoopCounter = 0;
>> 1161 do
>> 1162 {
>> 1163 RightHadron=G4ParticleTable::GetParticleTable()->FindParticle(
>> 1164 +Baryon[Di_q1-1][Di_q2-1][ProdQ-1][StateDiQ]);
>> 1165 G4double RightHadronMass=RightHadron->GetPDGMass();
932 1166
933 NumberOf_FS=0; << 1167 if(StringMass > LeftHadronMass + RightHadronMass)
934 for (G4int ProdQ=1; ProdQ < 6; ProdQ++) <<
935 { 1168 {
936 G4int StateADiQ=0; << 1169 if ( NumberOf_FS > 34 ) {
937 const G4int maxNumberOfLoops = << 1170 G4ExceptionDescription ed;
938 G4int loopCounter = 0; << 1171 ed << " NumberOf_FS exceeds its limit: NumberOf_FS=" << NumberOf_FS << G4endl;
939 do // while(Meson[AbsIDquark-1][ProdQ-1][ << 1172 G4Exception( "G4LundStringFragmentation::Diquark_AntiDiquark_aboveThreshold_lastSplitting ",
940 { << 1173 "HAD_LUND_001", JustWarning, ed );
941 LeftHadron=G4ParticleTable::GetParticleT << 1174 NumberOf_FS = 34;
942 -Baryon[ADi_q1-1][ADi_q2-1][Prod << 1175 }
943 << 1176
944 if (LeftHadron == NULL) continue; << 1177 G4double FS_Psqr=lambda(StringMassSqr,sqr(LeftHadronMass), sqr(RightHadronMass));
945 G4double LeftHadronMass=LeftHadron->GetP << 1178 //FS_Psqr=1.;
946 << 1179 FS_Weight[NumberOf_FS]=std::sqrt(FS_Psqr)*FS_Psqr*
947 G4int StateDiQ=0; << 1180 BaryonWeight[ADi_q1-1][ADi_q2-1][ProdQ-1][StateADiQ]*
948 const G4int maxNumberO << 1181 BaryonWeight[Di_q1-1][Di_q2-1][ProdQ-1][StateDiQ]*
949 G4int internalLoopCoun << 1182 Prob_QQbar[ProdQ-1];
950 do // while(Baryon[Di_q1-1][Di_q2-1][Pro << 1183
951 { << 1184 FS_LeftHadron[NumberOf_FS] = LeftHadron;
952 RightHadron=G4ParticleTable::GetPartic << 1185 FS_RightHadron[NumberOf_FS]= RightHadron;
953 +Baryon[Di_q1-1][Di_q2-1][Prod << 1186
954 << 1187 NumberOf_FS++;
955 if (RightHadron == NULL) continue; << 1188 }
956 G4double RightHadronMass=RightHadron-> << 1189
957 << 1190 StateDiQ++;
958 if (StringMass > LeftHadronMass + Righ << 1191
959 { << 1192 } while( (Baryon[Di_q1-1][Di_q2-1][ProdQ-1][StateDiQ]!=0) &&
960 if ( N << 1193 ++internalLoopCounter < maxNumberOfInternalLoops ); /* Loop checking, 07.08.2015, A.Ribon */
961 G4Ex << 1194 if ( internalLoopCounter >= maxNumberOfInternalLoops ) {
962 ed < << 1195 return false;
963 G4Ex << 1196 }
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 1197
989 StateADiQ++; << 1198 StateADiQ++;
990 } while( (Baryon[ADi_q1-1][ADi_q2-1][ProdQ << 1199 } while( (Baryon[ADi_q1-1][ADi_q2-1][ProdQ-1][StateADiQ]!=0) &&
991 ++loopCounter < maxNu << 1200 ++loopCounter < maxNumberOfLoops ); /* Loop checking, 07.08.2015, A.Ribon */
992 if ( loopCounter >= maxNumberO << 1201 if ( loopCounter >= maxNumberOfLoops ) {
993 return false; << 1202 return false;
994 } << 1203 }
995 } // End of for (G4int ProdQ=1; ProdQ < 4; P << 1204 }
996 1205
997 return true; << 1206 return true;
998 } 1207 }
999 1208
1000 //------------------------------------------- 1209 //----------------------------------------------------------------------------------------
1001 << 1210 G4bool G4LundStringFragmentation::
1002 G4bool G4LundStringFragmentation::Quark_Diqua << 1211 Quark_Diquark_lastSplitting(G4FragmentingString * & string,
1003 << 1212 G4ParticleDefinition * & LeftHadron,
1004 << 1213 G4ParticleDefinition * & RightHadron)
1005 { 1214 {
1006 G4double StringMass = string->Mass(); << 1215 G4double StringMass = string->Mass();
1007 G4double StringMassSqr= sqr(StringMass); << 1216 G4double StringMassSqr= sqr(StringMass);
1008 1217
1009 G4ParticleDefinition * Di_Quark; << 1218 G4ParticleDefinition * Di_Quark;
1010 G4ParticleDefinition * Quark; << 1219 G4ParticleDefinition * Quark;
1011 1220
1012 if (string->GetLeftParton()->GetParticleSub << 1221 if(string->GetLeftParton()->GetParticleSubType()== "quark") {
1013 { << 1222 Quark =string->GetLeftParton();
1014 Quark =string->GetLeftParton(); << 1223 Di_Quark=string->GetRightParton();
1015 Di_Quark=string->GetRightParton(); << 1224 } else {
1016 } else << 1225 Quark =string->GetRightParton();
1017 { << 1226 Di_Quark=string->GetLeftParton();
1018 Quark =string->GetRightParton(); << 1227 }
1019 Di_Quark=string->GetLeftParton(); << 1228
1020 } << 1229 G4int IDquark =Quark->GetPDGEncoding();
1021 << 1230 G4int AbsIDquark =std::abs(IDquark);
1022 G4int IDquark =Quark->GetPDGEncoding << 1231 G4int IDdi_quark =Di_Quark->GetPDGEncoding();
1023 G4int AbsIDquark =std::abs(IDquark); << 1232 G4int AbsIDdi_quark=std::abs(IDdi_quark);
1024 G4int IDdi_quark =Di_Quark->GetPDGEncodin << 1233 G4int Di_q1=AbsIDdi_quark/1000;
1025 G4int AbsIDdi_quark=std::abs(IDdi_quark); << 1234 G4int Di_q2=(AbsIDdi_quark-Di_q1*1000)/100;
1026 G4int Di_q1=AbsIDdi_quark/1000; << 1235
1027 G4int Di_q2=(AbsIDdi_quark-Di_q1*1000)/100; << 1236 G4int SignDiQ= 1;
1028 G4int SignDiQ= 1; << 1237 if(IDdi_quark < 0) SignDiQ=-1;
1029 if (IDdi_quark < 0) SignDiQ=-1; << 1238
1030 << 1239 NumberOf_FS=0;
1031 NumberOf_FS=0; << 1240 for(G4int ProdQ=1; ProdQ < 4; ProdQ++)
1032 for (G4int ProdQ=1; ProdQ < 4; ProdQ++) // << 1241 {
1033 { // << 1242 G4int SignQ;
1034 G4int SignQ; << 1243 if(IDquark > 0) {
1035 if (IDquark > 0) << 1244 SignQ=-1;
1036 { << 1245 if(IDquark == 2) SignQ= 1;
1037 SignQ=-1; << 1246 if((IDquark == 1) && (ProdQ == 3)) SignQ= 1; // K0
1038 if (IDquark == 2) << 1247 if((IDquark == 3) && (ProdQ == 1)) SignQ=-1; // K0bar
1039 if ((IDquark == 1) && (ProdQ == 3)) Sig << 1248 } else {
1040 if ((IDquark == 3) && (ProdQ == 1)) Sig << 1249 SignQ= 1;
1041 if (IDquark == 4) << 1250 if(IDquark == -2) SignQ=-1;
1042 if (IDquark == 5) Sig << 1251 if((IDquark ==-1) && (ProdQ == 3)) SignQ=-1; // K0bar
1043 } else << 1252 if((IDquark ==-3) && (ProdQ == 1)) SignQ= 1; // K0
1044 { << 1253 }
1045 SignQ= 1; <<
1046 if (IDquark == -2) Sig <<
1047 if ((IDquark ==-1) && (ProdQ == 3)) Sig <<
1048 if ((IDquark ==-3) && (ProdQ == 1)) Sig <<
1049 if (IDquark == -4) Sig <<
1050 if (IDquark == -5) Sig <<
1051 } <<
1052 <<
1053 if (AbsIDquark == ProdQ) SignQ <<
1054 <<
1055 G4int StateQ=0; <<
1056 const G4int maxNumberOfLoops <<
1057 G4int loopCounter = 0; <<
1058 do // while(Meson[AbsIDquark-1][ProdQ-1] <<
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 <<
1115 return true; <<
1116 } <<
1117 1254
1118 //------------------------------------------- << 1255 if(AbsIDquark == ProdQ) SignQ= 1;
1119 1256
1120 G4bool G4LundStringFragmentation::Quark_AntiQ << 1257 //G4cout<<G4endl;
1121 << 1258 //G4cout<<"Insert QQbar "<<ProdQ<<" Sign "<<SignQ<<G4endl;
1122 <<
1123 { <<
1124 G4double StringMass = string->Mass(); <<
1125 G4double StringMassSqr= sqr(StringMass); <<
1126 1259
1127 G4ParticleDefinition * Quark; << 1260 G4int StateQ=0;
1128 G4ParticleDefinition * Anti_Quark; << 1261 const G4int maxNumberOfLoops = 1000;
>> 1262 G4int loopCounter = 0;
>> 1263 do
>> 1264 {
>> 1265 LeftHadron=G4ParticleTable::GetParticleTable()->FindParticle(
>> 1266 SignQ*Meson[AbsIDquark-1][ProdQ-1][StateQ]);
>> 1267 G4double LeftHadronMass=LeftHadron->GetPDGMass();
>> 1268
>> 1269 G4int StateDiQ=0;
>> 1270 const G4int maxNumberOfInternalLoops = 1000;
>> 1271 G4int internalLoopCounter = 0;
>> 1272 do
>> 1273 {
>> 1274 RightHadron=G4ParticleTable::GetParticleTable()->FindParticle(
>> 1275 SignDiQ*Baryon[Di_q1-1][Di_q2-1][ProdQ-1][StateDiQ]);
>> 1276 G4double RightHadronMass=RightHadron->GetPDGMass();
1129 1277
1130 if (string->GetLeftParton()->GetPDGEncoding << 1278 if(StringMass > LeftHadronMass + RightHadronMass)
1131 { <<
1132 Quark =string->GetLeftParton(); <<
1133 Anti_Quark=string->GetRightParton(); <<
1134 } else <<
1135 { 1279 {
1136 Quark =string->GetRightParton(); << 1280 if ( NumberOf_FS > 34 ) {
1137 Anti_Quark=string->GetLeftParton(); << 1281 G4ExceptionDescription ed;
1138 } << 1282 ed << " NumberOf_FS exceeds its limit: NumberOf_FS=" << NumberOf_FS << G4endl;
1139 << 1283 G4Exception( "G4LundStringFragmentation::Quark_Diquark_lastSplitting ",
1140 G4int IDquark =Quark->GetPDGEncodin << 1284 "HAD_LUND_002", JustWarning, ed );
1141 G4int AbsIDquark =std::abs(IDquark); << 1285 NumberOf_FS = 34;
1142 G4int QuarkCharge =Qcharge[IDquar << 1286 }
1143 << 1287
1144 G4int IDanti_quark =Anti_Quark->GetPDGEn << 1288 G4double FS_Psqr=lambda(StringMassSqr,sqr(LeftHadronMass),sqr(RightHadronMass));
1145 G4int AbsIDanti_quark =std::abs(IDanti_quar << 1289 FS_Weight[NumberOf_FS]=std::sqrt(FS_Psqr)*MesonWeight[AbsIDquark-1][ProdQ-1][StateQ]*
1146 G4int AntiQuarkCharge =-Qcharge[AbsID << 1290 BaryonWeight[Di_q1-1][Di_q2-1][ProdQ-1][StateDiQ]*
1147 << 1291 Prob_QQbar[ProdQ-1];
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 1292
1168 if ((IDanti_quark ==-1) && (ProdQ == 3)) << 1293 FS_LeftHadron[NumberOf_FS] = LeftHadron;
1169 if ((IDanti_quark ==-3) && (ProdQ == 1)) << 1294 FS_RightHadron[NumberOf_FS]= RightHadron;
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 <<
1255 return true; <<
1256 } <<
1257 1295
1258 //------------------------------------------- << 1296 NumberOf_FS++;
1259 <<
1260 G4int G4LundStringFragmentation::SampleState( <<
1261 { <<
1262 if ( NumberOf_FS > 349 ) { <<
1263 G4ExceptionDescription ed; <<
1264 ed << " NumberOf_FS exceeds its lim <<
1265 G4Exception( "G4LundStringFragmenta <<
1266 NumberOf_FS = 349; <<
1267 } 1297 }
1268 1298
1269 G4double SumWeights=0.; << 1299 StateDiQ++;
1270 for (G4int i=0; i<NumberOf_FS; i++) {SumWei <<
1271 1300
1272 G4double ksi=G4UniformRand(); << 1301 } while( (Baryon[Di_q1-1][Di_q2-1][ProdQ-1][StateDiQ]!=0) &&
1273 G4double Sum=0.; << 1302 ++internalLoopCounter < maxNumberOfInternalLoops ); /* Loop checking, 07.08.2015, A.Ribon */
1274 G4int indexPosition = 0; << 1303 if ( internalLoopCounter >= maxNumberOfInternalLoops ) {
>> 1304 return false;
>> 1305 }
1275 1306
1276 for (G4int i=0; i<NumberOf_FS; i++) << 1307 StateQ++;
1277 { << 1308 } while( (Meson[AbsIDquark-1][ProdQ-1][StateQ]!=0) &&
1278 Sum+=(FS_Weight[i]/SumWeights); << 1309 ++loopCounter < maxNumberOfLoops ); /* Loop checking, 07.08.2015, A.Ribon */
1279 indexPosition=i; << 1310 if ( loopCounter >= maxNumberOfLoops ) {
1280 if (Sum >= ksi) break; << 1311 return false;
1281 } << 1312 }
1282 return indexPosition; << 1313 }
1283 } <<
1284 1314
1285 //------------------------------------------- << 1315 return true;
>> 1316 }
1286 1317
1287 void G4LundStringFragmentation::Sample4Moment << 1318 //----------------------------------------------------------------------------------------
1288 << 1319 G4bool G4LundStringFragmentation::
1289 << 1320 Quark_AntiQuark_lastSplitting(G4FragmentingString * & string,
>> 1321 G4ParticleDefinition * & LeftHadron,
>> 1322 G4ParticleDefinition * & RightHadron)
1290 { 1323 {
1291 // ------ Sampling of momenta of 2 last pro << 1324 G4double StringMass = string->Mass();
1292 G4ThreeVector Pt; << 1325 G4double StringMassSqr= sqr(StringMass);
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 1326
1300 G4double r_val = sqr(InitialMass*InitialMas << 1327 G4ParticleDefinition * Quark;
1301 sqr(2.*Mass*AntiMass); << 1328 G4ParticleDefinition * Anti_Quark;
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 1329
1318 G4int loopCounter = 0; << 1330 if(string->GetLeftParton()->GetPDGEncoding()>0) {
1319 do << 1331 Quark =string->GetLeftParton();
1320 { << 1332 Anti_Quark=string->GetRightParton();
1321 Pt=SampleQuarkPt(Pabs); Pt.setZ(0); G4dou << 1333 } else {
1322 MassMt = std::sqrt( Mass * Mass << 1334 Quark =string->GetRightParton();
1323 AntiMassMt= std::sqrt(AntiMass * AntiMass << 1335 Anti_Quark=string->GetLeftParton();
1324 } << 1336 }
1325 while ( (InitialMass < MassMt + AntiMassMt) << 1337
>> 1338 G4int IDquark =Quark->GetPDGEncoding();
>> 1339 G4int AbsIDquark =std::abs(IDquark);
>> 1340 G4int IDanti_quark =Anti_Quark->GetPDGEncoding();
>> 1341 G4int AbsIDanti_quark=std::abs(IDanti_quark);
>> 1342
>> 1343 NumberOf_FS=0;
>> 1344 for(G4int ProdQ=1; ProdQ < 4; ProdQ++)
>> 1345 {
>> 1346 G4int SignQ=-1;
>> 1347 if(IDquark == 2) SignQ= 1;
>> 1348 if((IDquark == 1) && (ProdQ == 3)) SignQ= 1; // K0
>> 1349 if((IDquark == 3) && (ProdQ == 1)) SignQ=-1; // K0bar
>> 1350 if(IDquark == ProdQ) SignQ= 1;
>> 1351
>> 1352 G4int SignAQ= 1;
>> 1353 if(IDanti_quark == -2) SignAQ=-1;
>> 1354 if((IDanti_quark ==-1) && (ProdQ == 3)) SignAQ=-1; // K0bar
>> 1355 if((IDanti_quark ==-3) && (ProdQ == 1)) SignAQ= 1; // K0
>> 1356 if(AbsIDanti_quark == ProdQ) SignAQ= 1;
>> 1357
>> 1358 G4int StateQ=0;
>> 1359 const G4int maxNumberOfLoops = 1000;
>> 1360 G4int loopCounter = 0;
>> 1361 do
>> 1362 {
>> 1363 LeftHadron=G4ParticleTable::GetParticleTable()->FindParticle(
>> 1364 SignQ*Meson[AbsIDquark-1][ProdQ-1][StateQ]);
>> 1365 G4double LeftHadronMass=LeftHadron->GetPDGMass();
>> 1366
>> 1367 G4int StateAQ=0;
>> 1368 const G4int maxNumberOfInternalLoops = 1000;
>> 1369 G4int internalLoopCounter = 0;
>> 1370 do
>> 1371 {
>> 1372 RightHadron=G4ParticleTable::GetParticleTable()->FindParticle(
>> 1373 SignAQ*Meson[AbsIDanti_quark-1][ProdQ-1][StateAQ]);
>> 1374 G4double RightHadronMass=RightHadron->GetPDGMass();
1326 1375
1327 SigmaQT = SigmaQTw; << 1376 if(StringMass > LeftHadronMass + RightHadronMass)
>> 1377 {
>> 1378 if ( NumberOf_FS > 34 ) {
>> 1379 G4ExceptionDescription ed;
>> 1380 ed << " NumberOf_FS exceeds its limit: NumberOf_FS=" << NumberOf_FS << G4endl;
>> 1381 G4Exception( "G4LundStringFragmentation::Quark_AntiQuark_lastSplitting ",
>> 1382 "HAD_LUND_003", JustWarning, ed );
>> 1383 NumberOf_FS = 34;
>> 1384 }
>> 1385
>> 1386 G4double FS_Psqr=lambda(StringMassSqr,sqr(LeftHadronMass),sqr(RightHadronMass));
>> 1387 //FS_Psqr=1.;
>> 1388 FS_Weight[NumberOf_FS]=std::sqrt(FS_Psqr)*MesonWeight[AbsIDquark-1][ProdQ-1][StateQ]*
>> 1389 MesonWeight[AbsIDanti_quark-1][ProdQ-1][StateAQ]*
>> 1390 Prob_QQbar[ProdQ-1];
>> 1391
>> 1392 if(string->GetLeftParton()->GetPDGEncoding()>0) {
>> 1393 FS_LeftHadron[NumberOf_FS] = RightHadron;
>> 1394 FS_RightHadron[NumberOf_FS]= LeftHadron;
>> 1395 } else {
>> 1396 FS_LeftHadron[NumberOf_FS] = LeftHadron;
>> 1397 FS_RightHadron[NumberOf_FS]= RightHadron;
>> 1398 }
>> 1399 NumberOf_FS++;
1328 1400
1329 AvailablePz2= sqr(InitialMass*InitialMass - << 1401 }
1330 4.*sqr(MassMt*AntiMassMt); <<
1331 1402
1332 AvailablePz2 /=(4.*InitialMass*InitialMass) << 1403 StateAQ++;
1333 AvailablePz = std::sqrt(AvailablePz2); << 1404 } while( (Meson[AbsIDanti_quark-1][ProdQ-1][StateAQ]!=0) &&
>> 1405 ++internalLoopCounter < maxNumberOfInternalLoops ); /* Loop checking, 07.08.2015, A.Ribon */
>> 1406 if ( internalLoopCounter >= maxNumberOfInternalLoops ) {
>> 1407 return false;
>> 1408 }
1334 1409
1335 G4double Px=Pt.getX(); << 1410 StateQ++;
1336 G4double Py=Pt.getY(); << 1411 } while( (Meson[AbsIDquark-1][ProdQ-1][StateQ]!=0) &&
>> 1412 ++loopCounter < maxNumberOfLoops ); /* Loop checking, 07.08.2015, A.Ribon */
>> 1413 if ( loopCounter >= maxNumberOfLoops ) {
>> 1414 return false;
>> 1415 }
>> 1416 }
>> 1417 return true;
>> 1418 }
1337 1419
1338 Mom->setPx(Px); Mom->setPy(Py); Mom->setPz( << 1420 //----------------------------------------------------------------------------------------------------------
1339 Mom->setE(std::sqrt(sqr(MassMt)+AvailablePz << 1421 G4int G4LundStringFragmentation::SampleState(void)
>> 1422 {
>> 1423 if ( NumberOf_FS > 34 ) {
>> 1424 G4ExceptionDescription ed;
>> 1425 ed << " NumberOf_FS exceeds its limit: NumberOf_FS=" << NumberOf_FS << G4endl;
>> 1426 G4Exception( "G4LundStringFragmentation::SampleState ", "HAD_LUND_004", JustWarning, ed );
>> 1427 NumberOf_FS = 34;
>> 1428 }
>> 1429
>> 1430 G4double SumWeights=0.;
1340 1431
1341 AntiMom->setPx(-Px); AntiMom->setPy(-Py); A << 1432 for(G4int i=0; i<NumberOf_FS; i++) {SumWeights+=FS_Weight[i];}// G4cout<<i<<" "<<FS_Weight[i]<<G4endl;}
1342 AntiMom->setE (std::sqrt(sqr(AntiMassMt)+Av <<
1343 1433
1344 #ifdef debug_LUNDfragmentation << 1434 G4double ksi=G4UniformRand();
1345 G4cout<<"Fmass Mom "<<Mom->getX()<<" << 1435 G4double Sum=0.;
1346 G4cout<<"Smass Mom "<<AntiMom->getX() << 1436 G4int indexPosition = 0;
1347 <<" "<<AntiMom->getT()<<G4endl; << 1437
1348 #endif << 1438 for(G4int i=0; i<NumberOf_FS; i++)
>> 1439 {
>> 1440 Sum+=(FS_Weight[i]/SumWeights);
>> 1441 indexPosition=i;
>> 1442 if(Sum >= ksi) break;
>> 1443 }
>> 1444 return indexPosition;
1349 } 1445 }
1350 1446
1351 //------------------------------------------- << 1447 // Uzhi June 2014 Insert from G4ExcitedStringDecay.cc
>> 1448 //-----------------------------------------------------------------------------
1352 1449
1353 G4double G4LundStringFragmentation::lambda(G4 << 1450 G4ParticleDefinition * G4LundStringFragmentation::
1354 { << 1451 DiQuarkSplitup( G4ParticleDefinition* decay, G4ParticleDefinition *&created )
1355 G4double lam = sqr(S - m1_Sqr - m2_Sqr) - 4 << 1452 {
1356 return lam; << 1453 //... can Diquark break or not?
>> 1454 if (G4UniformRand() < DiquarkBreakProb )
>> 1455 {
>> 1456 //... Diquark break
>> 1457
>> 1458 G4int stableQuarkEncoding = decay->GetPDGEncoding()/1000;
>> 1459 G4int decayQuarkEncoding = (decay->GetPDGEncoding()/100)%10;
>> 1460 if (G4UniformRand() < 0.5)
>> 1461 {
>> 1462 G4int Swap = stableQuarkEncoding;
>> 1463 stableQuarkEncoding = decayQuarkEncoding;
>> 1464 decayQuarkEncoding = Swap;
>> 1465 }
>> 1466
>> 1467 G4int IsParticle=(decayQuarkEncoding>0) ? -1 : +1; // if we have a quark, we need antiquark)
>> 1468
>> 1469 //G4cout<<"GetStrangeSuppress() "<<GetStrangeSuppress()<<G4endl;
>> 1470 //G4int Uzhi; G4cin>>Uzhi;
>> 1471
>> 1472 //G4double StrSup=GetStrangeSuppress();
>> 1473 //StrangeSuppress=0.34;
>> 1474 pDefPair QuarkPair = CreatePartonPair(IsParticle,false); // no diquarks wanted
>> 1475 //StrangeSuppress=StrSup;
>> 1476 //... Build new Diquark
>> 1477 G4int QuarkEncoding=QuarkPair.second->GetPDGEncoding();
>> 1478 G4int i10 = std::max(std::abs(QuarkEncoding), std::abs(stableQuarkEncoding));
>> 1479 G4int i20 = std::min(std::abs(QuarkEncoding), std::abs(stableQuarkEncoding));
>> 1480 G4int spin = (i10 != i20 && G4UniformRand() <= 0.5)? 1 : 3;
>> 1481 G4int NewDecayEncoding = -1*IsParticle*(i10 * 1000 + i20 * 100 + spin);
>> 1482 created = FindParticle(NewDecayEncoding);
>> 1483 G4ParticleDefinition * decayQuark=FindParticle(decayQuarkEncoding);
>> 1484 G4ParticleDefinition * had=hadronizer->Build(QuarkPair.first, decayQuark);
>> 1485 return had;
>> 1486 // return hadronizer->Build(QuarkPair.first, decayQuark);
>> 1487
>> 1488 } else {
>> 1489 //... Diquark does not break
>> 1490
>> 1491 G4int IsParticle=(decay->GetPDGEncoding()>0) ? +1 : -1; // if we have a diquark, we need quark)
>> 1492
>> 1493 G4double StrSup=GetStrangeSuppress();
>> 1494 StrangeSuppress=0.43; //0.42 0.38
>> 1495 pDefPair QuarkPair = CreatePartonPair(IsParticle,false); // no diquarks wanted
>> 1496 StrangeSuppress=StrSup;
>> 1497
>> 1498 created = QuarkPair.second;
>> 1499
>> 1500 G4ParticleDefinition * had=hadronizer->Build(QuarkPair.first, decay);
>> 1501 return had;
>> 1502 }
1357 } 1503 }
1358 1504
1359 // ------------------------------------------ << 1505 // Uzhi June 2014 End of the inserting
1360 G4LundStringFragmentation::~G4LundStringFragm <<
1361 {} <<
1362 1506
1363 1507