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