ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4HadLeadBias.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4HadLeadBias.cc
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  }