Geant4  10.03.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4CollisionOutput.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // $Id: G4CollisionOutput.cc 71954 2013-06-29 04:40:40Z mkelsey $
27 //
28 // 20100114 M. Kelsey -- Remove G4CascadeMomentum, use G4LorentzVector directly
29 // 20100309 M. Kelsey -- Introduced bug checking i3 for valid tuning pair
30 // 20100409 M. Kelsey -- Move non-inlinable code here out of .hh file, replace
31 // loop over push_back() with block insert().
32 // 20100418 M. Kelsey -- Add function to boost output lists to lab frame
33 // 20100520 M. Kelsey -- Add function to rotate Z axis, from G4Casc.Interface
34 // 20100620 M. Kelsey -- Add some diagnostics in setOnShell, simplify if's
35 // 20100630 M. Kelsey -- Use "getExcitationEnergyInGeV()" instead of ".001*"
36 // 20100701 M. Kelsey -- G4InuclNuclei now includes excitation energy as part
37 // of the reported mass and four-vector.
38 // 20100714 M. Kelsey -- Modify setOnShell() to avoid creating particles
39 // with negative kinetic energy.
40 // 20100715 M. Kelsey -- Add total charge and baryon number functions, and a
41 // combined "add()" function to put two of these together
42 // 20100924 M. Kelsey -- Use "OutgoingNuclei" name consistently, replacing
43 // old "TargetFragment". Add new (reusable) G4Fragment buffer
44 // and access functions for initial post-cascade processing.
45 // Move implementation of add() to .cc file.
46 // 20101019 M. Kelsey -- CoVerity report: unitialized constructor
47 // 20110214 M. Kelsey -- Follow G4InuclParticle::Model enumerator migration
48 // 20110225 M. Kelsey -- Add non-const functions to remove list elements
49 // 20110302 M. Kelsey -- Fix message in setOnShell() by moving ini_mom calc
50 // 20110307 M. Kelsey -- Need to discard existing ouput lists in trivialize()
51 // 20110311 M. Kelsey -- Include nuclear fragment in setOnShell balancing,
52 // including calculation of final-state momentum
53 // 20110519 M. Kelsey -- Drop unused "p2" variable from selectPairToTune()
54 // 20110801 M. Kelsey -- Use resize to avoid temporaries when copying from
55 // G4ReactionProductVector
56 // 20110922 M. Kelsey -- Follow G4InuclParticle::print(ostream&) migration,
57 // Add optional stream argument to printCollisionOutput
58 // 20121002 M. Kelsey -- Add strangeness calculation
59 // 20130628 M. Kelsey -- Support multiple recoil fragments (for G4Fissioner)
60 // 20141208 M. Kelsey -- Split function to do pair-wise "hard" tuning
61 
62 #include <algorithm>
63 
64 #include "G4CollisionOutput.hh"
65 #include "G4SystemOfUnits.hh"
66 #include "G4CascadParticle.hh"
67 #include "G4ParticleLargerEkin.hh"
68 #include "G4LorentzConvertor.hh"
69 #include "G4LorentzRotation.hh"
70 #include "G4LorentzVector.hh"
72 #include "G4ReactionProduct.hh"
73 #include "G4ThreeVector.hh"
74 
75 typedef std::vector<G4InuclElementaryParticle>::iterator particleIterator;
76 typedef std::vector<G4InuclNuclei>::iterator nucleiIterator;
77 typedef std::vector<G4Fragment>::iterator fragmentIterator;
78 
79 const G4Fragment G4CollisionOutput::emptyFragment; // All zeros
80 
81 
83  : verboseLevel(0), eex_rest(0), on_shell(false) {
84  if (verboseLevel > 1)
85  G4cout << " >>> G4CollisionOutput::G4CollisionOutput" << G4endl;
86 }
87 
88 
90 {
91  if (this != &right) {
92  verboseLevel = right.verboseLevel;
93  outgoingParticles = right.outgoingParticles;
94  outgoingNuclei = right.outgoingNuclei;
95  recoilFragments = right.recoilFragments;
96  eex_rest = right.eex_rest;
97  on_shell = right.on_shell;
98  }
99  return *this;
100 }
101 
103  outgoingNuclei.clear();
104  outgoingParticles.clear();
105  recoilFragments.clear();
106  eex_rest = 0.;
107  on_shell = false;
108 }
109 
110 
111 // Get requested recoil fragment from list, or return dummy version
112 
114  return ( (index >= 0 && index < numberOfFragments())
115  ? recoilFragments[index] : emptyFragment);
116 }
117 
118 
119 // Merge two complete objects
120 
122  addOutgoingParticles(right.outgoingParticles);
123  addOutgoingNuclei(right.outgoingNuclei);
124  recoilFragments = right.recoilFragments; // Replace, don't combine
125  eex_rest = 0.;
126  on_shell = false;
127 }
128 
129 
130 // Append to lists
131 
132 void G4CollisionOutput::addOutgoingParticles(const std::vector<G4InuclElementaryParticle>& particles) {
133  outgoingParticles.insert(outgoingParticles.end(),
134  particles.begin(), particles.end());
135 }
136 
137 void G4CollisionOutput::addOutgoingNuclei(const std::vector<G4InuclNuclei>& nuclea) {
138  outgoingNuclei.insert(outgoingNuclei.end(), nuclea.begin(), nuclea.end());
139 }
140 
141 // These are primarily for G4IntraNucleiCascader internal checks
142 
144  addOutgoingParticle(cparticle.getParticle());
145 }
146 
147 void G4CollisionOutput::addOutgoingParticles(const std::vector<G4CascadParticle>& cparticles) {
148  for (unsigned i=0; i<cparticles.size(); i++)
149  addOutgoingParticle(cparticles[i]);
150 }
151 
152 // This comes from PreCompound de-excitation, both particles and nuclei
153 
155  if (!rproducts) return; // Sanity check, no error if null
156 
157  if (verboseLevel) {
158  G4cout << " >>> G4CollisionOutput::addOutgoingParticles(G4RPVector)"
159  << G4endl;
160  }
161 
162  G4ReactionProductVector::const_iterator j;
163  for (j=rproducts->begin(); j!=rproducts->end(); ++j) {
164  const G4ParticleDefinition* pd = (*j)->GetDefinition();
166 
167  // FIXME: Momentum returned by value; extra copying here!
168  G4LorentzVector mom((*j)->GetMomentum(), (*j)->GetTotalEnergy());
169  mom /= GeV; // Convert from GEANT4 to Bertini units
170 
171  if (verboseLevel>1)
172  G4cout << " Processing " << pd->GetParticleName() << " (" << type
173  << "), momentum " << mom << " GeV" << G4endl;
174 
175  // Nucleons and nuclei are jumbled together in the list
176  // NOTE: Resize and set/fill avoid unnecessary temporary copies
177  if (type) {
178  outgoingParticles.resize(numberOfOutgoingParticles()+1);
179  outgoingParticles.back().fill(mom, pd, G4InuclParticle::PreCompound);
180 
181  if (verboseLevel>1) G4cout << outgoingParticles.back() << G4endl;
182  } else {
183  outgoingNuclei.resize(numberOfOutgoingNuclei()+1);
184  outgoingNuclei.back().fill(mom,pd->GetAtomicMass(),pd->GetAtomicNumber(),
186 
187  if (verboseLevel>1) G4cout << outgoingNuclei.back() << G4endl;
188  }
189  }
190 }
191 
192 
193 // Removing elements from lists by index
194 
196  if (index >= 0 && index < numberOfOutgoingParticles())
197  outgoingParticles.erase(outgoingParticles.begin()+(size_t)index);
198 }
199 
201  if (index >= 0 && index < numberOfOutgoingNuclei())
202  outgoingNuclei.erase(outgoingNuclei.begin()+(size_t)index);
203 }
204 
205 // Remove elements from list by reference, or by value comparison
206 
209  std::find(outgoingParticles.begin(), outgoingParticles.end(), particle);
210  if (pos != outgoingParticles.end()) outgoingParticles.erase(pos);
211 }
212 
215  std::find(outgoingNuclei.begin(), outgoingNuclei.end(), nuclei);
216  if (pos != outgoingNuclei.end()) outgoingNuclei.erase(pos);
217 }
218 
219 // Remove specified recoil fragment(s) from buffer
220 
222  if (index < 0) recoilFragments.clear();
223  else if (index < numberOfFragments())
224  recoilFragments.erase(recoilFragments.begin()+(size_t)index);
225 }
226 
227 
228 // Kinematics of final state, for recoil and conservation checks
229 
231  if (verboseLevel > 1)
232  G4cout << " >>> G4CollisionOutput::getTotalOutputMomentum" << G4endl;
233 
234  G4LorentzVector tot_mom;
235  G4int i(0);
236  for(i=0; i < numberOfOutgoingParticles(); i++) {
237  tot_mom += outgoingParticles[i].getMomentum();
238  }
239  for(i=0; i < numberOfOutgoingNuclei(); i++) {
240  tot_mom += outgoingNuclei[i].getMomentum();
241  }
242  for(i=0; i < numberOfFragments(); i++) {
243  tot_mom += recoilFragments[i].GetMomentum()/GeV; // Need Bertini units!
244  }
245 
246  return tot_mom;
247 }
248 
250  if (verboseLevel > 1)
251  G4cout << " >>> G4CollisionOutput::getTotalCharge" << G4endl;
252 
253  G4int charge = 0;
254  G4int i(0);
255  for(i=0; i < numberOfOutgoingParticles(); i++) {
256  charge += G4int(outgoingParticles[i].getCharge());
257  }
258  for(i=0; i < numberOfOutgoingNuclei(); i++) {
259  charge += G4int(outgoingNuclei[i].getCharge());
260  }
261  for(i=0; i < numberOfFragments(); i++) {
262  charge += recoilFragments[i].GetZ_asInt();
263  }
264 
265  return charge;
266 }
267 
269  if (verboseLevel > 1)
270  G4cout << " >>> G4CollisionOutput::getTotalBaryonNumber" << G4endl;
271 
272  G4int baryon = 0;
273  G4int i(0);
274  for(i=0; i < numberOfOutgoingParticles(); i++) {
275  baryon += outgoingParticles[i].baryon();
276  }
277  for(i=0; i < numberOfOutgoingNuclei(); i++) {
278  baryon += G4int(outgoingNuclei[i].getA());
279  }
280  for(i=0; i < numberOfFragments(); i++) {
281  baryon += recoilFragments[i].GetA_asInt();
282  }
283 
284  return baryon;
285 }
286 
288  if (verboseLevel > 1)
289  G4cout << " >>> G4CollisionOutput::getTotalStrangeness" << G4endl;
290 
291  G4int strange = 0;
292  G4int i(0);
293  for(i=0; i < numberOfOutgoingParticles(); i++) {
294  strange += outgoingParticles[i].getStrangeness();
295  }
296 
297  return strange;
298 }
299 
300 
301 void G4CollisionOutput::printCollisionOutput(std::ostream& os) const {
302  os << " Output: " << G4endl
303  << " Outgoing Particles: " << numberOfOutgoingParticles() << G4endl;
304 
305  G4int i(0);
306  for(i=0; i < numberOfOutgoingParticles(); i++)
307  os << outgoingParticles[i] << G4endl;
308 
309  os << " Outgoing Nuclei: " << numberOfOutgoingNuclei() << G4endl;
310  for(i=0; i < numberOfOutgoingNuclei(); i++)
311  os << outgoingNuclei[i] << G4endl;
312  for(i=0; i < numberOfFragments(); i++)
313  os << recoilFragments[i] << G4endl;
314 }
315 
316 
317 // Boost particles and fragment to LAB -- "convertor" must already be configured
318 
320  if (verboseLevel > 1)
321  G4cout << " >>> G4CollisionOutput::boostToLabFrame" << G4endl;
322 
323  particleIterator ipart = outgoingParticles.begin();
324  for(; ipart != outgoingParticles.end(); ipart++) {
325  ipart->setMomentum(boostToLabFrame(ipart->getMomentum(), convertor));
326  }
327 
328  std::sort(outgoingParticles.begin(), outgoingParticles.end(),
330 
331  nucleiIterator inuc = outgoingNuclei.begin();
332  for (; inuc != outgoingNuclei.end(); inuc++) {
333  inuc->setMomentum(boostToLabFrame(inuc->getMomentum(), convertor));
334  }
335 
336  // NOTE: Fragment momentum must be converted to and from Bertini units
337  G4LorentzVector fmom;
338  fragmentIterator ifrag = recoilFragments.begin();
339  for (; ifrag != recoilFragments.end(); ifrag++) {
340  fmom = ifrag->GetMomentum() / GeV;
341  ifrag->SetMomentum(boostToLabFrame(fmom, convertor)*GeV);
342  }
343 }
344 
345 // Pass by value to allow direct (internal) manipulation
348  const G4LorentzConvertor& convertor) const {
349  if (convertor.reflectionNeeded()) mom.setZ(-mom.z());
350  mom = convertor.rotate(mom);
351  mom = convertor.backToTheLab(mom);
352 
353  return mom;
354 }
355 
356 // Apply LorentzRotation to all particles in event
357 
359  if (verboseLevel > 1)
360  G4cout << " >>> G4CollisionOutput::rotateEvent" << G4endl;
361 
362  particleIterator ipart = outgoingParticles.begin();
363  for(; ipart != outgoingParticles.end(); ipart++)
364  ipart->setMomentum(ipart->getMomentum()*=rotate);
365 
366  nucleiIterator inuc = outgoingNuclei.begin();
367  for (; inuc != outgoingNuclei.end(); inuc++)
368  inuc->setMomentum(inuc->getMomentum()*=rotate);
369 
370  fragmentIterator ifrag = recoilFragments.begin();
371  for (; ifrag != recoilFragments.end(); ifrag++) {
372  G4LorentzVector mom = ifrag->GetMomentum(); // Need copy
373  ifrag->SetMomentum(mom*=rotate);
374  }
375 }
376 
377 
380  if (verboseLevel > 1)
381  G4cout << " >>> G4CollisionOutput::trivialize" << G4endl;
382 
383  reset(); // Discard existing output, replace with bullet/target
384 
385  if (G4InuclNuclei* nuclei_target = dynamic_cast<G4InuclNuclei*>(target)) {
386  outgoingNuclei.push_back(*nuclei_target);
387  } else {
388  G4InuclElementaryParticle* particle =
389  dynamic_cast<G4InuclElementaryParticle*>(target);
390  outgoingParticles.push_back(*particle);
391  }
392 
393  if (G4InuclNuclei* nuclei_bullet = dynamic_cast<G4InuclNuclei*>(bullet)) {
394  outgoingNuclei.push_back(*nuclei_bullet);
395  } else {
396  G4InuclElementaryParticle* particle =
397  dynamic_cast<G4InuclElementaryParticle*>(bullet);
398  outgoingParticles.push_back(*particle);
399  }
400 }
401 
402 
405  if (verboseLevel > 1)
406  G4cout << " >>> G4CollisionOutput::setOnShell" << G4endl;
407 
408  const G4double accuracy = 0.00001; // momentum concerves at the level of 10 keV
409 
410  on_shell = false;
411 
412  G4LorentzVector ini_mom = bullet->getMomentum();
413  G4LorentzVector momt = target->getMomentum();
415  if(verboseLevel > 2){
416  G4cout << " bullet momentum = " << ini_mom.e() <<", "<< ini_mom.x() <<", "<< ini_mom.y()<<", "<< ini_mom.z()<<G4endl;
417  G4cout << " target momentum = " << momt.e()<<", "<< momt.x()<<", "<< momt.y()<<", "<< momt.z()<<G4endl;
418  G4cout << " Fstate momentum = " << out_mom.e()<<", "<< out_mom.x()<<", "<< out_mom.y()<<", "<< out_mom.z()<<G4endl;
419  }
420 
421  ini_mom += momt;
422 
423  mom_non_cons = ini_mom - out_mom;
424  G4double pnc = mom_non_cons.rho();
425  G4double enc = mom_non_cons.e();
426 
428 
429  if(verboseLevel > 2) {
431  G4cout << " momentum non conservation: " << G4endl
432  << " e " << enc << " p " << pnc << G4endl
433  << " remaining exitation " << eex_rest << G4endl;
434  }
435 
436  if (std::fabs(enc) <= accuracy && pnc <= accuracy) {
437  on_shell = true;
438  return;
439  }
440 
441  // Adjust "last" particle's four-momentum to balance event
442  // ONLY adjust particles with sufficient e or p to remain physical!
443 
444  if (verboseLevel > 2) G4cout << " re-balancing four-momenta" << G4endl;
445 
447  G4int nnuc = numberOfOutgoingNuclei();
448  G4int nfrag = numberOfFragments();
449 
450  G4LorentzVector last_mom; // Buffer to reduce memory churn
451 
452  if (npart > 0) {
453  for (G4int ip=npart-1; ip>=0; ip--) {
454  if (outgoingParticles[ip].getKineticEnergy()+enc > 0.) {
455  last_mom = outgoingParticles[ip].getMomentum();
456  last_mom += mom_non_cons;
457  outgoingParticles[ip].setMomentum(last_mom);
458  break;
459  }
460  }
461  } else if (nnuc > 0) {
462  for (G4int in=nnuc-1; in>=0; in--) {
463  if (outgoingNuclei[in].getKineticEnergy()+enc > 0.) {
464  last_mom = outgoingNuclei[in].getMomentum();
465  last_mom += mom_non_cons;
466  outgoingNuclei[in].setMomentum(last_mom);
467  break;
468  }
469  }
470  } else if (nfrag > 0) {
471  for (G4int ifr=nfrag-1; ifr>=0; ifr--) {
472  // NOTE: G4Fragment does not use Bertini units!
473  last_mom = recoilFragments[ifr].GetMomentum()/GeV;
474  if ((last_mom.e()-last_mom.m())+enc > 0.) {
475  last_mom += mom_non_cons;
476  recoilFragments[ifr].SetMomentum(last_mom*GeV);
477  break;
478  }
479  }
480  }
481 
482  // Recompute momentum non-conservation parameters
483  out_mom = getTotalOutputMomentum();
484  mom_non_cons = ini_mom - out_mom;
485  pnc = mom_non_cons.rho();
486  enc = mom_non_cons.e();
487 
488  if(verboseLevel > 2){
490  G4cout << " momentum non conservation after (1): " << G4endl
491  << " e " << enc << " p " << pnc << G4endl;
492  }
493 
494  // Can energy be balanced just with nuclear excitation?
495  G4bool need_hard_tuning = true;
496 
497  G4double encMeV = mom_non_cons.e() / GeV; // Excitation below is in MeV
498  if (nfrag > 0) {
499  for (G4int i=0; i < nfrag; i++) {
500  G4double eex = recoilFragments[0].GetExcitationEnergy();
501  if (eex > 0.0 && eex + encMeV >= 0.0) {
502  // NOTE: G4Fragment doesn't have function to change excitation energy
503  // ==> theRecoilFragment.SetExcitationEnergy(eex+encMeV);
504 
505  G4LorentzVector fragMom = recoilFragments[i].GetMomentum();
506  G4double newMass = fragMom.m() + encMeV; // .m() includes eex
507  fragMom.setVectM(fragMom.vect(), newMass);
508  need_hard_tuning = false;
509  break;
510  }
511  }
512  } else if (nnuc > 0) {
513  for (G4int i=0; i < nnuc; i++) {
514  G4double eex = outgoingNuclei[i].getExitationEnergy();
515  if (eex > 0.0 && eex + encMeV >= 0.0) {
516  outgoingNuclei[i].setExitationEnergy(eex+encMeV);
517  need_hard_tuning = false;
518  break;
519  }
520  }
521  if (need_hard_tuning && encMeV > 0.) {
522  outgoingNuclei[0].setExitationEnergy(encMeV);
523  need_hard_tuning = false;
524  }
525  }
526 
527  if (!need_hard_tuning) {
528  on_shell = true;
529  return;
530  }
531 
532  // Momentum (hard) tuning required for energy conservation
533  if (verboseLevel > 2)
534  G4cout << " trying hard (particle-pair) tuning" << G4endl;
535 
536  /*****
537  // Hard tuning of quasielastic particle against nucleus
538  if (npart == 1) {
539  if (verboseLevel > 2)
540  G4cout << " tuning particle 0 against recoil nucleus" << G4endl;
541 
542  G4LorentzVector mom1 = outgoingParticles[0].getMomentum();
543  G4LorentzVector mom2 = recoilFragments[0].GetMomentum()/GeV;
544 
545  // Preserve momentum direction of outgoing particle
546  G4ThreeVector pdir = mom1.vect().unit();
547 
548  // Two-body kinematics (nucleon against nucleus) in overall CM system
549  G4double ecmsq = ini_mom.m2();
550  G4double a = 0.5 * (ecmsq - mom1.m2() - mom2.m2());
551  G4double pmod = std::sqrt((a*a - mom1.m2()*mom2.m2()) / ecmsq);
552 
553  mom1.setVectM(pdir*pmod, mom1.m());
554  mom2.setVectM(-pdir*pmod, mom2.m());
555 
556  // Boost CM momenta back into lab frame, assign back to particles
557  mom1.boost(-ini_mom.boostVector());
558  mom2.boost(-ini_mom.boostVector());
559 
560  outgoingParticles[0].setMomentum(mom1);
561  recoilFragments[0].SetMomentum(mom2*GeV);
562  } else { // Hard tuning using pair of outgoing particles
563  *****/
564  std::pair<std::pair<G4int, G4int>, G4int> tune_par = selectPairToTune(enc);
565  std::pair<G4int, G4int> tune_particles = tune_par.first;
566  G4int mom_ind = tune_par.second;
567 
568  G4bool tuning_possible =
569  (tune_particles.first >= 0 && tune_particles.second >= 0 &&
570  mom_ind >= G4LorentzVector::X);
571 
572  if (!tuning_possible) {
573  if (verboseLevel > 2) G4cout << " tuning impossible " << G4endl;
574  return;
575  }
576 
577  if(verboseLevel > 2) {
578  G4cout << " p1 " << tune_particles.first << " p2 " << tune_particles.second
579  << " ind " << mom_ind << G4endl;
580  }
581 
582  G4LorentzVector mom1 = outgoingParticles[tune_particles.first].getMomentum();
583  G4LorentzVector mom2 = outgoingParticles[tune_particles.second].getMomentum();
584 
585  if (tuneSelectedPair(mom1, mom2, mom_ind)) {
586  outgoingParticles[tune_particles.first ].setMomentum(mom1);
587  outgoingParticles[tune_particles.second].setMomentum(mom2);
588  }
589  else return;
590  //***** }
591 
592  out_mom = getTotalOutputMomentum();
593  std::sort(outgoingParticles.begin(), outgoingParticles.end(), G4ParticleLargerEkin());
594  mom_non_cons = ini_mom - out_mom;
595  pnc = mom_non_cons.rho();
596  enc = mom_non_cons.e();
597 
598  on_shell = (std::fabs(enc) < accuracy || pnc < accuracy);
599 
600  if(verboseLevel > 2) {
601  G4cout << " momentum non conservation tuning: " << G4endl
602  << " e " << enc << " p " << pnc
603  << (on_shell?" success":" FAILURE") << G4endl;
604  }
605 }
606 
607 
608 // Returns excitation energy in GeV
609 
611  eex_rest = 0.;
612  G4int i(0);
613  for (i=0; i < numberOfOutgoingNuclei(); i++)
614  eex_rest += outgoingNuclei[i].getExitationEnergyInGeV();
615  for (i=0; i < numberOfFragments(); i++)
616  eex_rest += recoilFragments[i].GetExcitationEnergy() / GeV;
617 }
618 
619 
620 std::pair<std::pair<G4int, G4int>, G4int>
621 G4CollisionOutput::selectPairToTune(G4double de) const {
622  if (verboseLevel > 2)
623  G4cout << " >>> G4CollisionOutput::selectPairToTune" << G4endl;
624 
625  std::pair<G4int, G4int> tup(-1, -1);
626  G4int i3 = -1;
627  std::pair<std::pair<G4int, G4int>, G4int> tuner(tup, i3); // Set invalid
628 
629  if (outgoingParticles.size() < 2) return tuner; // Nothing to do
630 
631  G4int ibest1 = -1;
632  G4int ibest2 = -1;
633  G4double pbest = 0.0;
634  G4double pcut = 0.3 * std::sqrt(1.88 * std::fabs(de));
635  G4double p1 = 0.0;
636 
637  for (G4int i = 0; i < G4int(outgoingParticles.size())-1; i++) {
638  G4LorentzVector mom1 = outgoingParticles[i].getMomentum();
639 
640  for (G4int j = i+1; j < G4int(outgoingParticles.size()); j++) {
641  G4LorentzVector mom2 = outgoingParticles[j].getMomentum();
642 
643  for (G4int l = G4LorentzVector::X; l<=G4LorentzVector::Z; l++) {
644  if (mom1[l]*mom2[l]<0.0) {
645  if (std::fabs(mom1[l])>pcut && std::fabs(mom2[l])>pcut) {
646  G4double psum = std::fabs(mom1[l]) + std::fabs(mom2[l]);
647 
648  if (psum > pbest) {
649  ibest1 = i;
650  ibest2 = j;
651  i3 = l;
652  p1 = mom1[l];
653  pbest = psum;
654  } // psum > pbest
655  } // mom1 and mom2 > pcut
656  } // mom1 ~ -mom2
657  } // for (G4int l ...
658  } // for (G4int j ...
659  } // for (G4int i ...
660 
661  if (i3 < 0) return tuner;
662 
663  tuner.second = i3; // Momentum axis for tuning
664 
665  // NOTE: Sign of de determines order for special case of p1==0.
666  if (de > 0.0) {
667  tuner.first.first = (p1>0.) ? ibest1 : ibest2;
668  tuner.first.second = (p1>0.) ? ibest2 : ibest1;
669  } else {
670  tuner.first.first = (p1<0.) ? ibest2 : ibest1;
671  tuner.first.second = (p1<0.) ? ibest1 : ibest2;
672  }
673 
674  return tuner;
675 }
676 
677 
678 G4bool G4CollisionOutput::tuneSelectedPair(G4LorentzVector& mom1,
679  G4LorentzVector& mom2,
680  G4int mom_ind) const {
681  if (verboseLevel > 2)
682  G4cout << " >>> G4CollisionOutput::tuneSelectedPair" << G4endl;
683 
684  G4double newE12 = mom1.e() + mom2.e() + mom_non_cons.e();
685  G4double R = 0.5 * (newE12*newE12 + mom2.e()*mom2.e() - mom1.e()*mom1.e()) / newE12;
686  G4double Q = -(mom1[mom_ind] + mom2[mom_ind]) / newE12;
687  G4double UDQ = 1.0 / (Q*Q - 1.0);
688  G4double W = (R * Q + mom2[mom_ind]) * UDQ;
689  G4double V = (mom2.e()*mom2.e() - R*R) * UDQ;
690  G4double DET = W*W + V;
691 
692  if (DET < 0.0) {
693  if (verboseLevel > 2) G4cout << " DET < 0 : tuning failed" << G4endl;
694  return false;
695  }
696 
697  // Tuning allowed only for non-negative determinant
698  G4double x1 = -(W + std::sqrt(DET));
699  G4double x2 = -(W - std::sqrt(DET));
700 
701  // choose the appropriate solution
702  G4bool xset = false;
703  G4double x = 0.0;
704 
705  if (mom_non_cons.e() > 0.0) { // x has to be > 0.0
706  if (x1 > 0.0) {
707  if (R + Q * x1 >= 0.0) {
708  x = x1;
709  xset = true;
710  }
711  }
712  if (!xset && x2 > 0.0) {
713  if (R + Q * x2 >= 0.0) {
714  x = x2;
715  xset = true;
716  }
717  }
718  } else {
719  if (x1 < 0.0) {
720  if (R + Q * x1 >= 0.) {
721  x = x1;
722  xset = true;
723  }
724  }
725  if (!xset && x2 < 0.0) {
726  if (R + Q * x2 >= 0.0) {
727  x = x2;
728  xset = true;
729  }
730  }
731  } // if(mom_non_cons.e() > 0.0)
732 
733  if (!xset) {
734  if (verboseLevel > 2) G4cout << " no appropriate solution found" << G4endl;
735  return false;
736  }
737 
738  mom1[mom_ind] += x; // retune momentums
739  mom2[mom_ind] -= x;
740  return true;
741 }
void trivialise(G4InuclParticle *bullet, G4InuclParticle *target)
G4CollisionOutput & operator=(const G4CollisionOutput &right)
void addOutgoingNuclei(const std::vector< G4InuclNuclei > &nuclea)
const XML_Char * target
Definition: expat.h:268
G4LorentzVector getMomentum() const
void removeRecoilFragment(G4int index=-1)
void printCollisionOutput(std::ostream &os=G4cout) const
static double Q[]
void addOutgoingParticle(const G4InuclElementaryParticle &particle)
tuple x
Definition: test.py:50
G4LorentzVector rotate(const G4LorentzVector &mom) const
G4LorentzVector backToTheLab(const G4LorentzVector &mom) const
int G4int
Definition: G4Types.hh:78
G4LorentzVector getTotalOutputMomentum() const
const G4InuclElementaryParticle & getParticle() const
const G4String & GetParticleName() const
G4int GetAtomicNumber() const
std::vector< G4ReactionProduct * > G4ReactionProductVector
void setVectM(const Hep3Vector &spatial, double mass)
void add(const G4CollisionOutput &right)
Hep3Vector vect() const
G4GLOB_DLL std::ostream G4cout
bool G4bool
Definition: G4Types.hh:79
G4int getTotalStrangeness() const
G4int numberOfOutgoingParticles() const
void boostToLabFrame(const G4LorentzConvertor &convertor)
void removeOutgoingParticle(G4int index)
double rho() const
void rotateEvent(const G4LorentzRotation &rotate)
G4int GetAtomicMass() const
G4int numberOfOutgoingNuclei() const
G4bool reflectionNeeded() const
void addOutgoingParticles(const std::vector< G4InuclElementaryParticle > &particles)
G4int getTotalBaryonNumber() const
static constexpr double GeV
Definition: G4SIunits.hh:217
std::vector< G4InuclNuclei >::const_iterator nucleiIterator
G4int getTotalCharge() const
std::vector< G4InuclElementaryParticle >::iterator particleIterator
Definition: G4BigBanger.cc:65
std::vector< G4Fragment >::iterator fragmentIterator
#define G4endl
Definition: G4ios.hh:61
G4int numberOfFragments() const
void removeOutgoingNucleus(G4int index)
const G4Fragment & getRecoilFragment(G4int index=0) const
double G4double
Definition: G4Types.hh:76
static const G4double pos
void setOnShell(G4InuclParticle *bullet, G4InuclParticle *target)