Geant4 Cross Reference |
1 // 1 // 2 // ******************************************* 2 // ******************************************************************** 3 // * License and Disclaimer 3 // * License and Disclaimer * 4 // * 4 // * * 5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. 9 // * include a list of copyright holders. * 10 // * 10 // * * 11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. * 17 // * 17 // * * 18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. * 24 // ******************************************* 24 // ******************************************************************** 25 // 25 // 26 // 20110906 M. Kelsey -- Use reference to G4H 26 // 20110906 M. Kelsey -- Use reference to G4HadSecondary instead of pointer 27 27 28 #include "G4HadLeadBias.hh" 28 #include "G4HadLeadBias.hh" 29 #include "G4Gamma.hh" 29 #include "G4Gamma.hh" 30 #include "G4PionZero.hh" 30 #include "G4PionZero.hh" 31 #include "Randomize.hh" 31 #include "Randomize.hh" 32 #include "G4HadFinalState.hh" 32 #include "G4HadFinalState.hh" 33 33 34 G4HadFinalState * G4HadLeadBias::Bias(G4HadF 34 G4HadFinalState * G4HadLeadBias::Bias(G4HadFinalState * result) 35 { 35 { 36 // G4cerr << "bias enter"<<G4endl; 36 // G4cerr << "bias enter"<<G4endl; 37 G4int nMeson(0), nBaryon(0), npi0(0), ngam 37 G4int nMeson(0), nBaryon(0), npi0(0), ngamma(0), nLepton(0); 38 G4int i(0); 38 G4int i(0); 39 G4int maxE = -1; 39 G4int maxE = -1; 40 G4double emax = 0; 40 G4double emax = 0; 41 if(result->GetStatusChange()==isAlive) 41 if(result->GetStatusChange()==isAlive) 42 { 42 { 43 emax = result->GetEnergyChange(); 43 emax = result->GetEnergyChange(); 44 } 44 } 45 //G4cout << "max energy "<<G4endl; 45 //G4cout << "max energy "<<G4endl; 46 for(i=0;i<static_cast<G4int>(result->GetNu << 46 for(i=0;i<result->GetNumberOfSecondaries();i++) 47 { 47 { 48 if(result->GetSecondary(i)->GetParticle( 48 if(result->GetSecondary(i)->GetParticle()->GetKineticEnergy()>emax) 49 { 49 { 50 maxE = i; 50 maxE = i; 51 emax = result->GetSecondary(i)->GetParticle( 51 emax = result->GetSecondary(i)->GetParticle()->GetKineticEnergy(); 52 } 52 } 53 } 53 } 54 //G4cout <<"loop1"<<G4endl; 54 //G4cout <<"loop1"<<G4endl; 55 for(i=0; i<static_cast<G4int>(result->GetN << 55 for(i=0; i<result->GetNumberOfSecondaries(); i++) 56 { 56 { 57 const G4DynamicParticle* aSecTrack = res 57 const G4DynamicParticle* aSecTrack = result->GetSecondary(i)->GetParticle(); 58 if(i==maxE) 58 if(i==maxE) 59 { 59 { 60 } 60 } 61 else if(aSecTrack->GetDefinition()->GetB 61 else if(aSecTrack->GetDefinition()->GetBaryonNumber()!=0) 62 { 62 { 63 nBaryon++; 63 nBaryon++; 64 } 64 } 65 else if(aSecTrack->GetDefinition()->GetL 65 else if(aSecTrack->GetDefinition()->GetLeptonNumber()!=0) 66 { 66 { 67 nLepton++; 67 nLepton++; 68 } 68 } 69 else if(aSecTrack->GetDefinition()==G4Ga 69 else if(aSecTrack->GetDefinition()==G4Gamma::Gamma()) 70 { 70 { 71 ngamma++; 71 ngamma++; 72 } 72 } 73 else if(aSecTrack->GetDefinition()==G4Pi 73 else if(aSecTrack->GetDefinition()==G4PionZero::PionZero()) 74 { 74 { 75 npi0++; 75 npi0++; 76 } 76 } 77 else 77 else 78 { 78 { 79 nMeson++; 79 nMeson++; 80 } 80 } 81 } 81 } 82 //G4cout << "BiasDebug 1 = "<<result->Get 82 //G4cout << "BiasDebug 1 = "<<result->GetNumberOfSecondaries()<<" " 83 // <<nMeson<<" "<< nBaryon<<" "<< n 83 // <<nMeson<<" "<< nBaryon<<" "<< npi0<<" "<< ngamma<<" "<< nLepton<<G4endl; 84 G4double mesonWeight = nMeson; 84 G4double mesonWeight = nMeson; 85 G4double baryonWeight = nBaryon; 85 G4double baryonWeight = nBaryon; 86 G4double gammaWeight = ngamma; 86 G4double gammaWeight = ngamma; 87 G4double npi0Weight = npi0; 87 G4double npi0Weight = npi0; 88 G4double leptonWeight = nLepton; 88 G4double leptonWeight = nLepton; 89 G4int randomMeson = static_cast<G4int>((nM 89 G4int randomMeson = static_cast<G4int>((nMeson+1)*G4UniformRand()); 90 G4int randomBaryon = static_cast<G4int>((n 90 G4int randomBaryon = static_cast<G4int>((nBaryon+1)*G4UniformRand()); 91 G4int randomGamma = static_cast<G4int>((ng 91 G4int randomGamma = static_cast<G4int>((ngamma+1)*G4UniformRand()); 92 G4int randomPi0 = static_cast<G4int>((npi0 92 G4int randomPi0 = static_cast<G4int>((npi0+1)*G4UniformRand()); 93 G4int randomLepton = static_cast<G4int>((n 93 G4int randomLepton = static_cast<G4int>((nLepton+1)*G4UniformRand()); 94 94 95 std::vector<G4HadSecondary> buffer; 95 std::vector<G4HadSecondary> buffer; 96 G4int cMeson(0), cBaryon(0), cpi0(0), cgam 96 G4int cMeson(0), cBaryon(0), cpi0(0), cgamma(0), cLepton(0); 97 for(i=0; i<static_cast<G4int>(result->GetN << 97 for(i=0; i<result->GetNumberOfSecondaries(); i++) 98 { 98 { 99 G4bool aCatch = false; 99 G4bool aCatch = false; 100 G4double weight = 1; 100 G4double weight = 1; 101 const G4HadSecondary* aSecTrack = result 101 const G4HadSecondary* aSecTrack = result->GetSecondary(i); 102 G4ParticleDefinition* aSecDef = aSecTrac 102 G4ParticleDefinition* aSecDef = aSecTrack->GetParticle()->GetDefinition(); 103 if(i==maxE) 103 if(i==maxE) 104 { 104 { 105 aCatch = true; 105 aCatch = true; 106 weight = 1; 106 weight = 1; 107 } 107 } 108 else if(aSecDef->GetBaryonNumber()!=0) 108 else if(aSecDef->GetBaryonNumber()!=0) 109 { 109 { 110 if(++cBaryon==randomBaryon) 110 if(++cBaryon==randomBaryon) 111 { 111 { 112 aCatch = true; 112 aCatch = true; 113 weight = baryonWeight; 113 weight = baryonWeight; 114 } 114 } 115 } 115 } 116 else if(aSecDef->GetLeptonNumber()!=0) 116 else if(aSecDef->GetLeptonNumber()!=0) 117 { 117 { 118 if(++cLepton==randomLepton) 118 if(++cLepton==randomLepton) 119 { 119 { 120 aCatch = true; 120 aCatch = true; 121 weight = leptonWeight; 121 weight = leptonWeight; 122 } 122 } 123 } 123 } 124 else if(aSecDef==G4Gamma::Gamma()) 124 else if(aSecDef==G4Gamma::Gamma()) 125 { 125 { 126 if(++cgamma==randomGamma) 126 if(++cgamma==randomGamma) 127 { 127 { 128 aCatch = true; 128 aCatch = true; 129 weight = gammaWeight; 129 weight = gammaWeight; 130 } 130 } 131 } 131 } 132 else if(aSecDef==G4PionZero::PionZero()) 132 else if(aSecDef==G4PionZero::PionZero()) 133 { 133 { 134 if(++cpi0==randomPi0) 134 if(++cpi0==randomPi0) 135 { 135 { 136 aCatch = true; 136 aCatch = true; 137 weight = npi0Weight; 137 weight = npi0Weight; 138 } 138 } 139 } 139 } 140 else 140 else 141 { 141 { 142 if(++cMeson==randomMeson) 142 if(++cMeson==randomMeson) 143 { 143 { 144 aCatch = true; 144 aCatch = true; 145 weight = mesonWeight; 145 weight = mesonWeight; 146 } 146 } 147 } 147 } 148 if(aCatch) 148 if(aCatch) 149 { 149 { 150 buffer.push_back(*aSecTrack); 150 buffer.push_back(*aSecTrack); 151 buffer.back().SetWeight(aSecTrack->GetWeight 151 buffer.back().SetWeight(aSecTrack->GetWeight()*weight); 152 } 152 } 153 else 153 else 154 { 154 { 155 delete aSecTrack; 155 delete aSecTrack; 156 } 156 } 157 } 157 } 158 result->ClearSecondaries(); 158 result->ClearSecondaries(); 159 // G4cerr << "pre"<<G4endl; 159 // G4cerr << "pre"<<G4endl; 160 result->AddSecondaries(buffer); 160 result->AddSecondaries(buffer); 161 // G4cerr << "bias exit"<<G4endl; 161 // G4cerr << "bias exit"<<G4endl; 162 162 163 return result; 163 return result; 164 } 164 } 165 165