Geant4_10
G4HadLeadBias.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 // 20110906 M. Kelsey -- Use reference to G4HadSecondary instead of pointer
27 
28 #include "G4HadLeadBias.hh"
29 #include "G4Gamma.hh"
30 #include "G4PionZero.hh"
31 #include "Randomize.hh"
32 #include "G4HadFinalState.hh"
33 
35  {
36  // G4cerr << "bias enter"<<G4endl;
37  G4int nMeson(0), nBaryon(0), npi0(0), ngamma(0), nLepton(0);
38  G4int i(0);
39  G4int maxE = -1;
40  G4double emax = 0;
41  if(result->GetStatusChange()==isAlive)
42  {
43  emax = result->GetEnergyChange();
44  }
45  //G4cout << "max energy "<<G4endl;
46  for(i=0;i<result->GetNumberOfSecondaries();i++)
47  {
48  if(result->GetSecondary(i)->GetParticle()->GetKineticEnergy()>emax)
49  {
50  maxE = i;
51  emax = result->GetSecondary(i)->GetParticle()->GetKineticEnergy();
52  }
53  }
54  //G4cout <<"loop1"<<G4endl;
55  for(i=0; i<result->GetNumberOfSecondaries(); i++)
56  {
57  const G4DynamicParticle* aSecTrack = result->GetSecondary(i)->GetParticle();
58  if(i==maxE)
59  {
60  }
61  else if(aSecTrack->GetDefinition()->GetBaryonNumber()!=0)
62  {
63  nBaryon++;
64  }
65  else if(aSecTrack->GetDefinition()->GetLeptonNumber()!=0)
66  {
67  nLepton++;
68  }
69  else if(aSecTrack->GetDefinition()==G4Gamma::Gamma())
70  {
71  ngamma++;
72  }
73  else if(aSecTrack->GetDefinition()==G4PionZero::PionZero())
74  {
75  npi0++;
76  }
77  else
78  {
79  nMeson++;
80  }
81  }
82  //G4cout << "BiasDebug 1 = "<<result->GetNumberOfSecondaries()<<" "
83  // <<nMeson<<" "<< nBaryon<<" "<< npi0<<" "<< ngamma<<" "<< nLepton<<G4endl;
84  G4double mesonWeight = nMeson;
85  G4double baryonWeight = nBaryon;
86  G4double gammaWeight = ngamma;
87  G4double npi0Weight = npi0;
88  G4double leptonWeight = nLepton;
89  G4int randomMeson = static_cast<G4int>((nMeson+1)*G4UniformRand());
90  G4int randomBaryon = static_cast<G4int>((nBaryon+1)*G4UniformRand());
91  G4int randomGamma = static_cast<G4int>((ngamma+1)*G4UniformRand());
92  G4int randomPi0 = static_cast<G4int>((npi0+1)*G4UniformRand());
93  G4int randomLepton = static_cast<G4int>((nLepton+1)*G4UniformRand());
94 
95  std::vector<G4HadSecondary> buffer;
96  G4int cMeson(0), cBaryon(0), cpi0(0), cgamma(0), cLepton(0);
97  for(i=0; i<result->GetNumberOfSecondaries(); i++)
98  {
99  G4bool aCatch = false;
100  G4double weight = 1;
101  const G4HadSecondary* aSecTrack = result->GetSecondary(i);
102  G4ParticleDefinition* aSecDef = aSecTrack->GetParticle()->GetDefinition();
103  if(i==maxE)
104  {
105  aCatch = true;
106  weight = 1;
107  }
108  else if(aSecDef->GetBaryonNumber()!=0)
109  {
110  if(++cBaryon==randomBaryon)
111  {
112  aCatch = true;
113  weight = baryonWeight;
114  }
115  }
116  else if(aSecDef->GetLeptonNumber()!=0)
117  {
118  if(++cLepton==randomLepton)
119  {
120  aCatch = true;
121  weight = leptonWeight;
122  }
123  }
124  else if(aSecDef==G4Gamma::Gamma())
125  {
126  if(++cgamma==randomGamma)
127  {
128  aCatch = true;
129  weight = gammaWeight;
130  }
131  }
132  else if(aSecDef==G4PionZero::PionZero())
133  {
134  if(++cpi0==randomPi0)
135  {
136  aCatch = true;
137  weight = npi0Weight;
138  }
139  }
140  else
141  {
142  if(++cMeson==randomMeson)
143  {
144  aCatch = true;
145  weight = mesonWeight;
146  }
147  }
148  if(aCatch)
149  {
150  buffer.push_back(*aSecTrack);
151  buffer.back().SetWeight(aSecTrack->GetWeight()*weight);
152  }
153  else
154  {
155  delete aSecTrack;
156  }
157  }
158  result->ClearSecondaries();
159  // G4cerr << "pre"<<G4endl;
160  result->AddSecondaries(buffer);
161  // G4cerr << "bias exit"<<G4endl;
162 
163  return result;
164  }
void AddSecondaries(const std::vector< G4HadSecondary > &addSecs)
G4HadSecondary * GetSecondary(size_t i)
G4double GetKineticEnergy() const
G4double G4NeutronHPJENDLHEData::G4double result
#define buffer
Definition: xmlparse.cc:611
G4double GetEnergyChange() const
G4ParticleDefinition * GetDefinition() const
double weight
Definition: plottest35.C:25
int G4int
Definition: G4Types.hh:78
virtual G4HadFinalState * Bias(G4HadFinalState *result)
#define G4UniformRand()
Definition: Randomize.hh:87
bool G4bool
Definition: G4Types.hh:79
static G4Gamma * Gamma()
Definition: G4Gamma.cc:86
static G4PionZero * PionZero()
Definition: G4PionZero.cc:104
G4DynamicParticle * GetParticle()
double G4double
Definition: G4Types.hh:76
G4int GetNumberOfSecondaries() const
G4HadFinalStateStatus GetStatusChange() const
G4double GetWeight() const