| 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 // 27 //
28 // ------------------------------------------- 28 // -----------------------------------------------------------------------------
29 // GEANT 4 class implementation file 29 // GEANT 4 class implementation file
30 // 30 //
31 // History: first implementation, Maxim K 31 // History: first implementation, Maxim Komogorov, 10-Jul-1998
32 // ------------------------------------------- 32 // -----------------------------------------------------------------------------
33 #include "G4QGSMFragmentation.hh" 33 #include "G4QGSMFragmentation.hh"
34 #include "G4PhysicalConstants.hh" 34 #include "G4PhysicalConstants.hh"
35 #include "G4SystemOfUnits.hh" 35 #include "G4SystemOfUnits.hh"
36 #include "Randomize.hh" 36 #include "Randomize.hh"
37 #include "G4ios.hh" 37 #include "G4ios.hh"
38 #include "G4FragmentingString.hh" 38 #include "G4FragmentingString.hh"
39 #include "G4DiQuarks.hh" 39 #include "G4DiQuarks.hh"
40 #include "G4Quarks.hh" 40 #include "G4Quarks.hh"
41 #include "G4HadronicParameters.hh" << 41
42 #include "G4Pow.hh" 42 #include "G4Pow.hh"
43 43
44 //#define debug_QGSMfragmentation 44 //#define debug_QGSMfragmentation
45 45
46 // Class G4QGSMFragmentation 46 // Class G4QGSMFragmentation
47 //******************************************** 47 //****************************************************************************************
48 48
49 G4QGSMFragmentation::G4QGSMFragmentation() << 49 G4QGSMFragmentation::G4QGSMFragmentation() :
>> 50 arho(0.5), aphi(0.), an(-0.5), ala(-0.75), aksi(-1.), alft(0.5)
50 { 51 {
51 SigmaQT = 0.45 * GeV; << 52 SetStrangenessSuppression((1.0 - 0.16)/2.);
>> 53 SetDiquarkSuppression(0.299);
>> 54 SetDiquarkBreakProbability(0.7);
52 55
53 MassCut = 0.35*GeV; << 56 //... pspin_meson is probability to create pseudo-scalar meson
>> 57 pspin_meson = 0.25; SetVectorMesonProbability(pspin_meson);
54 58
55 SetStrangenessSuppression((1.0 - 0.16)/2.) << 59 //... pspin_barion is probability to create 1/2 barion
>> 60 pspin_barion = 0.5; SetSpinThreeHalfBarionProbability(pspin_barion);
56 61
57 // Check if charmed and bottom hadrons are << 62 //... vectorMesonMix[] is quark mixing parameters for vector mesons (Variable spin = 3)
58 // set the non-zero probabilities for c-cb << 63 vectorMesonMix[0] = 0.; //0.5
59 // else set them to 0.0. If these probabil << 64 vectorMesonMix[1] = 0.375; //0.5;
60 // hadrons can't/can be created during the << 65 vectorMesonMix[2] = 0.0;
61 // (i.e. not heavy) projectile hadron nucl << 66 vectorMesonMix[3] = 0.375; //0.5;
62 if ( G4HadronicParameters::Instance()->Ena << 67 vectorMesonMix[4] = 1.0;
63 SetProbCCbar(0.0002); // According to O << 68 vectorMesonMix[5] = 1.0;
64 SetProbBBbar(5.0e-5); // According to O <<
65 } else { <<
66 SetProbCCbar(0.0); <<
67 SetProbBBbar(0.0); <<
68 } <<
69 <<
70 SetDiquarkSuppression(0.32); <<
71 SetDiquarkBreakProbability(0.7); <<
72 69
>> 70 SetVectorMesonMixings(vectorMesonMix);
73 SetMinMasses(); 71 SetMinMasses();
74 <<
75 arho = 0.5; // alpha_rho0 <<
76 aphi = 0.0; // alpha_fi <<
77 aJPs =-2.2; // alpha_J/Psi <<
78 aUps =-8.0; // alpha_Y ??? O. Pisk <<
79 <<
80 aksi =-1.0; <<
81 alft = 0.5; // 2 * alpha'_R *<Pt^2> <<
82 <<
83 an = -0.5 ; <<
84 ala = -0.75; // an - arho/2 + aphi/2 <<
85 alaC = an - arho/2.0 + aJPs/2.0; <<
86 alaB = an - arho/2.0 + aUps/2.0; <<
87 aXi = 0.0; // ?? <<
88 aXiC = 0.0; // ?? <<
89 aXiB = 0.0; // ?? <<
90 aXiCC = 0.0; // ?? <<
91 aXiCB = 0.0; // ?? <<
92 aXiBB = 0.0; // ?? <<
93 <<
94 SetFFq2q(); <<
95 SetFFq2qq(); <<
96 SetFFqq2q(); <<
97 SetFFqq2qq(); <<
98 // d u s c b <<
99 G4int Index[5][5] = { { 0, 1, 2, 3, 4 } <<
100 { 1, 5, 6, 7, 8 } <<
101 { 2, 6, 9, 10, 11 } <<
102 { 3, 7, 10, 12, 13 } <<
103 { 4, 8, 11, 13, 14 } <<
104 for (G4int i = 0; i < 5; i++ ) { <<
105 for ( G4int j = 0; j < 5; j++ ) { <<
106 IndexDiQ[i][j] = Index[i][j]; <<
107 } <<
108 }; <<
109 } 72 }
110 73
111 G4QGSMFragmentation::~G4QGSMFragmentation() 74 G4QGSMFragmentation::~G4QGSMFragmentation()
112 {} 75 {}
113 76
114 //-------------------------------------------- 77 //----------------------------------------------------------------------------------------------------------
115 78
116 G4KineticTrackVector* G4QGSMFragmentation::Fra 79 G4KineticTrackVector* G4QGSMFragmentation::FragmentString(const G4ExcitedString& theString)
117 { 80 {
118 <<
119 G4FragmentingString aString(theString); <<
120 SetMinimalStringMass(&aString); <<
121 <<
122 #ifdef debug_QGSMfragmentation 81 #ifdef debug_QGSMfragmentation
123 G4cout<<G4endl<<"QGSM StringFragm: String Ma 82 G4cout<<G4endl<<"QGSM StringFragm: String Mass "
124 <<theString.Get4M 83 <<theString.Get4Momentum().mag()<<" Pz "
125 <<theString.Get4M 84 <<theString.Get4Momentum().pz()
126 <<"-------------- 85 <<"------------------------------------"<<G4endl;
127 G4cout<<"String ends Direct "<<theString.Get 86 G4cout<<"String ends Direct "<<theString.GetLeftParton()->GetPDGcode()<<" "
128 <<theString.Get 87 <<theString.GetRightParton()->GetPDGcode()<<" "
129 <<theString.Get 88 <<theString.GetDirection()<< G4endl;
130 G4cout<<"Left mom "<<theString.GetLeftParto 89 G4cout<<"Left mom "<<theString.GetLeftParton()->Get4Momentum()<<G4endl;
131 G4cout<<"Right mom "<<theString.GetRightPart 90 G4cout<<"Right mom "<<theString.GetRightParton()->Get4Momentum()<<G4endl;
132 G4cout<<"Check for Fragmentation "<<G4endl; 91 G4cout<<"Check for Fragmentation "<<G4endl;
133 #endif 92 #endif
134 93
135 // Can no longer modify Parameters for Fragm 94 // Can no longer modify Parameters for Fragmentation.
136 PastInitPhase=true; 95 PastInitPhase=true;
137 96
138 // Check if string has enough mass to fragme 97 // Check if string has enough mass to fragment...
139 G4KineticTrackVector * LeftVector=NULL; << 98 G4KineticTrackVector * LeftVector=LightFragmentationTest(&theString);
140 <<
141 if ( !IsItFragmentable(&aString) ) { <<
142 LeftVector=ProduceOneHadron(&theString); <<
143 99
144 #ifdef debug_QGSMfragmentation << 100 #ifdef debug_QGSMfragmentation
145 if ( LeftVector != 0 ) G4cout<<"Non fragm << 101 if ( LeftVector != 0 ) G4cout<<"Non fragmentable - the string is converted to one hadron "<<G4endl;
146 #endif << 102 #endif
147 103
148 if ( LeftVector == nullptr ) LeftVector = << 104 if ( LeftVector != 0 ) return LeftVector;
149 return LeftVector; <<
150 } <<
151 105
152 #ifdef debug_QGSMfragmentation 106 #ifdef debug_QGSMfragmentation
153 G4cout<<"The string will be fragmented. "<<G 107 G4cout<<"The string will be fragmented. "<<G4endl;
154 #endif 108 #endif
155 109
156 LeftVector = new G4KineticTrackVector; 110 LeftVector = new G4KineticTrackVector;
157 G4KineticTrackVector * RightVector=new G4Kin 111 G4KineticTrackVector * RightVector=new G4KineticTrackVector;
158 112
159 G4ExcitedString *theStringInCMS=CopyExcited( 113 G4ExcitedString *theStringInCMS=CopyExcited(theString);
160 G4LorentzRotation toCms=theStringInCMS->Tran 114 G4LorentzRotation toCms=theStringInCMS->TransformToAlignedCms();
161 115
162 G4bool success=false, inner_sucess=true; 116 G4bool success=false, inner_sucess=true;
163 G4int attempt=0; 117 G4int attempt=0;
164 while ( !success && attempt++ < StringLoopIn 118 while ( !success && attempt++ < StringLoopInterrupt ) /* Loop checking, 07.08.2015, A.Ribon */
165 { 119 {
166 #ifdef debug_QGSMfragmentation 120 #ifdef debug_QGSMfragmentation
167 G4cout<<"Loop_toFrag "<<theStr 121 G4cout<<"Loop_toFrag "<<theStringInCMS->GetLeftParton()->GetPDGcode()<<" "
168 <<theStr 122 <<theStringInCMS->GetRightParton()->GetPDGcode()<<" "
169 <<theStr 123 <<theStringInCMS->GetDirection()<< G4endl;
170 #endif 124 #endif
171 125
172 G4FragmentingString *currentString=new G4F 126 G4FragmentingString *currentString=new G4FragmentingString(*theStringInCMS);
173 127
174 std::for_each(LeftVector->begin(), LeftVec 128 std::for_each(LeftVector->begin(), LeftVector->end(), DeleteKineticTrack());
175 LeftVector->clear(); 129 LeftVector->clear();
176 std::for_each(RightVector->begin(), RightV 130 std::for_each(RightVector->begin(), RightVector->end(), DeleteKineticTrack());
177 RightVector->clear(); 131 RightVector->clear();
178 132
179 inner_sucess=true; // set false on failur 133 inner_sucess=true; // set false on failure..
180 const G4int maxNumberOfLoops = 134 const G4int maxNumberOfLoops = 1000;
181 G4int loopCounter = -1; 135 G4int loopCounter = -1;
182 while (! StopFragmenting(currentString) && 136 while (! StopFragmenting(currentString) && ++loopCounter < maxNumberOfLoops ) /* Loop checking, 07.08.2015, A.Ribon */
183 { // Split current string into hadron + n 137 { // Split current string into hadron + new string
184 138
185 #ifdef debug_QGSMfragm 139 #ifdef debug_QGSMfragmentation
186 G4cout<<"The string ca 140 G4cout<<"The string can fragment. "<<G4endl;;
187 #endif 141 #endif
188 G4FragmentingString *newString=0; // us 142 G4FragmentingString *newString=0; // used as output from SplitUp...
189 G4KineticTrack * Hadron=Splitup(currentS 143 G4KineticTrack * Hadron=Splitup(currentString,newString);
190 144
191 if ( Hadron != 0 ) 145 if ( Hadron != 0 )
192 { 146 {
193 #ifdef debug_QGSMfr 147 #ifdef debug_QGSMfragmentation
194 G4cout<<"Hadron pro 148 G4cout<<"Hadron prod at fragm. "<<Hadron->GetDefinition()->GetParticleName()<<G4endl;
195 #endif 149 #endif
196 // To close the pro 150 // To close the production of hadrons at fragmentation stage
197 if ( currentString->GetDecayDirection 151 if ( currentString->GetDecayDirection() > 0 )
198 LeftVector->push_back(Hadron); 152 LeftVector->push_back(Hadron);
199 else 153 else
200 RightVector->push_back(Hadron); 154 RightVector->push_back(Hadron);
201 155
202 delete currentString; 156 delete currentString;
203 currentString=newString; 157 currentString=newString;
204 158
205 } else { 159 } else {
206 160
207 #ifdef debug_QGSMfr 161 #ifdef debug_QGSMfragmentation
208 G4cout<<"abandon .. 162 G4cout<<"abandon ... start from the beginning ---------------"<<G4endl;
209 #endif 163 #endif
210 164
211 // Abandon ... start from the beginni 165 // Abandon ... start from the beginning
212 if (newString) delete newString; 166 if (newString) delete newString;
213 inner_sucess=false; 167 inner_sucess=false;
214 break; 168 break;
215 } 169 }
216 } 170 }
217 if ( loopCounter >= maxNumberO 171 if ( loopCounter >= maxNumberOfLoops ) {
218 inner_sucess=false; 172 inner_sucess=false;
219 } 173 }
220 174
221 // Split current string into 2 final Hadro 175 // Split current string into 2 final Hadrons
222 #ifdef debug_QGSMfragmentation 176 #ifdef debug_QGSMfragmentation
223 if( inner_sucess ) { << 177 G4cout<<"Split remaining string into 2 final hadrons."<<G4endl;
224 G4cout<<"Split remaining str <<
225 } else { <<
226 G4cout<<" New attempt to fragment string <<
227 } <<
228 #endif 178 #endif
229 // To the close production of 179 // To the close production of hadrons at last string decay
230 if ( inner_sucess && 180 if ( inner_sucess &&
231 SplitLast(currentString,LeftVector, R 181 SplitLast(currentString,LeftVector, RightVector) )
232 { 182 {
233 success=true; 183 success=true;
234 } 184 }
235 delete currentString; 185 delete currentString;
236 } 186 }
237 187
238 delete theStringInCMS; 188 delete theStringInCMS;
239 189
240 if ( ! success ) 190 if ( ! success )
241 { 191 {
242 std::for_each(LeftVector->begin(), LeftVec 192 std::for_each(LeftVector->begin(), LeftVector->end(), DeleteKineticTrack());
243 LeftVector->clear(); 193 LeftVector->clear();
244 std::for_each(RightVector->begin(), RightV 194 std::for_each(RightVector->begin(), RightVector->end(), DeleteKineticTrack());
245 delete RightVector; 195 delete RightVector;
246 return LeftVector; 196 return LeftVector;
247 } 197 }
248 198
249 // Join Left- and RightVector into LeftVecto 199 // Join Left- and RightVector into LeftVector in correct order.
250 while(!RightVector->empty()) /* Loop checki 200 while(!RightVector->empty()) /* Loop checking, 07.08.2015, A.Ribon */
251 { 201 {
252 LeftVector->push_back(RightVector->back( 202 LeftVector->push_back(RightVector->back());
253 RightVector->erase(RightVector->end()-1) 203 RightVector->erase(RightVector->end()-1);
254 } 204 }
255 delete RightVector; 205 delete RightVector;
256 206
257 CalculateHadronTimePosition(theString.Get4Mo 207 CalculateHadronTimePosition(theString.Get4Momentum().mag(), LeftVector);
258 208
259 G4LorentzRotation toObserverFrame(toCms.inve 209 G4LorentzRotation toObserverFrame(toCms.inverse());
260 210
261 for (size_t C1 = 0; C1 < LeftVector->size(); << 211 for(size_t C1 = 0; C1 < LeftVector->size(); C1++)
262 { 212 {
263 G4KineticTrack* Hadron = LeftVector->oper 213 G4KineticTrack* Hadron = LeftVector->operator[](C1);
264 G4LorentzVector Momentum = Hadron->Get4Mo 214 G4LorentzVector Momentum = Hadron->Get4Momentum();
265 Momentum = toObserverFrame*Momentum; 215 Momentum = toObserverFrame*Momentum;
266 Hadron->Set4Momentum(Momentum); 216 Hadron->Set4Momentum(Momentum);
267 G4LorentzVector Coordinate(Hadron->GetPos 217 G4LorentzVector Coordinate(Hadron->GetPosition(), Hadron->GetFormationTime());
268 Momentum = toObserverFrame*Coordinate; 218 Momentum = toObserverFrame*Coordinate;
269 Hadron->SetFormationTime(Momentum.e()); 219 Hadron->SetFormationTime(Momentum.e());
270 G4ThreeVector aPosition(Momentum.vect()); 220 G4ThreeVector aPosition(Momentum.vect());
271 Hadron->SetPosition(theString.GetPosition 221 Hadron->SetPosition(theString.GetPosition()+aPosition);
272 } 222 }
273 return LeftVector; 223 return LeftVector;
274 } 224 }
275 225
276 //-------------------------------------------- 226 //----------------------------------------------------------------------------------------------------------
277 227
278 G4bool G4QGSMFragmentation::IsItFragmentable(c << 228 G4double G4QGSMFragmentation::GetLightConeZ(G4double zmin, G4double zmax, G4int PartonEncoding,
279 { << 229 G4ParticleDefinition* pHadron, G4double , G4double )
280 return sqr( MinimalStringMass + MassCut ) < << 230 {
281 } << 231 #ifdef debug_QGSMfragmentation
282 << 232 G4cout<<"GetLightConeZ zmin zmax Parton pHadron "<<zmin<<" "<<zmax<<" "<<PartonEncoding<<" "<<pHadron->GetParticleName()<<G4endl;
283 //-------------------------------------------- << 233 #endif
284 <<
285 G4bool G4QGSMFragmentation::StopFragmenting(co <<
286 { <<
287 SetMinimalStringMass(string); <<
288 if ( MinimalStringMass < 0.0 ) return <<
289 <<
290 G4double smass = string->Mass(); <<
291 G4double x = (string->IsAFourQuarkString()) <<
292 : 0.66e-6*(smass - MinimalStringMass)*(sma <<
293 <<
294 G4bool res = true; <<
295 if(x > 0.0) { <<
296 res = (x < 200.) ? (G4UniformRand() <<
297 } <<
298 return res; <<
299 } <<
300 <<
301 //-------------------------------------------- <<
302 <<
303 G4KineticTrack * G4QGSMFragmentation::Splitup( <<
304 G4FragmentingStri <<
305 { <<
306 #ifdef debug_QGSMfragmentation <<
307 G4cout<<G4endl; <<
308 G4cout<<"Start SplitUP (G4VLongitudinal <<
309 G4cout<<"String partons: " <<string->Ge <<
310 <<string->Ge <<
311 <<"Direction " <<string->Ge <<
312 #endif <<
313 <<
314 //... random choice of string end to us <<
315 G4int SideOfDecay = (G4UniformRand() < <<
316 if (SideOfDecay < 0) <<
317 { <<
318 string->SetLeftPartonStable(); <<
319 } else <<
320 { <<
321 string->SetRightPartonStable(); <<
322 } <<
323 <<
324 G4ParticleDefinition *newStringEnd; <<
325 G4ParticleDefinition * HadronDefinition <<
326 if (string->DecayIsQuark()) <<
327 { <<
328 G4double ProbDqADq = GetDiquarkSuppress(); <<
329 <<
330 G4int NumberOfpossibleBaryons = 2; <<
331 <<
332 if (string->GetLeftParton()->GetParticleSu <<
333 if (string->GetRightParton()->GetParticleS <<
334 <<
335 G4double ActualProb = ProbDqADq ; <<
336 ActualProb *= (1.0-G4Exp(2.0*(1.0 - <<
337 <<
338 SetDiquarkSuppression(ActualProb); <<
339 <<
340 HadronDefinition= QuarkSplitup(strin <<
341 <<
342 SetDiquarkSuppression(ProbDqADq); <<
343 } else { <<
344 HadronDefinition= DiQuarkSplitup(str <<
345 } <<
346 <<
347 if ( HadronDefinition == NULL ) return <<
348 <<
349 #ifdef debug_QGSMfragmentation <<
350 G4cout<<"The parton "<<string->GetDecay <<
351 <<" produces hadron "<<HadronDefi <<
352 <<" and is transformed to "<<newS <<
353 G4cout<<"The side of the string decay L <<
354 #endif <<
355 // create new String from old, ie. keep <<
356 234
357 newString=new G4FragmentingString(*stri << 235 G4double z;
358 << 236 G4double d1, d2, yf;
>> 237 G4double invD1(0.),invD2(0.), r1(0.),r2(0.),r12(0.);
359 238
360 #ifdef debug_QGSMfragmentation << 239 G4int absCode = std::abs( PartonEncoding );
361 G4cout<<"An attempt to determine its en << 240 G4int absHadronCode=std::abs(pHadron->GetPDGEncoding());
362 #endif <<
363 G4LorentzVector* HadronMomentum=SplitEa <<
364 241
365 delete newString; newString=0; << 242 G4int q1, q2, q3;
366 << 243 q1 = absHadronCode/1000; q2 = (absHadronCode % 1000)/100; q3 = (absHadronCode % 100)/10;
367 G4KineticTrack * Hadron =0; <<
368 if ( HadronMomentum != 0 ) { <<
369 244
370 #ifdef debug_QGSMfragmentation << 245 G4bool StrangeHadron = (q1 == 3) || (q2 == 3) || (q3 == 3);
371 G4cout<<"The attempt was successful <<
372 #endif <<
373 G4ThreeVector Pos; <<
374 Hadron = new G4KineticTrack(HadronDefinit <<
375 246
376 newString=new G4FragmentingString(*st << 247 if (absCode < 10)
>> 248 { // A quark fragmentation ----------------------------
>> 249 if (absCode == 1 || absCode == 2)
>> 250 {
>> 251 if (absHadronCode < 1000)
>> 252 { // Meson produced
>> 253 if ( !StrangeHadron ) {d1=2.0; d2 = -arho + alft;} // d1=2.0;
>> 254 else {d1=1.0; d2 = -aphi + alft;} // d1=1.0;
>> 255 }
>> 256 else
>> 257 { // Baryon produced
>> 258 if ( !StrangeHadron ) {d1=0.0; d2 = arho - 2.0*an + alft;}
>> 259 else {d1=0.0; d2 = 2.0*arho - 2.0*an - aphi + alft;}
>> 260 }
>> 261 }
>> 262 else if (absCode == 3)
>> 263 {
>> 264 if (absHadronCode < 1000) {d1=1.0 - aphi; d2 = -arho + alft;} // Meson produced s->K + u/d d1=1.0
>> 265 else {d1=1.0 - aphi; d2 = arho - 2.0*an + alft;} // Baryon produced
377 266
378 delete HadronMomentum; << 267 } else throw G4HadronicException(__FILE__, __LINE__, "Unknown PDGencoding in G4QGSMFragmentation::G4LightConeZ");
379 } <<
380 else <<
381 { <<
382 #ifdef debug_QGSMfragmentation <<
383 G4cout<<"The attempt was not successf <<
384 #endif <<
385 } <<
386 268
387 #ifdef debug_VStringDecay << 269 #ifdef debug_QGSMfragmentation
388 G4cout<<"End SplitUP (G4VLongitudinalSt << 270 G4cout<<"d1 d2 "<<d1<<" "<<d2<<G4endl;
389 #endif << 271 #endif
390 272
391 return Hadron; << 273 d1+=1.0; d2+=1.0;
392 } <<
393 274
394 //-------------------------------------------- << 275 invD1=1./d1; invD2=1./d2;
395 276
396 G4ParticleDefinition *G4QGSMFragmentation::DiQ << 277 const G4int maxNumberOfLoops = 10000;
397 << 278 G4int loopCounter = 0;
398 { << 279 do
399 //... can Diquark break or not? << 280 {
400 if (G4UniformRand() < DiquarkBreakProb ) / << 281 r1=G4Pow::GetInstance()->powA(G4UniformRand(),invD1);
401 { << 282 r2=G4Pow::GetInstance()->powA(G4UniformRand(),invD2);
402 G4int stableQuarkEncoding = decay->GetPD << 283 r12=r1+r2;
403 G4int decayQuarkEncoding = (decay->GetPD << 284 z=r1/r12;
>> 285 } while ( ( (r12 > 1.0) || !((zmin <= z)&&(z <= zmax))) && ++loopCounter < maxNumberOfLoops ); /* Loop checking, 07.08.2015, A.Ribon */
>> 286 if ( loopCounter >= maxNumberOfLoops ) {
>> 287 z = 0.5*(zmin + zmax); // Just a value between zmin and zmax, no physics considerations at all!
>> 288 }
404 289
405 if (G4UniformRand() < 0.5) << 290 return z;
>> 291 }
>> 292 else
>> 293 { // A di-quark fragmentation -------------------------
>> 294 if (absCode == 1103 || absCode == 2101 ||
>> 295 absCode == 2203 || absCode == 2103)
>> 296 {
>> 297 if(absHadronCode < 1000) // Meson production
406 { 298 {
407 G4int Swap = stableQuarkEncoding; << 299 if ( !StrangeHadron ) {d1=1.0; d2= arho - 2.0*an + alft;}
408 stableQuarkEncoding = decayQuarkEncod << 300 else {d1=0.0; d2 = 2.*arho - 2.0*an - aphi + alft;} // d1=1.0;
409 decayQuarkEncoding = Swap; << 301 }
>> 302 else // Baryon production
>> 303 {
>> 304 if ( !StrangeHadron ) {d1=2.0*(arho - an); d2= -arho + alft;}
>> 305 else {d1=2.0*(arho - an); d2 =-aphi + alft;}
410 } 306 }
411 307
412 G4int IsParticle=(decayQuarkEncoding>0) << 308 #ifdef debug_QGSMfragmentation
413 << 309 G4cout<<"d1 d2 "<<d1<<" "<<d2<<G4endl;
414 G4double StrSup=GetStrangeSuppress(); << 310 #endif
415 SetStrangenessSuppression((1.0 - 0.07)/2 << 311
416 pDefPair QuarkPair = CreatePartonPair(Is << 312 d1+=1.0; d2+=1.0;
417 SetStrangenessSuppression(StrSup); << 313 invD1=1./d1; invD2=1./d2;
418 << 314
419 //... Build new Diquark << 315 const G4int maxNumberOfLoops = 10000;
420 G4int QuarkEncoding=QuarkPair.second->Ge << 316 G4int loopCounter = 0;
421 G4int i10 = std::max(std::abs(QuarkEnco << 317 do
422 G4int i20 = std::min(std::abs(QuarkEnco << 318 {
423 G4int spin = (i10 != i20 && G4UniformRan << 319 r1=G4Pow::GetInstance()->powA(G4UniformRand(),invD1);
424 G4int NewDecayEncoding = -1*IsParticle*( << 320 r2=G4Pow::GetInstance()->powA(G4UniformRand(),invD2);
425 << 321 r12=r1+r2;
426 created = FindParticle(NewDecayEncoding) << 322 z=r1/r12;
427 G4ParticleDefinition * decayQuark=FindPa << 323 } while ( ( (r12 > 1.0) || !((zmin <= z)&&(z <= zmax))) && ++loopCounter < maxNumberOfLoops ); /* Loop checking, 07.08.2015, A.Ribon */
428 G4ParticleDefinition * had=hadronizer->B << 324 if ( loopCounter >= maxNumberOfLoops ) {
429 << 325 z = 0.5*(zmin + zmax); // Just a value between zmin and zmax, no physics considerations at all!
430 DecayQuark = decay->GetPDGEncoding(); << 326 }
431 NewQuark = NewDecayEncoding; <<
432 <<
433 return had; <<
434 327
435 } else { //... Diquark does not break << 328 return z;
>> 329 }
>> 330 else if (absCode == 3101 || absCode == 3103 || // For strange d-quarks
>> 331 absCode == 3201 || absCode == 3203)
>> 332 {
>> 333 // For future improvements
>> 334 // if (absHadronCode < 1000) {d1=1.0;} // Meson production
>> 335 // else {d1=2.0;} // Baryon production
>> 336 d2 = (alft - (2.*ala - arho));
>> 337 }
>> 338 else
>> 339 {
>> 340 // if (absHadronCode < 1000) {d1=1.0;} // Meson production
>> 341 // else {d1=2.0;} // Baryon production
>> 342 d2 = (alft - (2.*aksi - arho));
>> 343 }
436 344
437 G4int IsParticle=(decay->GetPDGEncoding( << 345 const G4int maxNumberOfLoops = 1000;
438 G4double StrSup=GetStrangeSuppress(); / << 346 G4int loopCounter = 0;
439 SetStrangenessSuppression((1.0 - 0.07)/2 << 347 do
440 pDefPair QuarkPair = CreatePartonPair(Is << 348 {
441 SetStrangenessSuppression(StrSup); << 349 z = zmin + G4UniformRand() * (zmax - zmin);
>> 350 d1 = (1. - z);
>> 351 yf = G4Pow::GetInstance()->powA(d1, d2);
>> 352 }
>> 353 while ( (G4UniformRand() > yf) && ++loopCounter < maxNumberOfLoops ); // Loop checking, 07.08.2015, A.Ribon //
>> 354 /* For future improvements
>> 355 d1+=1.0; d2+=1.0;
>> 356 invD1=1./d1; invD2=1./d2;
>> 357 const G4int maxNumberOfLoops = 10000;
>> 358 G4int loopCounter = 0;
>> 359 do
>> 360 {
>> 361 r1=G4Pow::GetInstance()->powA(G4UniformRand(),invD1);
>> 362 r2=G4Pow::GetInstance()->powA(G4UniformRand(),invD2);
>> 363 r12=r1+r2;
>> 364 z=r1/r12;
>> 365 } while( ( (r12 > 1.0) || !((zmin <= z)&&(z <= zmax))) && ++loopCounter < maxNumberOfLoops );
>> 366 */
>> 367 if ( loopCounter >= maxNumberOfLoops ) {
>> 368 z = 0.5*(zmin + zmax); // Just a value between zmin and zmax, no physics considerations at all!
>> 369 }
442 370
443 created = QuarkPair.second; << 371 return z;
444 372
445 DecayQuark = decay->GetPDGEncoding(); << 373 }
446 NewQuark = created->GetPDGEncoding(); <<
447 374
448 G4ParticleDefinition * had=hadronizer->B << 375 return z;
449 return had; <<
450 } <<
451 } 376 }
452 <<
453 //-------------------------------------------- 377 //-----------------------------------------------------------------------------------------
454 378
455 G4LorentzVector * G4QGSMFragmentation::SplitEa 379 G4LorentzVector * G4QGSMFragmentation::SplitEandP(G4ParticleDefinition * pHadron,
456 380 G4FragmentingString * string,
457 381 G4FragmentingString * NewString)
458 { 382 {
459 G4double HadronMass = pHadron->GetPDGMa 383 G4double HadronMass = pHadron->GetPDGMass();
460 384
461 SetMinimalStringMass(NewString); 385 SetMinimalStringMass(NewString);
462 386
463 if ( MinimalStringMass < 0.0 ) return n <<
464 <<
465 #ifdef debug_QGSMfragmentation 387 #ifdef debug_QGSMfragmentation
466 G4cout<<"G4QGSMFragmentation::SplitEand 388 G4cout<<"G4QGSMFragmentation::SplitEandP "<<pHadron->GetParticleName()<<G4endl;
467 G4cout<<"String 4 mom, String M "<<stri 389 G4cout<<"String 4 mom, String M "<<string->Get4Momentum()<<" "<<string->Mass()<<G4endl;
468 G4cout<<"HadM MinimalStringMassLeft Str 390 G4cout<<"HadM MinimalStringMassLeft StringM hM+sM "<<HadronMass<<" "<<MinimalStringMass<<" "
469 <<string->Mass()<<" "<<HadronMass 391 <<string->Mass()<<" "<<HadronMass+MinimalStringMass<<G4endl;
470 #endif 392 #endif
471 393
472 if (HadronMass + MinimalStringMass > st 394 if (HadronMass + MinimalStringMass > string->Mass())
473 { 395 {
474 #ifdef debug_QGSMfragmentation 396 #ifdef debug_QGSMfragmentation
475 G4cout<<"Mass of the string is not su 397 G4cout<<"Mass of the string is not sufficient to produce the hadron!"<<G4endl;
476 #endif 398 #endif
477 return 0; 399 return 0;
478 } // have to start all over! 400 } // have to start all over!
479 401
480 // calculate and assign hadron transver 402 // calculate and assign hadron transverse momentum component HadronPx andHadronPy
481 G4double StringMT2 = string->MassT2(); 403 G4double StringMT2 = string->MassT2();
482 G4double StringMT = std::sqrt(StringMT 404 G4double StringMT = std::sqrt(StringMT2);
483 405
484 G4LorentzVector String4Momentum = strin 406 G4LorentzVector String4Momentum = string->Get4Momentum();
485 String4Momentum.setPz(0.); 407 String4Momentum.setPz(0.);
486 G4ThreeVector StringPt = String4Momentu 408 G4ThreeVector StringPt = String4Momentum.vect();
487 409
488 G4ThreeVector HadronPt , RemSysPt; 410 G4ThreeVector HadronPt , RemSysPt;
489 G4double HadronMassT2, ResidualMas 411 G4double HadronMassT2, ResidualMassT2;
490 412
491 // Mt distribution is implemented <<
492 G4double HadronMt, Pt, Pt2, phi; <<
493 <<
494 //... sample Pt of the hadron 413 //... sample Pt of the hadron
495 G4int attempt=0; 414 G4int attempt=0;
496 do 415 do
497 { 416 {
498 attempt++; if (attempt > StringLoopIn 417 attempt++; if (attempt > StringLoopInterrupt) return 0;
499 418
500 HadronMt = HadronMass - 200.0*G4Log(G << 419 HadronPt =SampleQuarkPt() + string->DecayPt();
501 Pt2 = sqr(HadronMt)-sqr(HadronMass); <<
502 phi = 2.*pi*G4UniformRand(); <<
503 G4ThreeVector SampleQuarkPtw= G4Three <<
504 HadronPt =SampleQuarkPtw + string->D <<
505 HadronPt.setZ(0); 420 HadronPt.setZ(0);
506 RemSysPt = StringPt - HadronPt; 421 RemSysPt = StringPt - HadronPt;
507 422
508 HadronMassT2 = sqr(HadronMass) + Hadr 423 HadronMassT2 = sqr(HadronMass) + HadronPt.mag2();
509 ResidualMassT2=sqr(MinimalStringMass) 424 ResidualMassT2=sqr(MinimalStringMass) + RemSysPt.mag2();
510 425
511 } while (std::sqrt(HadronMassT2) + std: 426 } while (std::sqrt(HadronMassT2) + std::sqrt(ResidualMassT2) > StringMT); /* Loop checking, 07.08.2015, A.Ribon */
512 427
513 //... sample z to define hadron longit 428 //... sample z to define hadron longitudinal momentum and energy
514 //... but first check the available pha 429 //... but first check the available phase space
515 430
516 G4double Pz2 = (sqr(StringMT2 - HadronM 431 G4double Pz2 = (sqr(StringMT2 - HadronMassT2 - ResidualMassT2) -
517 4*HadronMassT2 * ResidualMassT2)/4./ 432 4*HadronMassT2 * ResidualMassT2)/4./StringMT2;
518 433
519 if ( Pz2 < 0 ) {return 0;} // 434 if ( Pz2 < 0 ) {return 0;} // have to start all over!
520 435
521 //... then compute allowed z region z_ 436 //... then compute allowed z region z_min <= z <= z_max
522 437
523 G4double Pz = std::sqrt(Pz2); 438 G4double Pz = std::sqrt(Pz2);
524 G4double zMin = (std::sqrt(HadronMassT2 439 G4double zMin = (std::sqrt(HadronMassT2+Pz2) - Pz)/std::sqrt(StringMT2);
525 G4double zMax = (std::sqrt(HadronMassT2 440 G4double zMax = (std::sqrt(HadronMassT2+Pz2) + Pz)/std::sqrt(StringMT2);
526 441
527 if (zMin >= zMax) return 0; // have 442 if (zMin >= zMax) return 0; // have to start all over!
528 443
529 G4double z = GetLightConeZ(zMin, zMax, 444 G4double z = GetLightConeZ(zMin, zMax,
530 string->GetDecayParton() 445 string->GetDecayParton()->GetPDGEncoding(), pHadron,
531 HadronPt.x(), HadronPt.y 446 HadronPt.x(), HadronPt.y());
532 447
533 //... now compute hadron longitudinal m 448 //... now compute hadron longitudinal momentum and energy
534 // longitudinal hadron momentum compone 449 // longitudinal hadron momentum component HadronPz
535 450
536 HadronPt.setZ( 0.5* string->GetDecayDir 451 HadronPt.setZ( 0.5* string->GetDecayDirection() *
537 (z * string->LightConeDecay() - 452 (z * string->LightConeDecay() -
538 HadronMassT2/(z * string->LightCone 453 HadronMassT2/(z * string->LightConeDecay())) );
539 G4double HadronE = 0.5* (z * string->L 454 G4double HadronE = 0.5* (z * string->LightConeDecay() +
540 HadronMassT2/(z * string->LightConeDe 455 HadronMassT2/(z * string->LightConeDecay()) );
541 456
542 G4LorentzVector * a4Momentum= new G4Lor 457 G4LorentzVector * a4Momentum= new G4LorentzVector(HadronPt,HadronE);
543 458
544 #ifdef debug_QGSMfragmentation 459 #ifdef debug_QGSMfragmentation
545 G4cout<<"string->GetDecayDirection() st 460 G4cout<<"string->GetDecayDirection() string->LightConeDecay() "
546 <<string->GetDecayDirection()<<" 461 <<string->GetDecayDirection()<<" "<<string->LightConeDecay()<<G4endl;
547 G4cout<<"HadronPt,HadronE "<<HadronPt<< 462 G4cout<<"HadronPt,HadronE "<<HadronPt<<" "<<HadronE<<G4endl;
548 G4cout<<"Out of QGSM SplitEandP "<<G4en 463 G4cout<<"Out of QGSM SplitEandP "<<G4endl;
549 #endif 464 #endif
550 465
551 return a4Momentum; 466 return a4Momentum;
552 } 467 }
553 468
554 //-------------------------------------------- <<
555 <<
556 G4double G4QGSMFragmentation::GetLightConeZ(G4 <<
557 G4 <<
558 { <<
559 G4double lambda = 2.0*(sqr(ptx)+sqr(pty))/sq <<
560 <<
561 #ifdef debug_QGSMfragmentation <<
562 G4cout<<"GetLightConeZ zmin zmax Parton pHad <<
563 <<" "/*<< pHadron->GetParticleName() * <<
564 #endif <<
565 <<
566 G4double z(0.); <<
567 G4int DiQold(0), DiQnew(0); <<
568 G4double d1(-1.0), d2(0.); <<
569 G4double invD1(0.),invD2(0.), r1(0.),r2(0.), <<
570 <<
571 G4int absDecayQuarkCode = std::abs( DecayQua <<
572 G4int absNewQuarkCode = std::abs( NewQuark <<
573 <<
574 G4int q1(0), q2(0); <<
575 // q1 = absDecayQuarkCode/1000; q2 = (absDe <<
576 <<
577 G4int qA(0), qB(0); <<
578 // qA = absNewQuarkCode/1000; qB = (absNe <<
579 <<
580 if ( (absDecayQuarkCode < 6) && (absNewQuark <<
581 d1 = FFq2q[absDecayQuarkCode-1][absNewQuar <<
582 } <<
583 <<
584 if ( (absDecayQuarkCode < 6) && (absNewQuark <<
585 qA = absNewQuarkCode/1000; qB = (absNewQu <<
586 d1 = FFq2qq[absDecayQuarkCode-1][DiQnew][0] <<
587 } <<
588 <<
589 if ( (absDecayQuarkCode > 6) && (absNewQuark <<
590 q1 = absDecayQuarkCode/1000; q2 = (absDeca <<
591 d1 = FFqq2q[DiQold][absNewQuarkCode-1][0]; <<
592 } <<
593 <<
594 if ( d1 < 0. ) { <<
595 q1 = absDecayQuarkCode/1000; q2 = (absDeca <<
596 qA = absNewQuarkCode/1000; qB = (absNewQ <<
597 d1 = FFqq2qq[DiQold][DiQnew][0]; d2 = FFqq <<
598 } <<
599 <<
600 d2 +=lambda; <<
601 d1+=1.0; d2+=1.0; <<
602 <<
603 invD1=1./d1; invD2=1./d2; <<
604 <<
605 const G4int maxNumberOfLoops = 10000; <<
606 G4int loopCounter = 0; <<
607 do // Jong's algorithm <<
608 { <<
609 r1=G4Pow::GetInstance()->powA(G4UniformRan <<
610 r2=G4Pow::GetInstance()->powA(G4UniformRan <<
611 r12=r1+r2; <<
612 z=r1/r12; <<
613 } while ( ( (r12 > 1.0) || !((zmin <= z)&&(z <<
614 ++loopCounter < maxNumberOfLoops ) <<
615 <<
616 if ( loopCounter >= maxNumberOfLoops ) { <<
617 z = 0.5*(zmin + zmax); // Just a value be <<
618 } <<
619 <<
620 return z; <<
621 } <<
622 469
623 //-------------------------------------------- 470 //-----------------------------------------------------------------------------------------
624 471
625 G4bool G4QGSMFragmentation::SplitLast(G4Fragme 472 G4bool G4QGSMFragmentation::SplitLast(G4FragmentingString * string,
626 G4KineticTrackVector * Lef 473 G4KineticTrackVector * LeftVector,
627 G4KineticTrackVector * Right 474 G4KineticTrackVector * RightVector)
628 { 475 {
629 //... perform last cluster decay 476 //... perform last cluster decay
630 477
631 G4ThreeVector ClusterVel =string->Get4Mome 478 G4ThreeVector ClusterVel =string->Get4Momentum().boostVector();
632 G4double ResidualMass =string->Mass(); 479 G4double ResidualMass =string->Mass();
633 480
634 #ifdef debug_QGSMfragmentation 481 #ifdef debug_QGSMfragmentation
635 G4cout<<"Split last----------------------- 482 G4cout<<"Split last-----------------------------------------"<<G4endl;
636 G4cout<<"StrMass "<<ResidualMass<<" q's " 483 G4cout<<"StrMass "<<ResidualMass<<" q's "
637 <<string->GetLeftParton()->GetPartic 484 <<string->GetLeftParton()->GetParticleName()<<" "
638 <<string->GetRightParton()->GetParti 485 <<string->GetRightParton()->GetParticleName()<<G4endl;
639 #endif 486 #endif
640 487
641 G4int cClusterInterrupt = 0; 488 G4int cClusterInterrupt = 0;
642 G4ParticleDefinition *LeftHadron = nullptr << 489 G4ParticleDefinition * LeftHadron, * RightHadron;
643 G4ParticleDefinition *RightHadron = nullpt <<
644 const G4int maxNumberOfLoops = 1000; 490 const G4int maxNumberOfLoops = 1000;
645 G4int loopCounter = 0; 491 G4int loopCounter = 0;
646 <<
647 G4double LeftHadronMass(0.); G4double Righ <<
648 do 492 do
649 { 493 {
650 if (cClusterInterrupt++ >= ClusterLoop << 494 if (cClusterInterrupt++ >= ClusterLoopInterrupt)
651 LeftHadronMass = -MaxMass; RightHadron << 495 {
>> 496 return false;
>> 497 }
652 498
653 G4ParticleDefinition * quark = nullptr; << 499 G4ParticleDefinition * quark = NULL;
654 string->SetLeftPartonStable(); // to query q 500 string->SetLeftPartonStable(); // to query quark contents..
655 501
656 if (string->DecayIsQuark() && string->Stable 502 if (string->DecayIsQuark() && string->StableIsQuark() )
657 { 503 {
658 //... there are quarks on cluster ends 504 //... there are quarks on cluster ends
659 505
660 G4int IsParticle=(string->GetLeftParton()- << 506 G4int IsParticle=(string->GetLeftParton()->GetPDGEncoding()>0) ? -1 : +1; // if we have a quark,
661 // if we have a quark, we need << 507 // we need antiquark or diquark
662 <<
663 pDefPair QuarkPair = CreatePartonPair(IsPa 508 pDefPair QuarkPair = CreatePartonPair(IsParticle);
664 quark = QuarkPair.second; 509 quark = QuarkPair.second;
665 510
666 LeftHadron= hadronizer->Build(QuarkPair.fi << 511 LeftHadron= hadronizer->BuildLowSpin(QuarkPair.first, string->GetLeftParton());
667 if ( LeftHadron == NULL ) continue; << 512 } else {
668 RightHadron = hadronizer->Build(stri << 513 //... there is a Diquark on cluster ends
669 if ( RightHadron == NULL ) continue; <<
670 } else if( (!string->DecayIsQuark() && stri <<
671 ( string->DecayIsQuark() && !stri <<
672 //... there is a Diquark on one of cluster <<
673 G4int IsParticle; 514 G4int IsParticle;
674 if ( string->StableIsQuark() ) { 515 if ( string->StableIsQuark() ) {
675 IsParticle=(string->GetLeftParton()->Get 516 IsParticle=(string->GetLeftParton()->GetPDGEncoding()>0) ? -1 : +1;
676 } else { 517 } else {
677 IsParticle=(string->GetLeftParton()->Get 518 IsParticle=(string->GetLeftParton()->GetPDGEncoding()>0) ? +1 : -1;
678 } 519 }
679 520
680 //G4double ProbSaS = 1.0 - 2.0 * G 521 //G4double ProbSaS = 1.0 - 2.0 * GetStrangeSuppress();
681 //G4double ActualProb = ProbSaS * 1. 522 //G4double ActualProb = ProbSaS * 1.4;
682 //SetStrangenessSuppression((1.0-Act 523 //SetStrangenessSuppression((1.0-ActualProb)/2.0);
683 524
684 pDefPair QuarkPair = CreatePartonPai 525 pDefPair QuarkPair = CreatePartonPair(IsParticle,false); // no diquarks wanted
685 //SetStrangenessSuppression((1.0-Pro 526 //SetStrangenessSuppression((1.0-ProbSaS)/2.0);
686 quark = QuarkPair.second; 527 quark = QuarkPair.second;
687 LeftHadron=hadronizer->Build(QuarkPa << 528 LeftHadron=hadronizer->BuildLowSpin(QuarkPair.first, string->GetLeftParton());
688 if ( LeftHadron == NULL ) continue; <<
689 RightHadron = hadronizer->Build(stri <<
690 if ( RightHadron == NULL ) continue; <<
691 } else { // Diquark and anti-diquark <<
692 if (cClusterInterrupt++ >= ClusterLo <<
693 G4int LeftQuark1= string->GetLeftPar <<
694 G4int LeftQuark2=(string->GetLeftPar <<
695 G4int RightQuark1= string->GetRightP <<
696 G4int RightQuark2=(string->GetRightP <<
697 if (G4UniformRand()<0.5) { <<
698 LeftHadron =hadronizer->Build(Fin <<
699 RightHadron =hadronizer->Build(Fin <<
700 } else { <<
701 LeftHadron =hadronizer->Build(Fin <<
702 RightHadron =hadronizer->Build(Fin <<
703 } <<
704 if ( (LeftHadron == NULL) || (RightHadron <<
705 } 529 }
706 LeftHadronMass = LeftHadron->GetPDGMa << 530
707 RightHadronMass = RightHadron->GetPDGM << 531 RightHadron = hadronizer->BuildLowSpin(string->GetRightParton(), quark);
708 //... repeat procedure, if mass of clu << 532
709 } while ( ( ResidualMass <= LeftHadronMass << 533 } while ( ( ResidualMass <= LeftHadron->GetPDGMass() + RightHadron->GetPDGMass() )
710 && ++loopCounter < maxNumberOfLo 534 && ++loopCounter < maxNumberOfLoops ); /* Loop checking, 07.08.2015, A.Ribon */
711 535
712 if ( loopCounter >= maxNumberOfLoops ) { 536 if ( loopCounter >= maxNumberOfLoops ) {
713 return false; 537 return false;
714 } 538 }
715 539
716 //... compute hadron momenta and energies 540 //... compute hadron momenta and energies
717 G4LorentzVector LeftMom, RightMom; 541 G4LorentzVector LeftMom, RightMom;
718 G4ThreeVector Pos; 542 G4ThreeVector Pos;
719 Sample4Momentum(&LeftMom , LeftHadron->Get 543 Sample4Momentum(&LeftMom , LeftHadron->GetPDGMass() ,
720 &RightMom, RightHadron->Ge 544 &RightMom, RightHadron->GetPDGMass(), ResidualMass);
721 LeftMom.boost(ClusterVel); 545 LeftMom.boost(ClusterVel);
722 RightMom.boost(ClusterVel); 546 RightMom.boost(ClusterVel);
723 547
724 #ifdef debug_QGSMfragmentation 548 #ifdef debug_QGSMfragmentation
725 G4cout<<LeftHadron->GetParticleName()<<" " 549 G4cout<<LeftHadron->GetParticleName()<<" "<<RightHadron->GetParticleName()<<G4endl;
726 G4cout<<"Left Hadrom P M "<<LeftMom<<" "< 550 G4cout<<"Left Hadrom P M "<<LeftMom<<" "<<LeftMom.mag()<<G4endl;
727 G4cout<<"Right Hadrom P M "<<RightMom<<" " 551 G4cout<<"Right Hadrom P M "<<RightMom<<" "<<RightMom.mag()<<G4endl;
728 #endif 552 #endif
729 553
730 LeftVector->push_back(new G4KineticTrack(L 554 LeftVector->push_back(new G4KineticTrack(LeftHadron, 0, Pos, LeftMom));
731 RightVector->push_back(new G4KineticTrack( 555 RightVector->push_back(new G4KineticTrack(RightHadron, 0, Pos, RightMom));
732 556
733 return true; 557 return true;
>> 558
>> 559 }
>> 560
>> 561 //----------------------------------------------------------------------------------------------------------
>> 562
>> 563 G4bool G4QGSMFragmentation::IsFragmentable(const G4FragmentingString * const string)
>> 564 {
>> 565 return sqr( FragmentationMass(string) + MassCut ) < string->Mass2();
>> 566 }
>> 567
>> 568 //----------------------------------------------------------------------------------------------------------
>> 569
>> 570 G4bool G4QGSMFragmentation::StopFragmenting(const G4FragmentingString * const string)
>> 571 {
>> 572 SetMinimalStringMass(string);
>> 573 if (string->FourQuarkString())
>> 574 {
>> 575 return G4UniformRand() < G4Exp(-0.0005*(string->Mass() - MinimalStringMass));
>> 576 } else {
>> 577 G4bool Result = G4UniformRand() <
>> 578 G4Exp(-0.66e-6*(string->Mass()*string->Mass() - MinimalStringMass*MinimalStringMass));
>> 579 // G4bool Result = string->Mass() < MinimalStringMass + 150.*MeV*G4UniformRand(); // a'la LUND
>> 580
>> 581 #ifdef debug_QGSMfragmentation
>> 582 G4cout<<"StopFragmenting MinimalStringMass string->Mass() "<<MinimalStringMass<<" "<<string->Mass()<<G4endl;
>> 583 G4cout<<"StopFragmenting - Yes/No "<<Result<<G4endl;
>> 584 #endif
>> 585 return Result;
>> 586 }
734 } 587 }
735 588
>> 589
736 //-------------------------------------------- 590 //----------------------------------------------------------------------------------------------------------
737 591
738 void G4QGSMFragmentation::Sample4Momentum(G4Lo 592 void G4QGSMFragmentation::Sample4Momentum(G4LorentzVector* Mom , G4double Mass ,
739 G4Lo 593 G4LorentzVector* AntiMom, G4double AntiMass, G4double InitialMass)
740 { 594 {
741 #ifdef debug_QGSMfragmentation <<
742 G4cout<<"Sample4Momentum Last------------- <<
743 G4cout<<" StrMass "<<InitialMass<<" Mass1 <<
744 G4cout<<" SumMass "<<Mass+AntiMass<<G4end <<
745 #endif <<
746 <<
747 G4double r_val = sqr(InitialMass*InitialMa 595 G4double r_val = sqr(InitialMass*InitialMass - Mass*Mass - AntiMass*AntiMass) - sqr(2.*Mass*AntiMass);
748 G4double Pabs = (r_val > 0.)? std::sqrt(r_ 596 G4double Pabs = (r_val > 0.)? std::sqrt(r_val)/(2.*InitialMass) : 0;
749 597
750 //... sample unit vector 598 //... sample unit vector
751 G4double pz = 1. - 2.*G4UniformRand(); 599 G4double pz = 1. - 2.*G4UniformRand();
752 G4double st = std::sqrt(1. - pz * pz)* 600 G4double st = std::sqrt(1. - pz * pz)*Pabs;
753 G4double phi = 2.*pi*G4UniformRand(); 601 G4double phi = 2.*pi*G4UniformRand();
754 G4double px = st*std::cos(phi); 602 G4double px = st*std::cos(phi);
755 G4double py = st*std::sin(phi); 603 G4double py = st*std::sin(phi);
756 pz *= Pabs; 604 pz *= Pabs;
757 605
758 Mom->setPx(px); Mom->setPy(py); Mom->setPz 606 Mom->setPx(px); Mom->setPy(py); Mom->setPz(pz);
759 Mom->setE(std::sqrt(Pabs*Pabs + Mass*Mass) 607 Mom->setE(std::sqrt(Pabs*Pabs + Mass*Mass));
760 608
761 AntiMom->setPx(-px); AntiMom->setPy(-py); 609 AntiMom->setPx(-px); AntiMom->setPy(-py); AntiMom->setPz(-pz);
762 AntiMom->setE (std::sqrt(Pabs*Pabs + AntiM 610 AntiMom->setE (std::sqrt(Pabs*Pabs + AntiMass*AntiMass));
763 } 611 }
764 612
765 //-------------------------------------------- << 613 //-----------------------------------------------------------------------------
766 614
767 void G4QGSMFragmentation::SetFFq2q() // q-> q << 615 G4ParticleDefinition *G4QGSMFragmentation::DiQuarkSplitup( G4ParticleDefinition* decay,
>> 616 G4ParticleDefinition *&created )
768 { 617 {
769 for (G4int i=0; i < 5; i++) { << 618 //... can Diquark break or not?
770 FFq2q[i][0][0] = 2.0 ; FFq2q[i][0][1] = -a << 619 if (G4UniformRand() < DiquarkBreakProb ) //... Diquark break
771 FFq2q[i][1][0] = 2.0 ; FFq2q[i][1][1] = -a << 620 {
772 FFq2q[i][2][0] = 1.0 ; FFq2q[i][2][1] = -a << 621 G4int stableQuarkEncoding = decay->GetPDGEncoding()/1000;
773 FFq2q[i][3][0] = 1.0 ; FFq2q[i][3][1] = -a << 622 G4int decayQuarkEncoding = (decay->GetPDGEncoding()/100)%10;
774 FFq2q[i][4][0] = 1.0 ; FFq2q[i][4][1] = -a << 623
775 } << 624 if (G4UniformRand() < 0.5)
>> 625 {
>> 626 G4int Swap = stableQuarkEncoding;
>> 627 stableQuarkEncoding = decayQuarkEncoding;
>> 628 decayQuarkEncoding = Swap;
>> 629 }
>> 630
>> 631 G4int IsParticle=(decayQuarkEncoding>0) ? -1 : +1; // if we have a quark, we need antiquark
>> 632
>> 633 G4double StrSup=GetStrangeSuppress();
>> 634 SetStrangenessSuppression((1.0 - 0.07)/2.);
>> 635 pDefPair QuarkPair = CreatePartonPair(IsParticle,false); // no diquarks wanted
>> 636 SetStrangenessSuppression(StrSup);
>> 637
>> 638 //... Build new Diquark
>> 639 G4int QuarkEncoding=QuarkPair.second->GetPDGEncoding();
>> 640 G4int i10 = std::max(std::abs(QuarkEncoding), std::abs(stableQuarkEncoding));
>> 641 G4int i20 = std::min(std::abs(QuarkEncoding), std::abs(stableQuarkEncoding));
>> 642 G4int spin = (i10 != i20 && G4UniformRand() <= 0.5)? 1 : 3;
>> 643 G4int NewDecayEncoding = -1*IsParticle*(i10 * 1000 + i20 * 100 + spin);
>> 644 created = FindParticle(NewDecayEncoding);
>> 645 G4ParticleDefinition * decayQuark=FindParticle(decayQuarkEncoding);
>> 646 G4ParticleDefinition * had=hadronizer->Build(QuarkPair.first, decayQuark);
>> 647
>> 648 return had;
>> 649
>> 650 } else { //... Diquark does not break
>> 651
>> 652 G4int IsParticle=(decay->GetPDGEncoding()>0) ? +1 : -1; // if we have a diquark, we need quark)
>> 653 G4double StrSup=GetStrangeSuppress(); // for changing s-sbar production
>> 654 SetStrangenessSuppression((1.0 - 0.07)/2.);
>> 655 pDefPair QuarkPair = CreatePartonPair(IsParticle,false); // no diquarks wanted
>> 656 SetStrangenessSuppression(StrSup);
>> 657
>> 658 created = QuarkPair.second;
>> 659
>> 660 G4ParticleDefinition * had=hadronizer->Build(QuarkPair.first, decay);
>> 661 return had;
>> 662 }
776 } 663 }
777 664
778 //-------------------------------------------- << 665 //-----------------------------------------------------------------------------
779 666
780 void G4QGSMFragmentation::SetFFq2qq() // q-> << 667 G4KineticTrack * G4QGSMFragmentation::Splitup( G4FragmentingString *string,
>> 668 G4FragmentingString *&newString)
781 { 669 {
782 for (G4int i=0; i < 5; i++) { << 670 #ifdef debug_QGSMfragmentation
783 FFq2qq[i][ 0][0] = 0.0 ; FFq2qq[i][ 0][1] << 671 G4cout<<G4endl;
784 FFq2qq[i][ 1][0] = 0.0 ; FFq2qq[i][ 1][1] << 672 G4cout<<"Start SplitUP (G4VLongitudinalStringDecay) ========================="<<G4endl;
785 FFq2qq[i][ 2][0] = 0.0 ; FFq2qq[i][ 2][1] << 673 G4cout<<"String partons: " <<string->GetLeftParton()->GetPDGEncoding()<<" "
786 FFq2qq[i][ 3][0] = 0.0 ; FFq2qq[i][ 3][1] << 674 <<string->GetRightParton()->GetPDGEncoding()<<" "
787 FFq2qq[i][ 4][0] = 0.0 ; FFq2qq[i][ 4][1] << 675 <<"Direction " <<string->GetDecayDirection()<<G4endl;
788 FFq2qq[i][ 5][0] = 0.0 ; FFq2qq[i][ 5][1] << 676 #endif
789 FFq2qq[i][ 6][0] = 0.0 ; FFq2qq[i][ 6][1] <<
790 FFq2qq[i][ 7][0] = 0.0 ; FFq2qq[i][ 7][1] <<
791 FFq2qq[i][ 8][0] = 0.0 ; FFq2qq[i][ 8][1] <<
792 FFq2qq[i][ 9][0] = 0.0 ; FFq2qq[i][ 9][1] <<
793 FFq2qq[i][10][0] = 0.0 ; FFq2qq[i][10][1] <<
794 FFq2qq[i][11][0] = 0.0 ; FFq2qq[i][11][1] <<
795 FFq2qq[i][12][0] = 0.0 ; FFq2qq[i][12][1] <<
796 FFq2qq[i][13][0] = 0.0 ; FFq2qq[i][13][1] <<
797 FFq2qq[i][14][0] = 0.0 ; FFq2qq[i][14][1] <<
798 } <<
799 } <<
800 677
801 //-------------------------------------------- << 678 //... random choice of string end to use for creating the hadron (decay)
>> 679 G4int SideOfDecay = (G4UniformRand() < 0.5)? 1: -1;
>> 680 if (SideOfDecay < 0)
>> 681 {
>> 682 string->SetLeftPartonStable();
>> 683 } else
>> 684 {
>> 685 string->SetRightPartonStable();
>> 686 }
>> 687
>> 688 G4ParticleDefinition *newStringEnd;
>> 689 G4ParticleDefinition * HadronDefinition;
>> 690 if (string->DecayIsQuark())
>> 691 {
>> 692 G4double ProbDqADq = GetDiquarkSuppress();
>> 693
>> 694 G4int NumberOfpossibleBaryons = 2;
>> 695
>> 696 if (string->GetLeftParton()->GetParticleSubType() != "quark") NumberOfpossibleBaryons++;
>> 697 if (string->GetRightParton()->GetParticleSubType() != "quark") NumberOfpossibleBaryons++;
>> 698
>> 699 G4double ActualProb = ProbDqADq ;
>> 700 ActualProb *= (1.0-G4Exp(2.0*(1.0 - string->Mass()/(NumberOfpossibleBaryons*1400.0))));
>> 701
>> 702 SetDiquarkSuppression(ActualProb);
>> 703
>> 704 HadronDefinition= QuarkSplitup(string->GetDecayParton(), newStringEnd);
>> 705
>> 706 SetDiquarkSuppression(ProbDqADq);
>> 707 } else {
>> 708 HadronDefinition= DiQuarkSplitup(string->GetDecayParton(), newStringEnd);
>> 709 }
802 710
803 void G4QGSMFragmentation::SetFFqq2q() // q1q2 << 711 #ifdef debug_QGSMfragmentation
>> 712 G4cout<<"The parton "<<string->GetDecayParton()->GetPDGEncoding()<<" "
>> 713 <<" produces hadron "<<HadronDefinition->GetParticleName()
>> 714 <<" and is transformed to "<<newStringEnd->GetPDGEncoding()<<G4endl;
>> 715 G4cout<<"The side of the string decay Left/Right (1/-1) "<<SideOfDecay<<G4endl;
>> 716 #endif
>> 717 // create new String from old, ie. keep Left and Right order, but replace decay
>> 718
>> 719 newString=new G4FragmentingString(*string,newStringEnd); // To store possible
>> 720 // quark containt of new string
>> 721
>> 722 #ifdef debug_QGSMfragmentation
>> 723 G4cout<<"An attempt to determine its energy (SplitEandP)"<<G4endl;
>> 724 #endif
>> 725 G4LorentzVector* HadronMomentum=SplitEandP(HadronDefinition, string, newString);
>> 726
>> 727 delete newString; newString=0;
>> 728
>> 729 G4KineticTrack * Hadron =0;
>> 730 if ( HadronMomentum != 0 ) {
>> 731
>> 732 #ifdef debug_QGSMfragmentation
>> 733 G4cout<<"The attempt was successful"<<G4endl;
>> 734 #endif
>> 735 G4ThreeVector Pos;
>> 736 Hadron = new G4KineticTrack(HadronDefinition, 0,Pos, *HadronMomentum);
>> 737
>> 738 newString=new G4FragmentingString(*string,newStringEnd,HadronMomentum);
>> 739
>> 740 delete HadronMomentum;
>> 741 }
>> 742 else
>> 743 {
>> 744 #ifdef debug_QGSMfragmentation
>> 745 G4cout<<"The attempt was not successful !!!"<<G4endl;
>> 746 #endif
>> 747 }
>> 748
>> 749 #ifdef debug_VStringDecay
>> 750 G4cout<<"End SplitUP (G4VLongitudinalStringDecay) ====================="<<G4endl;
>> 751 #endif
>> 752
>> 753 return Hadron;
>> 754 }
>> 755
>> 756 //---------------------------------------------------------------
>> 757 void G4QGSMFragmentation::SetMinMasses()
804 { 758 {
805 for (G4int i=0; i < 15; i++) { << 759 // ------ For estimation of a minimal string mass ---------------
806 FFqq2q[i][0][0] = 2.0*(arho - an); FFqq2q[ << 760 Mass_of_light_quark =140.*MeV;
807 FFqq2q[i][1][0] = 2.0*(arho - an); FFqq2q[ << 761 Mass_of_heavy_quark =500.*MeV;
808 FFqq2q[i][2][0] = 2.0*(arho - an); FFqq2q[ << 762 Mass_of_string_junction=720.*MeV;
809 FFqq2q[i][3][0] = 2.0*(arho - an); FFqq2q[ << 763
810 FFqq2q[i][4][0] = 2.0*(arho - an); FFqq2q[ << 764 G4double minMQQbarStr[3][3] ={ {350.0, 350.0, 710.0}, //DDbar, DUbar, DSbar in MeV
811 } << 765 {350.0, 350.0, 710.0}, //UDbar, UUbar, USbar in Mev
>> 766 {710.0, 710.0,1070.0 }};//SDbar, SUbar, SSbar in MeV
>> 767 for (G4int i=0; i<3; i++){ for (G4int j=0; j<3; j++){minMassQQbarStr[i][j]=minMQQbarStr[i][j];};};
>> 768
>> 769
>> 770 G4double minMQDiQStr[3][3][3] = {{{1160., 1160., 1340.}, {1160., 1160., 1340.}, {1340., 1340., 1540.},}, //d-dd, d-du, d-ds, d-ud, d-uu, d-us, d-sd, d-su, d-ss
>> 771 {{1160., 1160., 1340.}, {1160., 1160., 1340.}, {1340., 1340., 1540.},}, //u-dd, u-du, u-ds, u-ud, u-uu, u-us, u-sd, u-su, u-ss
>> 772 {{1520., 1520., 1690.}, {1520., 1520., 1690.}, {1690., 1690., 1890. }}};//s-dd, s-du, s-ds, s-ud, s-uu, s-us, s-sd, s-su, s-ss
>> 773 for (G4int i=0; i<3; i++){ for (G4int j=0; j<3; j++){ for (G4int k=0; k<3; k++){minMassQDiQStr[i][j][k]=minMQDiQStr[i][j][k];};};};
>> 774
>> 775 // ------ An estimated minimal string mass ----------------------
>> 776 MinimalStringMass = 0.;
>> 777 MinimalStringMass2 = 0.;
812 } 778 }
813 779
814 //-------------------------------------------- << 780 //--------------------------------------------------------------------------------------
>> 781 void G4QGSMFragmentation::SetMinimalStringMass(const G4FragmentingString * const string)
>> 782 {
>> 783 G4double EstimatedMass=0.;
>> 784
>> 785 G4int Qleft =std::abs(string->GetLeftParton()->GetPDGEncoding());
>> 786 G4int Qright=std::abs(string->GetRightParton()->GetPDGEncoding());
815 787
816 void G4QGSMFragmentation::SetFFqq2qq() // q1( << 788 if ((Qleft < 4) && (Qright < 4)) { // Q-Qbar string
>> 789 EstimatedMass=minMassQQbarStr[Qleft-1][Qright-1];
>> 790 MinimalStringMass=EstimatedMass;
>> 791 SetMinimalStringMass2(EstimatedMass);
>> 792 return;
>> 793 }
>> 794
>> 795 if ((Qleft < 4) && (Qright > 1000)) { // Q - DiQ string
>> 796 G4int q1=Qright/1000;
>> 797 G4int q2=(Qright/100)%10;
>> 798 EstimatedMass=minMassQDiQStr[Qleft-1][q1-1][q2-1];
>> 799 MinimalStringMass=EstimatedMass;
>> 800 SetMinimalStringMass2(EstimatedMass);
>> 801 return;
>> 802 }
>> 803
>> 804 if ((Qleft > 1000) && (Qright < 4)) { // DiQ - Q string
>> 805 G4int q1=Qleft/1000;
>> 806 G4int q2=(Qleft/100)%10;
>> 807 EstimatedMass=minMassQDiQStr[Qright-1][q1-1][q2-1];
>> 808 MinimalStringMass=EstimatedMass;
>> 809 SetMinimalStringMass2(EstimatedMass);
>> 810 return;
>> 811 }
>> 812
>> 813 // DiQuark - Anti DiQuark string -----------------
>> 814 G4int Number_of_quarks=0;
>> 815 G4int Number_of_squarks=0;
>> 816
>> 817 G4double StringM=string->Get4Momentum().mag();
>> 818
>> 819 #ifdef debug_QGSMfragmentation
>> 820 // G4cout<<"MinStringMass// Input String mass "<<string->Get4Momentum().mag()<<" Qleft "<<Qleft<<G4endl;
>> 821 #endif
>> 822
>> 823 if ( Qleft > 1000)
>> 824 {
>> 825 Number_of_quarks+=2;
>> 826 G4int q1=Qleft/1000;
>> 827 if ( q1 < 3) {EstimatedMass +=Mass_of_light_quark;}
>> 828 if ( q1 > 2) {EstimatedMass +=Mass_of_heavy_quark; Number_of_squarks++;}
>> 829
>> 830 G4int q2=(Qleft/100)%10;
>> 831 if ( q2 < 3) {EstimatedMass +=Mass_of_light_quark;}
>> 832 if ( q2 > 2) {EstimatedMass +=Mass_of_heavy_quark; Number_of_squarks++;}
>> 833 }
>> 834
>> 835 #ifdef debug_QGSMfragmentation
>> 836 // G4cout<<"Min mass with Qleft "<<Qleft<<" "<<EstimatedMass<<G4endl;
>> 837 #endif
>> 838
>> 839 if ( Qright > 1000)
>> 840 {
>> 841 Number_of_quarks+=2;
>> 842 G4int q1=Qright/1000;
>> 843 if ( q1 < 3) {EstimatedMass +=Mass_of_light_quark;}
>> 844 if ( q1 > 2) {EstimatedMass +=Mass_of_heavy_quark; Number_of_squarks++;}
>> 845
>> 846 G4int q2=(Qright/100)%10;
>> 847 if ( q2 < 3) {EstimatedMass +=Mass_of_light_quark;}
>> 848 if ( q2 > 2) {EstimatedMass +=Mass_of_heavy_quark; Number_of_squarks++;}
>> 849 //EstimatedMass +=Mass_of_string_junction;
>> 850 }
>> 851
>> 852 #ifdef debug_QGSMfragmentation
>> 853 // G4cout<<"Min mass with Qleft and Qright "<<Qright<<" "<<EstimatedMass<<G4endl;
>> 854 // G4cout<<"Number_of_quarks "<<Number_of_quarks<<" Number_of_squarks "<<Number_of_squarks<<G4endl;
>> 855 #endif
>> 856
>> 857 if (Number_of_quarks==4)
>> 858 {
>> 859 if (StringM > 1880.) { // 2*Mn = 1880
>> 860 if (Number_of_squarks==0) {EstimatedMass += 1320.*MeV;} //560+1320=1880=2*Mn
>> 861 else if (Number_of_squarks==1) {EstimatedMass += 1150.*MeV;} //920+1150=2070=M(Lam+N)
>> 862 else if (Number_of_squarks==2) {EstimatedMass += 960.*MeV;} //1280+960=2240= 2*M Lam
>> 863 else if (Number_of_squarks==3) {EstimatedMass += 800.*MeV;} //1640+800=2440=Mxi+Mlam
>> 864 else if (Number_of_squarks==4) {EstimatedMass += 640.*MeV;} //2000+640=2640=2*Mxi
>> 865 else {}
>> 866 }
>> 867 else
>> 868 {
>> 869 if (Number_of_squarks < 3) {EstimatedMass -= 200.*MeV;}
>> 870 else if (Number_of_squarks==3) {EstimatedMass -= 50.*MeV;}
>> 871 else if (Number_of_squarks==4) {EstimatedMass -= 40.*MeV;}
>> 872 else {}
>> 873 }
>> 874 }
>> 875
>> 876 #ifdef debug_QGSMfragmentation
>> 877 // G4cout<<"EstimatedMass -------------------- "<<EstimatedMass <<G4endl;
>> 878 #endif
>> 879
>> 880 MinimalStringMass=EstimatedMass;
>> 881 SetMinimalStringMass2(EstimatedMass);
>> 882 }
>> 883
>> 884 //--------------------------------------------------------------------------------------
>> 885 void G4QGSMFragmentation::SetMinimalStringMass2(const G4double aValue)
817 { 886 {
818 for (G4int i=0; i < 15; i++) { << 887 MinimalStringMass2=aValue * aValue;
819 FFqq2qq[i][0][0] = 0. ; FFqq2qq[i][0][1] <<
820 FFqq2qq[i][1][0] = 0. ; FFqq2qq[i][1][1] <<
821 FFqq2qq[i][2][0] = 0. ; FFqq2qq[i][2][1] <<
822 FFqq2qq[i][3][0] = 0. ; FFqq2qq[i][3][1] <<
823 FFqq2qq[i][4][0] = 0. ; FFqq2qq[i][4][1] <<
824 } <<
825 } 888 }
826 889
827 890