| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 // Historic fragment from M.Komogorov; clean-u
27
28 #include "G4ExcitedStringDecay.hh"
29 #include "G4SystemOfUnits.hh"
30 #include "G4KineticTrack.hh"
31 #include "G4LundStringFragmentation.hh"
32 #include "G4HadronicInteractionRegistry.hh"
33
34 #include "G4SampleResonance.hh"
35
36 //#define debug_G4ExcitedStringDecay
37 //#define debug_G4ExcitedStringCorr
38
39 G4ExcitedStringDecay::G4ExcitedStringDecay(G4V
40 : G4VStringFragmentation(), theStringDecay(p
41 {
42 if(!ptr) {
43 G4HadronicInteraction* p =
44 G4HadronicInteractionRegistry::Instance(
45 theStringDecay = static_cast<G4VLongitudin
46 if(!theStringDecay) { theStringDecay = new
47 }
48 SetModelName(theStringDecay->GetModelName())
49 }
50
51 G4ExcitedStringDecay::~G4ExcitedStringDecay()
52 {}
53
54 G4KineticTrackVector *G4ExcitedStringDecay::Fr
55 {
56 return theStringDecay->FragmentString(theStr
57 }
58
59 G4KineticTrackVector *G4ExcitedStringDecay::Fr
60 {
61 G4LorentzVector KTsum(0.,0.,0.,0.);
62
63 #ifdef debug_G4ExcitedStringDecay
64 G4cout<<G4endl;
65 G4cout<<"--------------------------- G4Excit
66 G4cout<<"Hadronization of Excited Strings: t
67 #endif
68
69 for ( unsigned int astring=0; astring < theS
70 // for ( unsigned int astring=0; astring < 1
71 {
72 // G4cout<<"theStrings->operator[](astring
73 if ( theStrings->operator[](astring)->IsEx
74 {KTsum+= theStrings->operator[](astri
75 else {KTsum+=theStrings->operator[](astrin
76 }
77
78 G4LorentzRotation toCms( -1 * KTsum.boostVec
79 G4LorentzRotation toLab(toCms.inverse());
80 G4LorentzVector Ptmp;
81 KTsum=G4LorentzVector(0.,0.,0.,0.);
82
83 for ( unsigned int astring=0; astring < theS
84 // for ( unsigned int astring=0; astring < 1
85 {
86 if ( theStrings->operator[](astring)->IsEx
87 {
88 Ptmp=toCms * theStrings->operator[](astr
89 theStrings->operator[](astring)->GetLeft
90
91 Ptmp=toCms * theStrings->operator[](astr
92 theStrings->operator[](astring)->GetRigh
93
94 KTsum+= theStrings->operator[](astring)-
95 }
96 else
97 {
98 Ptmp=toCms * theStrings->operator[](astr
99 theStrings->operator[](astring)->GetKine
100 KTsum+= theStrings->operator[](astring)-
101 }
102 }
103
104 G4SampleResonance BrW;
105 const G4ParticleDefinition* TrackDefinition=
106
107 G4KineticTrackVector * theResult = new G4Kin
108 G4int attempts(0);
109 G4bool success=false;
110 G4bool NeedEnergyCorrector=false;
111 do {
112 #ifdef debug_G4ExcitedStringDecay
113 G4cout<<"New try No "<<attempts<<" to
114 #endif
115
116 std::for_each(theResult->begin() , theResult
117 theResult->clear();
118
119 attempts++;
120
121 G4LorentzVector KTsecondaries(0.,0.,0.,0.);
122 NeedEnergyCorrector=false;
123
124 for ( unsigned int astring=0; astring < theS
125 // for ( unsigned int astring=0; astri
126 // for ( unsigned int astring=1; astri
127 // for ( unsigned int astring=0; astri
128 {
129 #ifdef debug_G4ExcitedStringDecay
130 G4cout<<"String No "<<astring+1<<" E
131
132 G4cout<<"String No "<<astring+1<<" 4
133 <<" "<<theStrings->operator[](astrin
134 #endif
135
136 G4KineticTrackVector * generatedKine
137 if ( theStrings->operator[](astring)->IsEx
138 {
139 #ifdef debug_G4ExcitedStringDecay
140 G4cout<<"Fragment String with par
141 <<theStrings->operator[](as
142 <<theStrings->operator[](as
143 <<"Direction "<<theStrings-
144 #endif
145 generatedKineticTracks=FragmentString
146 #ifdef debug_G4ExcitedStringDecay
147 G4cout<<"(G4ExcitedStringDecay) N
148 <<generatedKineticTracks->s
149 #endif
150 } else {
151 #ifdef debug_G4ExcitedStringDecay
152 G4cout<<" GetTrack from the Str
153 #endif
154 G4LorentzVector Mom=theStrings->o
155 G4KineticTrack * aTrack= new G4Ki
156 theStrings->operat
157 theStrings->operat
158 G4ThreeVector(0),
159
160 aTrack->SetPosition(theStrings->o
161
162 #ifdef debug_G4ExcitedStringDecay
163 G4cout<<" A particle stored in
164 #endif
165
166 generatedKineticTracks = new G4KineticT
167 generatedKineticTracks->push_back(aTrac
168 }
169
170 if (generatedKineticTracks == nullptr || g
171 {
172 // G4cerr << "G4VPartonStringMode
173 delete generatedKineticTracks;
174 generatedKineticTracks = nullptr;
175 continue;
176 }
177
178 G4LorentzVector KTsum1(0.,0.,0.,0.);
179 for ( unsigned int aTrack=0; aTrack<
180 {
181 #ifdef debug_G4ExcitedStringDecay
182 G4cout<<"Prod part No. "<<aTrack+
183 <<(*generatedKineticTracks)
184 <<(*generatedKineticTracks)
185 <<(*generatedKineticTracks)
186 #endif
187 // --------------- Sampling mass
188 TrackDefinition = (*generatedKine
189
190 if (TrackDefinition->IsShortLived())
191 {
192 G4double NewTrackMass =
193 BrW.SampleMass( TrackDefiniti
194 BrW.GetMinimu
195 TrackDefiniti
196 G4LorentzVector Tmp=G4LorentzVe
197 Tmp.setE(std::sqrt(sqr(NewTrack
198
199 (*generatedKineticTracks)[aTrac
200
201 #ifdef debug_G4ExcitedStringDec
202 G4cout<<"Resonance *** "<<aTrac
203 <<(*generatedKineticTrack
204 <<(*generatedKineticTrack
205 <<(*generatedKineticTrack
206 #endif
207 }
208 //-------------------------------
209
210 theResult->push_back(generatedKin
211 KTsum1+= (*generatedKineticTracks
212 }
213 KTsecondaries+=KTsum1;
214
215 #ifdef debug_G4ExcitedStringDecay
216 G4cout << "String secondaries(" <<ge
217 <<"Init string momentum: "<
218 <<"Final hadrons momentum: "<
219 #endif
220
221 if ( KTsum1.e() > 0 &&
222 std::abs((KTsum1.e()-theString
223 {
224 NeedEnergyCorrector=true;
225 }
226
227 #ifdef debug_G4ExcitedStringDecay
228 G4cout<<"NeedEnergyCorrection yes/no
229 #endif
230
231 // clean up
232 delete generatedKineticTracks;
233 success=true;
234 }
235
236 if ( NeedEnergyCorrector ) success=EnergyAnd
237 } while (!success && (attempts < 100)); /*
238
239 for ( unsigned int aTrack=0; aTrack<theResul
240 {
241 Ptmp=(*theResult)[aTrack]->Get4Momentum();
242 Ptmp.transform( toLab);
243 (*theResult)[aTrack]->Set4Momentum(Ptmp);
244 }
245
246 #ifdef debug_G4ExcitedStringDecay
247 G4cout<<"End of the strings fragmentation (G
248
249 G4LorentzVector KTsum1(0.,0.,0.,0.);
250
251 for ( unsigned int aTrack=0; aTrack<theResul
252 {
253 G4cout << " corrected tracks .. " << (*t
254 <<" " << (*theResult)[aTrack]->Get4Mome
255 <<" " << (*theResult)[aTrack]->Get4Mome
256 KTsum1+= (*theResult)[aTrack]->Get4Momen
257 }
258
259 G4cout << "Needcorrector/success " << NeedEn
260 << ", Corrected total 4 momentum " <
261 if ( ! success ) G4cout << "failed to correc
262
263 G4cout<<"End of the Hadronization (G4Excited
264 #endif
265
266 if (!success)
267 {
268 if (theResult->size() != 0)
269 {
270 std::for_each(theResult->begin() , theRe
271 theResult->clear();
272 delete theResult; theResult=0;
273 }
274 for ( unsigned int astring=0; astring < th
275 // for ( unsigned int astring=0; astring <
276 {
277 if ( theStrings->operator[](astring)->Is
278 {
279 Ptmp=theStrings->operator[](astring)->
280 Ptmp.transform( toLab);
281 theStrings->operator[](astring)->GetLe
282
283 Ptmp=theStrings->operator[](astring)->
284 Ptmp.transform( toLab);
285 theStrings->operator[](astring)->GetRi
286 }
287 else
288 {
289 Ptmp=theStrings->operator[](astring)->
290 Ptmp.transform( toLab);
291 theStrings->operator[](astring)->GetKi
292 }
293 }
294 }
295 return theResult;
296 }
297
298
299 G4bool G4ExcitedStringDecay::
300 EnergyAndMomentumCorrector(G4KineticTrackVecto
301 {
302 const int nAttemptScale = 500;
303 const double ErrLimit = 1.E-5;
304 if (Output->empty()) return TRUE;
305 G4LorentzVector SumMom;
306 G4double SumMass = 0;
307 G4double TotalCollisionMass = Total
308
309 std::vector<G4double> HadronMass; G4double
310
311 #ifdef debug_G4ExcitedStringCorr
312 G4cout<<G4endl<<"EnergyAndMomentumCorrecto
313 #endif
314 // Calculate sum hadron 4-momenta and summ
315 unsigned int cHadron;
316 for (cHadron = 0; cHadron < Output->size()
317 {
318 SumMom += Output->operator[](cHadron)
319 HadronM=Output->operator[](cHadron)->G
320 SumMass += Output->operator[](cHadron)
321 }
322
323 #ifdef debug_G4ExcitedStringCorr
324 G4cout<<"Sum part mom "<<SumMom<<" "<<SumM
325 <<"Sum str mom "<<TotalCollisionMom
326 G4cout<<"SumMass TotalCollisionMass "<<Sum
327 #endif
328
329 // Cannot correct a single particle
330 if (Output->size() < 2) return FALSE;
331
332 if (SumMass > TotalCollisionMass) return F
333 SumMass = SumMom.m2();
334 if (SumMass < 0) return FALSE;
335 SumMass = std::sqrt(SumMass);
336
337 // Compute c.m.s. hadron velocity and boos
338 // G4ThreeVector Beta = -SumMom.boostVecto
339 G4ThreeVector Beta = -TotalCollisionMom.bo
340 Output->Boost(Beta);
341
342 // Scale total c.m.s. hadron energy (hadro
343 // It should be equal interaction mass
344 G4double Scale = 1;
345 G4int cAttempt = 0;
346 G4double Sum = 0;
347 G4bool success = false;
348 for (cAttempt = 0; cAttempt < nAttemptScal
349 {
350 Sum = 0;
351 for (cHadron = 0; cHadron < Output->size
352 {
353 HadronM = HadronMass.at(cHadron);
354 G4LorentzVector HadronMom = Output->op
355 HadronMom.setVect(Scale*HadronMom.vect
356 G4double E = std::sqrt(HadronMom.vect(
357
358 HadronMom.setE(E);
359 Output->operator[](cHadron)->Set4Momen
360 Sum += E;
361 }
362 Scale = TotalCollisionMass/Sum;
363 #ifdef debug_G4ExcitedStringCorr
364 G4cout << "Scale-1=" << Scale -1
365 << ", TotalCollisionMass=" << To
366 << ", Sum=" << Sum
367 << G4endl;
368 #endif
369 if (std::fabs(Scale - 1) <= ErrLimit)
370 {
371 success = true;
372 break;
373 }
374 }
375
376 #ifdef debug_G4ExcitedStringCorr
377 if (!success)
378 {
379 G4cout << "G4ExcitedStringDecay::EnergyA
380 G4cout << " Scale not unity at end of
381 G4cout << " Number of secondaries: " <
382 G4cout << " Wanted total energy: " <<
383 G4cout << " Increase number of attempt
384 // throw G4HadronicException(__FILE__, _
385 }
386 #endif
387
388 // Compute c.m.s. interaction velocity and
389 Beta = TotalCollisionMom.boostVector();
390 Output->Boost(Beta);
391
392 return success;
393 }
394
395