Geant4  10.00.p02
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 
61 #include <algorithm>
62 
63 #include "G4CollisionOutput.hh"
64 #include "G4SystemOfUnits.hh"
65 #include "G4CascadParticle.hh"
66 #include "G4ParticleLargerEkin.hh"
67 #include "G4LorentzConvertor.hh"
68 #include "G4LorentzRotation.hh"
69 #include "G4LorentzVector.hh"
71 #include "G4ReactionProduct.hh"
72 #include "G4ThreeVector.hh"
73 
74 typedef std::vector<G4InuclElementaryParticle>::iterator particleIterator;
75 typedef std::vector<G4InuclNuclei>::iterator nucleiIterator;
76 typedef std::vector<G4Fragment>::iterator fragmentIterator;
77 
79 
80 
82  : verboseLevel(0), eex_rest(0), on_shell(false) {
83  if (verboseLevel > 1)
84  G4cout << " >>> G4CollisionOutput::G4CollisionOutput" << G4endl;
85 }
86 
87 
89 {
90  if (this != &right) {
91  verboseLevel = right.verboseLevel;
95  eex_rest = right.eex_rest;
96  on_shell = right.on_shell;
97  }
98  return *this;
99 }
100 
102  outgoingNuclei.clear();
103  outgoingParticles.clear();
104  recoilFragments.clear();
105  eex_rest = 0.;
106  on_shell = false;
107 }
108 
109 
110 // Get requested recoil fragment from list, or return dummy version
111 
113  return ( (index >= 0 && index < numberOfFragments())
114  ? recoilFragments[index] : emptyFragment);
115 }
116 
117 
118 // Merge two complete objects
119 
123  recoilFragments = right.recoilFragments; // Replace, don't combine
124  eex_rest = 0.;
125  on_shell = false;
126 }
127 
128 
129 // Append to lists
130 
131 void G4CollisionOutput::addOutgoingParticles(const std::vector<G4InuclElementaryParticle>& particles) {
133  particles.begin(), particles.end());
134 }
135 
136 void G4CollisionOutput::addOutgoingNuclei(const std::vector<G4InuclNuclei>& nuclea) {
137  outgoingNuclei.insert(outgoingNuclei.end(), nuclea.begin(), nuclea.end());
138 }
139 
140 // These are primarily for G4IntraNucleiCascader internal checks
141 
143  addOutgoingParticle(cparticle.getParticle());
144 }
145 
146 void G4CollisionOutput::addOutgoingParticles(const std::vector<G4CascadParticle>& cparticles) {
147  for (unsigned i=0; i<cparticles.size(); i++)
148  addOutgoingParticle(cparticles[i]);
149 }
150 
151 // This comes from PreCompound de-excitation, both particles and nuclei
152 
154  if (!rproducts) return; // Sanity check, no error if null
155 
156  if (verboseLevel) {
157  G4cout << " >>> G4CollisionOutput::addOutgoingParticles(G4RPVector)"
158  << G4endl;
159  }
160 
161  G4ReactionProductVector::const_iterator j;
162  for (j=rproducts->begin(); j!=rproducts->end(); ++j) {
163  G4ParticleDefinition* pd = (*j)->GetDefinition();
165 
166  // FIXME: Momentum returned by value; extra copying here!
167  G4LorentzVector mom((*j)->GetMomentum(), (*j)->GetTotalEnergy());
168  mom /= GeV; // Convert from GEANT4 to Bertini units
169 
170  if (verboseLevel>1)
171  G4cout << " Processing " << pd->GetParticleName() << " (" << type
172  << "), momentum " << mom << " GeV" << G4endl;
173 
174  // Nucleons and nuclei are jumbled together in the list
175  // NOTE: Resize and set/fill avoid unnecessary temporary copies
176  if (type) {
178  outgoingParticles.back().fill(mom, pd, G4InuclParticle::PreCompound);
179 
180  if (verboseLevel>1) G4cout << outgoingParticles.back() << G4endl;
181  } else {
183  outgoingNuclei.back().fill(mom,pd->GetAtomicMass(),pd->GetAtomicNumber(),
185 
186  if (verboseLevel>1) G4cout << outgoingNuclei.back() << G4endl;
187  }
188  }
189 }
190 
191 
192 // Removing elements from lists by index
193 
195  if (index >= 0 && index < numberOfOutgoingParticles())
196  outgoingParticles.erase(outgoingParticles.begin()+(size_t)index);
197 }
198 
200  if (index >= 0 && index < numberOfOutgoingNuclei())
201  outgoingNuclei.erase(outgoingNuclei.begin()+(size_t)index);
202 }
203 
204 // Remove elements from list by reference, or by value comparison
205 
208  std::find(outgoingParticles.begin(), outgoingParticles.end(), particle);
209  if (pos != outgoingParticles.end()) outgoingParticles.erase(pos);
210 }
211 
214  std::find(outgoingNuclei.begin(), outgoingNuclei.end(), nuclei);
215  if (pos != outgoingNuclei.end()) outgoingNuclei.erase(pos);
216 }
217 
218 // Remove specified recoil fragment(s) from buffer
219 
221  if (index < 0) recoilFragments.clear();
222  else if (index < numberOfFragments())
223  recoilFragments.erase(recoilFragments.begin()+(size_t)index);
224 }
225 
226 
227 // Kinematics of final state, for recoil and conservation checks
228 
230  if (verboseLevel > 1)
231  G4cout << " >>> G4CollisionOutput::getTotalOutputMomentum" << G4endl;
232 
233  G4LorentzVector tot_mom;
234  G4int i(0);
235  for(i=0; i < numberOfOutgoingParticles(); i++) {
236  tot_mom += outgoingParticles[i].getMomentum();
237  }
238  for(i=0; i < numberOfOutgoingNuclei(); i++) {
239  tot_mom += outgoingNuclei[i].getMomentum();
240  }
241  for(i=0; i < numberOfFragments(); i++) {
242  tot_mom += recoilFragments[i].GetMomentum()/GeV; // Need Bertini units!
243  }
244 
245  return tot_mom;
246 }
247 
249  if (verboseLevel > 1)
250  G4cout << " >>> G4CollisionOutput::getTotalCharge" << G4endl;
251 
252  G4int charge = 0;
253  G4int i(0);
254  for(i=0; i < numberOfOutgoingParticles(); i++) {
255  charge += G4int(outgoingParticles[i].getCharge());
256  }
257  for(i=0; i < numberOfOutgoingNuclei(); i++) {
258  charge += G4int(outgoingNuclei[i].getCharge());
259  }
260  for(i=0; i < numberOfFragments(); i++) {
261  charge += recoilFragments[i].GetZ_asInt();
262  }
263 
264  return charge;
265 }
266 
268  if (verboseLevel > 1)
269  G4cout << " >>> G4CollisionOutput::getTotalBaryonNumber" << G4endl;
270 
271  G4int baryon = 0;
272  G4int i(0);
273  for(i=0; i < numberOfOutgoingParticles(); i++) {
274  baryon += outgoingParticles[i].baryon();
275  }
276  for(i=0; i < numberOfOutgoingNuclei(); i++) {
277  baryon += G4int(outgoingNuclei[i].getA());
278  }
279  for(i=0; i < numberOfFragments(); i++) {
280  baryon += recoilFragments[i].GetA_asInt();
281  }
282 
283  return baryon;
284 }
285 
287  if (verboseLevel > 1)
288  G4cout << " >>> G4CollisionOutput::getTotalStrangeness" << G4endl;
289 
290  G4int strange = 0;
291  G4int i(0);
292  for(i=0; i < numberOfOutgoingParticles(); i++) {
293  strange += outgoingParticles[i].getStrangeness();
294  }
295 
296  return strange;
297 }
298 
299 
300 void G4CollisionOutput::printCollisionOutput(std::ostream& os) const {
301  os << " Output: " << G4endl
302  << " Outgoing Particles: " << numberOfOutgoingParticles() << G4endl;
303 
304  G4int i(0);
305  for(i=0; i < numberOfOutgoingParticles(); i++)
306  os << outgoingParticles[i] << G4endl;
307 
308  os << " Outgoing Nuclei: " << numberOfOutgoingNuclei() << G4endl;
309  for(i=0; i < numberOfOutgoingNuclei(); i++)
310  os << outgoingNuclei[i] << G4endl;
311  for(i=0; i < numberOfFragments(); i++)
312  os << recoilFragments[i] << G4endl;
313 }
314 
315 
316 // Boost particles and fragment to LAB -- "convertor" must already be configured
317 
319  if (verboseLevel > 1)
320  G4cout << " >>> G4CollisionOutput::boostToLabFrame" << G4endl;
321 
322  particleIterator ipart = outgoingParticles.begin();
323  for(; ipart != outgoingParticles.end(); ipart++) {
324  ipart->setMomentum(boostToLabFrame(ipart->getMomentum(), convertor));
325  }
326 
327  std::sort(outgoingParticles.begin(), outgoingParticles.end(),
329 
330  nucleiIterator inuc = outgoingNuclei.begin();
331  for (; inuc != outgoingNuclei.end(); inuc++) {
332  inuc->setMomentum(boostToLabFrame(inuc->getMomentum(), convertor));
333  }
334 
335  // NOTE: Fragment momentum must be converted to and from Bertini units
336  G4LorentzVector fmom;
337  fragmentIterator ifrag = recoilFragments.begin();
338  for (; ifrag != recoilFragments.end(); ifrag++) {
339  fmom = ifrag->GetMomentum() / GeV;
340  ifrag->SetMomentum(boostToLabFrame(fmom, convertor)*GeV);
341  }
342 }
343 
344 // Pass by value to allow direct (internal) manipulation
347  const G4LorentzConvertor& convertor) const {
348  if (convertor.reflectionNeeded()) mom.setZ(-mom.z());
349  mom = convertor.rotate(mom);
350  mom = convertor.backToTheLab(mom);
351 
352  return mom;
353 }
354 
355 // Apply LorentzRotation to all particles in event
356 
358  if (verboseLevel > 1)
359  G4cout << " >>> G4CollisionOutput::rotateEvent" << G4endl;
360 
361  particleIterator ipart = outgoingParticles.begin();
362  for(; ipart != outgoingParticles.end(); ipart++)
363  ipart->setMomentum(ipart->getMomentum()*=rotate);
364 
365  nucleiIterator inuc = outgoingNuclei.begin();
366  for (; inuc != outgoingNuclei.end(); inuc++)
367  inuc->setMomentum(inuc->getMomentum()*=rotate);
368 
369  fragmentIterator ifrag = recoilFragments.begin();
370  for (; ifrag != recoilFragments.end(); ifrag++) {
371  G4LorentzVector mom = ifrag->GetMomentum(); // Need copy
372  ifrag->SetMomentum(mom*=rotate);
373  }
374 }
375 
376 
378  G4InuclParticle* target) {
379  if (verboseLevel > 1)
380  G4cout << " >>> G4CollisionOutput::trivialize" << G4endl;
381 
382  reset(); // Discard existing output, replace with bullet/target
383 
384  if (G4InuclNuclei* nuclei_target = dynamic_cast<G4InuclNuclei*>(target)) {
385  outgoingNuclei.push_back(*nuclei_target);
386  } else {
387  G4InuclElementaryParticle* particle =
388  dynamic_cast<G4InuclElementaryParticle*>(target);
389  outgoingParticles.push_back(*particle);
390  }
391 
392  if (G4InuclNuclei* nuclei_bullet = dynamic_cast<G4InuclNuclei*>(bullet)) {
393  outgoingNuclei.push_back(*nuclei_bullet);
394  } else {
395  G4InuclElementaryParticle* particle =
396  dynamic_cast<G4InuclElementaryParticle*>(bullet);
397  outgoingParticles.push_back(*particle);
398  }
399 }
400 
401 
403  G4InuclParticle* target) {
404  if (verboseLevel > 1)
405  G4cout << " >>> G4CollisionOutput::setOnShell" << G4endl;
406 
407  const G4double accuracy = 0.00001; // momentum concerves at the level of 10 keV
408 
409  on_shell = false;
410 
411  G4LorentzVector ini_mom = bullet->getMomentum();
412  G4LorentzVector momt = target->getMomentum();
414  if(verboseLevel > 2){
415  G4cout << " bullet momentum = " << ini_mom.e() <<", "<< ini_mom.x() <<", "<< ini_mom.y()<<", "<< ini_mom.z()<<G4endl;
416  G4cout << " target momentum = " << momt.e()<<", "<< momt.x()<<", "<< momt.y()<<", "<< momt.z()<<G4endl;
417  G4cout << " Fstate momentum = " << out_mom.e()<<", "<< out_mom.x()<<", "<< out_mom.y()<<", "<< out_mom.z()<<G4endl;
418  }
419 
420  ini_mom += momt;
421  G4LorentzVector mom_non_cons = ini_mom - out_mom;
422 
423  G4double pnc = mom_non_cons.rho();
424  G4double enc = mom_non_cons.e();
425 
427 
428  if(verboseLevel > 2) {
430  G4cout << " momentum non conservation: " << G4endl
431  << " e " << enc << " p " << pnc << G4endl
432  << " remaining exitation " << eex_rest << G4endl;
433  }
434 
435  if (std::fabs(enc) <= accuracy && pnc <= accuracy) {
436  on_shell = true;
437  return;
438  }
439 
440  // Adjust "last" particle's four-momentum to balance event
441  // ONLY adjust particles with sufficient e or p to remain physical!
442 
443  if (verboseLevel > 2) G4cout << " re-balancing four-momenta" << G4endl;
444 
446  G4int nnuc = numberOfOutgoingNuclei();
447  G4int nfrag = numberOfFragments();
448 
449  G4LorentzVector last_mom; // Buffer to reduce memory churn
450 
451  if (npart > 0) {
452  for (G4int ip=npart-1; ip>=0; ip--) {
453  if (outgoingParticles[ip].getKineticEnergy()+enc > 0.) {
454  last_mom = outgoingParticles[ip].getMomentum();
455  last_mom += mom_non_cons;
456  outgoingParticles[ip].setMomentum(last_mom);
457  break;
458  }
459  }
460  } else if (nnuc > 0) {
461  for (G4int in=nnuc-1; in>=0; in--) {
462  if (outgoingNuclei[in].getKineticEnergy()+enc > 0.) {
463  last_mom = outgoingNuclei[in].getMomentum();
464  last_mom += mom_non_cons;
465  outgoingNuclei[in].setMomentum(last_mom);
466  break;
467  }
468  }
469  } else if (nfrag > 0) {
470  for (G4int ifr=nfrag-1; ifr>=0; ifr--) {
471  // NOTE: G4Fragment does not use Bertini units!
472  last_mom = recoilFragments[ifr].GetMomentum()/GeV;
473  if ((last_mom.e()-last_mom.m())+enc > 0.) {
474  last_mom += mom_non_cons;
475  recoilFragments[ifr].SetMomentum(last_mom*GeV);
476  break;
477  }
478  }
479  }
480 
481  // Recompute momentum non-conservation parameters
482  out_mom = getTotalOutputMomentum();
483  mom_non_cons = ini_mom - out_mom;
484  pnc = mom_non_cons.rho();
485  enc = mom_non_cons.e();
486 
487  if(verboseLevel > 2){
489  G4cout << " momentum non conservation after (1): " << G4endl
490  << " e " << enc << " p " << pnc << G4endl;
491  }
492 
493  // Can energy be balanced just with nuclear excitation?
494  G4bool need_hard_tuning = true;
495 
496  G4double encMeV = mom_non_cons.e() / GeV; // Excitation below is in MeV
497  if (nfrag > 0) {
498  for (G4int i=0; i < nfrag; i++) {
499  G4double eex = recoilFragments[0].GetExcitationEnergy();
500  if (eex > 0.0 && eex + encMeV >= 0.0) {
501  // NOTE: G4Fragment doesn't have function to change excitation energy
502  // ==> theRecoilFragment.SetExcitationEnergy(eex+encMeV);
503 
504  G4LorentzVector fragMom = recoilFragments[i].GetMomentum();
505  G4double newMass = fragMom.m() + encMeV; // .m() includes eex
506  fragMom.setVectM(fragMom.vect(), newMass);
507  need_hard_tuning = false;
508  break;
509  }
510  }
511  } else if (nnuc > 0) {
512  for (G4int i=0; i < nnuc; i++) {
513  G4double eex = outgoingNuclei[i].getExitationEnergy();
514  if (eex > 0.0 && eex + encMeV >= 0.0) {
515  outgoingNuclei[i].setExitationEnergy(eex+encMeV);
516  need_hard_tuning = false;
517  break;
518  }
519  }
520  if (need_hard_tuning && encMeV > 0.) {
521  outgoingNuclei[0].setExitationEnergy(encMeV);
522  need_hard_tuning = false;
523  }
524  }
525 
526  if (!need_hard_tuning) {
527  on_shell = true;
528  return;
529  }
530 
531  // Momentum (hard) tuning required for energy conservation
532  if (verboseLevel > 2)
533  G4cout << " trying hard (particle-pair) tuning" << G4endl;
534 
535  std::pair<std::pair<G4int, G4int>, G4int> tune_par = selectPairToTune(mom_non_cons.e());
536  std::pair<G4int, G4int> tune_particles = tune_par.first;
537  G4int mom_ind = tune_par.second;
538 
539  if(verboseLevel > 2) {
540  G4cout << " p1 " << tune_particles.first << " p2 " << tune_particles.second
541  << " ind " << mom_ind << G4endl;
542  }
543 
544  G4bool tuning_possible =
545  (tune_particles.first >= 0 && tune_particles.second >= 0 &&
546  mom_ind >= G4LorentzVector::X);
547 
548  if (!tuning_possible) {
549  if (verboseLevel > 2) G4cout << " tuning impossible " << G4endl;
550  return;
551  }
552 
553  G4LorentzVector mom1 = outgoingParticles[tune_particles.first].getMomentum();
554  G4LorentzVector mom2 = outgoingParticles[tune_particles.second].getMomentum();
555  G4double newE12 = mom1.e() + mom2.e() + mom_non_cons.e();
556  G4double R = 0.5 * (newE12 * newE12 + mom2.e() * mom2.e() - mom1.e() * mom1.e()) / newE12;
557  G4double Q = -(mom1[mom_ind] + mom2[mom_ind]) / newE12;
558  G4double UDQ = 1.0 / (Q * Q - 1.0);
559  G4double W = (R * Q + mom2[mom_ind]) * UDQ;
560  G4double V = (mom2.e() * mom2.e() - R * R) * UDQ;
561  G4double DET = W * W + V;
562 
563  if (DET < 0.0) {
564  if (verboseLevel > 2) G4cout << " DET < 0 " << G4endl;
565  return;
566  }
567 
568  // Tuning allowed only for non-negative determinant
569  G4double x1 = -(W + std::sqrt(DET));
570  G4double x2 = -(W - std::sqrt(DET));
571 
572  // choose the appropriate solution
573  G4bool xset = false;
574  G4double x = 0.0;
575 
576  if(mom_non_cons.e() > 0.0) { // x has to be > 0.0
577  if(x1 > 0.0) {
578  if(R + Q * x1 >= 0.0) {
579  x = x1;
580  xset = true;
581  };
582  };
583  if(!xset && x2 > 0.0) {
584  if(R + Q * x2 >= 0.0) {
585  x = x2;
586  xset = true;
587  };
588  };
589  } else {
590  if(x1 < 0.0) {
591  if(R + Q * x1 >= 0.) {
592  x = x1;
593  xset = true;
594  };
595  };
596  if(!xset && x2 < 0.0) {
597  if(R + Q * x2 >= 0.0) {
598  x = x2;
599  xset = true;
600  };
601  };
602  } // if(mom_non_cons.e() > 0.0)
603 
604  if(!xset) {
605  if(verboseLevel > 2)
606  G4cout << " no appropriate solution found " << G4endl;
607  return;
608  } // if(xset)
609 
610  // retune momentums
611  mom1[mom_ind] += x;
612  mom2[mom_ind] -= x;
613  outgoingParticles[tune_particles.first ].setMomentum(mom1);
614  outgoingParticles[tune_particles.second].setMomentum(mom2);
615  out_mom = getTotalOutputMomentum();
616  std::sort(outgoingParticles.begin(), outgoingParticles.end(), G4ParticleLargerEkin());
617  mom_non_cons = ini_mom - out_mom;
618  pnc = mom_non_cons.rho();
619  enc = mom_non_cons.e();
620 
621  on_shell = (std::fabs(enc) < accuracy || pnc < accuracy);
622 
623  if(verboseLevel > 2) {
624  G4cout << " momentum non conservation tuning: " << G4endl
625  << " e " << enc << " p " << pnc
626  << (on_shell?" success":" FAILURE") << G4endl;
627  }
628 }
629 
630 
631 // Returns excitation energy in GeV
632 
634  eex_rest = 0.;
635  G4int i(0);
636  for (i=0; i < numberOfOutgoingNuclei(); i++)
637  eex_rest += outgoingNuclei[i].getExitationEnergyInGeV();
638  for (i=0; i < numberOfFragments(); i++)
639  eex_rest += recoilFragments[i].GetExcitationEnergy() / GeV;
640 }
641 
642 
643 std::pair<std::pair<G4int, G4int>, G4int>
645  if (verboseLevel > 2)
646  G4cout << " >>> G4CollisionOutput::selectPairToTune" << G4endl;
647 
648  std::pair<G4int, G4int> tup(-1, -1);
649  G4int i3 = -1;
650  std::pair<std::pair<G4int, G4int>, G4int> tuner(tup, i3); // Set invalid
651 
652  if (outgoingParticles.size() < 2) return tuner; // Nothing to do
653 
654  G4int ibest1 = -1;
655  G4int ibest2 = -1;
656  G4double pbest = 0.0;
657  G4double pcut = 0.3 * std::sqrt(1.88 * std::fabs(de));
658  G4double p1 = 0.0;
659 
660  for (G4int i = 0; i < G4int(outgoingParticles.size())-1; i++) {
661  G4LorentzVector mom1 = outgoingParticles[i].getMomentum();
662 
663  for (G4int j = i+1; j < G4int(outgoingParticles.size()); j++) {
664  G4LorentzVector mom2 = outgoingParticles[j].getMomentum();
665 
666  for (G4int l = G4LorentzVector::X; l<=G4LorentzVector::Z; l++) {
667  if (mom1[l]*mom2[l]<0.0) {
668  if (std::fabs(mom1[l])>pcut && std::fabs(mom2[l])>pcut) {
669  G4double psum = std::fabs(mom1[l]) + std::fabs(mom2[l]);
670 
671  if(psum > pbest) {
672  ibest1 = i;
673  ibest2 = j;
674  i3 = l;
675  p1 = mom1[l];
676  pbest = psum;
677  } // psum > pbest
678  } // mom1 and mom2 > pcut
679  } // mom1 ~ -mom2
680  } // for (G4int l ...
681  } // for (G4int j ...
682  } // for (G4int i ...
683 
684  if (i3 < 0) return tuner;
685 
686  tuner.second = i3; // Momentum axis for tuning
687 
688  // NOTE: Sign of de determines order for special case of p1==0.
689  if (de > 0.0) {
690  tuner.first.first = (p1>0.) ? ibest1 : ibest2;
691  tuner.first.second = (p1>0.) ? ibest2 : ibest1;
692  } else {
693  tuner.first.first = (p1<0.) ? ibest2 : ibest1;
694  tuner.first.second = (p1<0.) ? ibest1 : ibest2;
695  }
696 
697  return tuner;
698 }
void trivialise(G4InuclParticle *bullet, G4InuclParticle *target)
G4CollisionOutput & operator=(const G4CollisionOutput &right)
std::vector< G4InuclNuclei >::iterator nucleiIterator
void addOutgoingNuclei(const std::vector< G4InuclNuclei > &nuclea)
std::vector< G4Fragment > recoilFragments
CLHEP::HepLorentzRotation G4LorentzRotation
G4LorentzVector getMomentum() const
void removeRecoilFragment(G4int index=-1)
void printCollisionOutput(std::ostream &os=G4cout) const
void addOutgoingParticle(const G4InuclElementaryParticle &particle)
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 add(const G4CollisionOutput &right)
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)
std::pair< std::pair< G4int, G4int >, G4int > selectPairToTune(G4double de) const
static const double GeV
Definition: G4SIunits.hh:196
void rotateEvent(const G4LorentzRotation &rotate)
G4int GetAtomicMass() const
std::vector< G4InuclElementaryParticle > outgoingParticles
G4int numberOfOutgoingNuclei() const
std::vector< G4InuclElementaryParticle >::iterator particleIterator
G4bool reflectionNeeded() const
void addOutgoingParticles(const std::vector< G4InuclElementaryParticle > &particles)
G4int getTotalBaryonNumber() const
static const G4Fragment emptyFragment
std::vector< G4InuclNuclei >::const_iterator nucleiIterator
G4int getTotalCharge() const
std::vector< G4InuclElementaryParticle >::iterator particleIterator
Definition: G4BigBanger.cc:64
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)
std::vector< G4InuclNuclei > outgoingNuclei
CLHEP::HepLorentzVector G4LorentzVector