Geant4  10.01.p01
G4CascadeHistory.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: G4CascadeHistory.cc 67796 2013-03-08 06:18:39Z mkelsey $
27 //
28 // G4CascadeHistory: Container to record all particles produced during
29 // cascade, with daughters; printout is formatted hierarchically.
30 
31 #include "G4CascadeHistory.hh"
32 #include "G4CascadParticle.hh"
34 #include "G4InuclParticle.hh"
35 #include <iostream>
36 #include <iomanip>
37 #include <ios>
38 
39 
40 // Constructors for individual entries (vertices)
41 
43  for (G4int i=0; i<10; i++) dId[i] = -1;
44  n = 0;
45 }
46 
47 
48 // Initialize buffers for new event
49 
51  if (verboseLevel>1) G4cout << " >>> G4CascadeHistory::Clear" << G4endl;
52  theHistory.clear();
53  entryPrinted.clear();
54 }
55 
56 
57 // Record a new cascade vertex (particle and daughters)
58 
60  std::vector<G4CascadParticle>& daug) {
61  if (verboseLevel>1) G4cout << " >>> G4CascadeHistory::AddVertex" << G4endl;
62 
63  // Create new entry for vertex or update particle kinematics
64  G4int id = AddEntry(cpart);
65  FillDaughters(id, daug);
66 
67  if (verboseLevel>3) {
68  G4cout << " entry " << id << " " << &theHistory[id] << " got "
69  << theHistory[id].n << " daughters:";
70  for (G4int i=0; i<theHistory[id].n; i++) {
71  G4cout << " " << theHistory[id].dId[i];
72  }
73  G4cout << G4endl;
74  }
75 
76  return id;
77 }
78 
79 void
80 G4CascadeHistory::FillDaughters(G4int iEntry, std::vector<G4CascadParticle>& daug) {
81  G4int nDaug = (G4int)daug.size();
82 
83  if (verboseLevel>1)
84  G4cout << " >>> G4CascadeHistory::FillDaughters " << iEntry << G4endl;
85 
86  // NOTE: Cannot use reference to element, as push_back can invalidate refs!
87  theHistory[iEntry].clear();
88 
89  theHistory[iEntry].n = nDaug;
90  for (G4int i=0; i<nDaug; i++) {
91  G4int id = AddEntry(daug[i]);
92  theHistory[iEntry].dId[i] = id;
93  }
94 
95  if (verboseLevel>3) {
96  G4cout << " got " << theHistory[iEntry].n << " daughters:";
97  for (G4int i=0; i<theHistory[iEntry].n; i++) {
98  G4cout << " " << theHistory[iEntry].dId[i];
99  }
100  G4cout << G4endl;
101  }
102 }
103 
104 // Add particle to history list, or update if already recorded
105 
107  AssignHistoryID(cpart); // Make sure particle has index
108 
109  G4int id = cpart.getHistoryId();
110  if (id < size()) {
111  if (verboseLevel>2)
112  G4cout << " AddEntry updating " << id << " " << &theHistory[id] << G4endl;
113  theHistory[id].cpart = cpart; // Copies kinematics
114  } else {
115  theHistory.push_back(HistoryEntry(cpart));
116  if (verboseLevel>2)
117  G4cout << " AddEntry creating " << id << " " << &theHistory.back() << G4endl;
118  }
119 
120  if (verboseLevel>3) G4cout << theHistory[id].cpart << G4endl; // Sanity check
121 
122  return id;
123 }
124 
125 // Discard particle reabsorbed during cascade
126 
128  if (verboseLevel>1) G4cout << " >>> G4CascadeHistory::DropEntry" << G4endl;
129 
130 
131  G4int id = cpart.getHistoryId(); // Particle must appear in history
132  if (id>=0) theHistory[id].n = -1; // Special flag for absorbed particle
133 }
134 
135 // Check if particle already in history, assign ID if not there
136 
138  if (cpart.getHistoryId() >= 0) return; // ID already assigned
139 
140  if (verboseLevel>2) {
141  G4cout << " >>> G4CascadeHistory::NewHistoryID assigning ID "
142  << size() << G4endl;
143  }
144 
145  cpart.setHistoryId(size());
146 }
147 
148 
149 // Generate hierarchical (indented) report of history
150 
151 std::ostream& operator<<(std::ostream& os, const G4CascadeHistory& history) {
152  history.Print(os);
153  return os;
154 }
155 
156 void G4CascadeHistory::Print(std::ostream& os) const {
157  if (verboseLevel) os << " >>> G4CascadeHistory::Print" << std::endl;
158 
159  os << " Cascade structure: vertices, (-O-) exciton, (***) outgoing"
160  << std::endl;
161 
162  for (G4int i=0; i<size(); i++) {
163  if (!PrintingDone(i)) PrintEntry(os, i);
164  }
165 }
166 
167 // Add single-line report for particle, along with daughters
168 
169 void G4CascadeHistory::PrintEntry(std::ostream& os, G4int iEntry) const {
170  if (iEntry >= size()) return; // Skip nonexistent entry
171  if (PrintingDone(iEntry)) return; // Skip entry already reported
172 
173  entryPrinted.insert(iEntry);
174 
175  const HistoryEntry& entry = theHistory[iEntry]; // For convenience
176  const G4CascadParticle& cpart = entry.cpart;
177 
178  G4int indent = cpart.getGeneration()*2;
179 
180  // Index and indentation of cascade vertex
181  std::ios::fmtflags osFlags = os.flags();
182  os.setf(std::ios::left); // Pushes all blanks to right end of output
183  os << "#" << std::setw(3+indent) << iEntry;
184  os.flags(osFlags);
185 
186  os << cpart.getParticle().getDefinition()->GetParticleName()
187  << " p " << cpart.getMomentum() << " (cosTh "
188  << cpart.getMomentum().vect().unit().z() << ")"
189  << " @ " << cpart.getPosition()
190  << " zone " << cpart.getCurrentZone();
191 
192  // Flag as final-state particle or report daughters iteratively
193  os << " (" << GuessTarget(entry) << ")";
194  if (entry.n > 0) {
195  os << " -> N=" << entry.n << std::endl;
196  for (G4int i=0; i<entry.n; i++) {
197  PrintEntry(os, entry.dId[i]);
198  }
199  } else os << std::endl;
200 }
201 
202 // Derive target of cascade step from particle and daughters
203 
204 const char* G4CascadeHistory::
206  if (verboseLevel>2) G4cout << " >>> G4CascadeHistory::GuessTarget" << G4endl;
207 
208  if (entry.n < 0) return "-O-"; // Exciton or trapped-decay
209  if (entry.n == 0) return "***"; // Outgoing (final state) particle
210 
211  const G4CascadParticle& cpart = entry.cpart; // For convenience
212  if (verboseLevel>3) G4cout << "cpart: " << cpart;
213 
214  // Compute baryon number and charge from daughters minus projectile
215  G4int targetB = -cpart.getParticle().baryon();
216  G4int targetQ = (G4int)-cpart.getParticle().getCharge();
217 
218  for (G4int i=0; i<entry.n; i++) {
219  const G4CascadParticle& cdaug = theHistory[entry.dId[i]].cpart;
220  if (verboseLevel>3)
221  G4cout << "cdaug " << i << " ID " << entry.dId[i] << ": " << cdaug;
222 
223  targetB += cdaug.getParticle().baryon();
224  targetQ += (G4int)cdaug.getParticle().getCharge();
225  }
226 
227  // Target possibilities are proton, neutron or dibaryon (pp, nn, pn)
228  if (targetB==1 && targetQ==0) return "n";
229  if (targetB==1 && targetQ==1) return "p";
230  if (targetB==2 && targetQ==0) return "nn";
231  if (targetB==2 && targetQ==1) return "pn";
232  if (targetB==2 && targetQ==2) return "pp";
233 
234  if (verboseLevel>2) {
235  G4cout << " ERROR identifying target: deltaB " << targetB
236  << " deltaQ " << targetQ << " from\n" << cpart << " to" << G4endl;
237  for (G4int j=0; j<entry.n; j++) {
238  G4cout << theHistory[entry.dId[j]].cpart;
239  }
240  }
241 
242  return "BAD TARGET"; // Should not get here if EPCollider worked right
243 }
G4LorentzVector getMomentum() const
G4bool PrintingDone(G4int iEntry) const
G4int dId[10]
G4int getGeneration() const
void setHistoryId(G4int id)
void Print(std::ostream &os) const
const G4ParticleDefinition * getDefinition() const
void AssignHistoryID(G4CascadParticle &cpart)
int G4int
Definition: G4Types.hh:78
const G4InuclElementaryParticle & getParticle() const
const G4String & GetParticleName() const
void PrintEntry(std::ostream &os, G4int iEntry) const
G4GLOB_DLL std::ostream G4cout
std::vector< HistoryEntry > theHistory
void DropEntry(const G4CascadParticle &cpart)
G4int AddEntry(G4CascadParticle &cpart)
G4int n
std::ostream & operator<<(std::ostream &os, const G4CascadeHistory &history)
G4int getCurrentZone() const
G4double getCharge() const
#define G4endl
Definition: G4ios.hh:61
G4int getHistoryId() const
const G4ThreeVector & getPosition() const
G4int size() const
void clear()
const char * GuessTarget(const HistoryEntry &entry) const
G4int AddVertex(G4CascadParticle &cpart, std::vector< G4CascadParticle > &daug)
std::set< G4int > entryPrinted
G4CascadParticle cpart
void FillDaughters(G4int iEntry, std::vector< G4CascadParticle > &daug)