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