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 // $Id: G4PenelopeSamplingData.cc 66241 2012-12-13 18:34:42Z gunter $ 26 // 27 // 27 // Author: Luciano Pandola 28 // Author: Luciano Pandola 28 // 29 // 29 // History: 30 // History: 30 // -------- 31 // -------- 31 // 09 Dec 2009 L Pandola First implementa 32 // 09 Dec 2009 L Pandola First implementation 32 // 33 // 33 #include "G4PenelopeSamplingData.hh" 34 #include "G4PenelopeSamplingData.hh" 34 35 35 //....oooOO0OOooo........oooOO0OOooo........oo 36 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo... 36 G4PenelopeSamplingData::G4PenelopeSamplingData 37 G4PenelopeSamplingData::G4PenelopeSamplingData(G4int nPoints) : 37 fNP(nPoints) << 38 np(nPoints) 38 { 39 { 39 //create vectors 40 //create vectors 40 fX = new G4DataVector(); << 41 x = new G4DataVector(); 41 fPAC = new G4DataVector(); << 42 pac = new G4DataVector(); 42 fA = new G4DataVector(); << 43 a = new G4DataVector(); 43 fB = new G4DataVector(); << 44 b = new G4DataVector(); 44 fITTL = new std::vector<size_t>; << 45 ITTL = new std::vector<size_t>; 45 fITTU = new std::vector<size_t>; << 46 ITTU = new std::vector<size_t>; 46 } 47 } 47 48 48 //....oooOO0OOooo........oooOO0OOooo........oo 49 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo... 49 G4PenelopeSamplingData::~G4PenelopeSamplingDat 50 G4PenelopeSamplingData::~G4PenelopeSamplingData() 50 { 51 { 51 if (fX) delete fX; << 52 if (x) delete x; 52 if (fPAC) delete fPAC; << 53 if (pac) delete pac; 53 if (fA) delete fA; << 54 if (a) delete a; 54 if (fB) delete fB; << 55 if (b) delete b; 55 if (fITTL) delete fITTL; << 56 if (ITTL) delete ITTL; 56 if (fITTU) delete fITTU; << 57 if (ITTU) delete ITTU; 57 } 58 } 58 59 59 //....oooOO0OOooo........oooOO0OOooo........oo 60 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.......oooOO0OOooo... 60 << 61 size_t G4PenelopeSamplingData::GetNumberOfStor 61 size_t G4PenelopeSamplingData::GetNumberOfStoredPoints() 62 { 62 { 63 size_t points = fX->size(); << 63 size_t points = x->size(); 64 64 65 //check everything is all right 65 //check everything is all right 66 if (fPAC->size() != points || fA->size() != << 66 if (pac->size() != points || a->size() != points || 67 fB->size() != points || fITTL->size() != << 67 b->size() != points || ITTL->size() != points || 68 fITTU->size() != points) << 68 ITTU->size() != points) 69 { 69 { 70 G4ExceptionDescription ed; 70 G4ExceptionDescription ed; 71 ed << "Data vectors look to have differe 71 ed << "Data vectors look to have different dimensions !" << G4endl; 72 G4Exception("G4PenelopeSamplingData::Get 72 G4Exception("G4PenelopeSamplingData::GetNumberOfStoredPoints()","em2040", 73 FatalException,ed); 73 FatalException,ed); 74 } 74 } 75 return points; 75 return points; 76 } 76 } 77 77 78 //....oooOO0OOooo........oooOO0OOooo........oo 78 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo... 79 void G4PenelopeSamplingData::Clear() 79 void G4PenelopeSamplingData::Clear() 80 { 80 { 81 if (fX) delete fX; << 81 if (x) delete x; 82 if (fPAC) delete fPAC; << 82 if (pac) delete pac; 83 if (fA) delete fA; << 83 if (a) delete a; 84 if (fB) delete fB; << 84 if (b) delete b; 85 if (fITTL) delete fITTL; << 85 if (ITTL) delete ITTL; 86 if (fITTU) delete fITTU; << 86 if (ITTU) delete ITTU; 87 //create vectors 87 //create vectors 88 fX = new G4DataVector(); << 88 x = new G4DataVector(); 89 fPAC = new G4DataVector(); << 89 pac = new G4DataVector(); 90 fA = new G4DataVector(); << 90 a = new G4DataVector(); 91 fB = new G4DataVector(); << 91 b = new G4DataVector(); 92 fITTL = new std::vector<size_t>; << 92 ITTL = new std::vector<size_t>; 93 fITTU = new std::vector<size_t>; << 93 ITTU = new std::vector<size_t>; 94 } 94 } 95 95 96 //....oooOO0OOooo........oooOO0OOooo........oo 96 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo... 97 void G4PenelopeSamplingData::AddPoint(G4double 97 void G4PenelopeSamplingData::AddPoint(G4double x0,G4double pac0,G4double a0,G4double b0, 98 size_t ITTL0,size_t ITTU0) 98 size_t ITTL0,size_t ITTU0) 99 { 99 { 100 fX->push_back(x0); << 100 x->push_back(x0); 101 fPAC->push_back(pac0); << 101 pac->push_back(pac0); 102 fA->push_back(a0); << 102 a->push_back(a0); 103 fB->push_back(b0); << 103 b->push_back(b0); 104 fITTL->push_back(ITTL0); << 104 ITTL->push_back(ITTL0); 105 fITTU->push_back(ITTU0); << 105 ITTU->push_back(ITTU0); 106 106 107 //check how many points we do have now 107 //check how many points we do have now 108 size_t nOfPoints = GetNumberOfStoredPoints() 108 size_t nOfPoints = GetNumberOfStoredPoints(); 109 109 110 if (nOfPoints > ((size_t)fNP)) << 110 if (nOfPoints > ((size_t)np)) 111 { 111 { 112 G4cout << "G4PenelopeSamplingData::AddPo 112 G4cout << "G4PenelopeSamplingData::AddPoint() " << G4endl; 113 G4cout << "WARNING: Up to now there are 113 G4cout << "WARNING: Up to now there are " << nOfPoints << " points in the table" << G4endl; 114 G4cout << "while the anticipated (declar << 114 G4cout << "while the anticipated (declared) number is " << np << G4endl; 115 } 115 } 116 return; 116 return; 117 } 117 } 118 118 119 //....oooOO0OOooo........oooOO0OOooo........oo 119 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.. 120 void G4PenelopeSamplingData::DumpTable() 120 void G4PenelopeSamplingData::DumpTable() 121 { 121 { 122 122 123 G4cout << "********************************* 123 G4cout << "*************************************************************************" << G4endl; 124 G4cout << GetNumberOfStoredPoints() << " poi 124 G4cout << GetNumberOfStoredPoints() << " points" << G4endl; 125 G4cout << "********************************* 125 G4cout << "*************************************************************************" << G4endl; 126 for (size_t i=0;i<GetNumberOfStoredPoints(); 126 for (size_t i=0;i<GetNumberOfStoredPoints();i++) 127 { 127 { 128 G4cout << i << " " << (*fX)[i] << " " << << 128 G4cout << i << " " << (*x)[i] << " " << (*pac)[i] << " " << (*a)[i] << " " << 129 (*fB)[i] << " " << (*fITTL)[i] << " " << (*f << 129 (*b)[i] << " " << (*ITTL)[i] << " " << (*ITTU)[i] << G4endl; 130 } 130 } 131 G4cout << "********************************* 131 G4cout << "*************************************************************************" << G4endl; 132 } 132 } 133 133 134 //....oooOO0OOooo........oooOO0OOooo........oo 134 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.. 135 G4double G4PenelopeSamplingData::GetX(size_t i 135 G4double G4PenelopeSamplingData::GetX(size_t index) 136 { 136 { 137 if (index < fX->size()) << 137 if (index < x->size()) 138 return (*fX)[index]; << 138 return (*x)[index]; 139 else 139 else 140 return 0; 140 return 0; 141 } 141 } 142 142 143 //....oooOO0OOooo........oooOO0OOooo........oo 143 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.. 144 G4double G4PenelopeSamplingData::GetPAC(size_t 144 G4double G4PenelopeSamplingData::GetPAC(size_t index) 145 { 145 { 146 if (index < fPAC->size()) << 146 if (index < pac->size()) 147 return (*fPAC)[index]; << 147 return (*pac)[index]; 148 else 148 else 149 return 0; 149 return 0; 150 } 150 } 151 151 152 //....oooOO0OOooo........oooOO0OOooo........oo 152 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.. 153 G4double G4PenelopeSamplingData::GetA(size_t i 153 G4double G4PenelopeSamplingData::GetA(size_t index) 154 { 154 { 155 if (index < fA->size()) << 155 if (index < a->size()) 156 return (*fA)[index]; << 156 return (*a)[index]; 157 else 157 else 158 return 0; 158 return 0; 159 } 159 } 160 160 161 //....oooOO0OOooo........oooOO0OOooo........oo 161 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.. 162 G4double G4PenelopeSamplingData::GetB(size_t i 162 G4double G4PenelopeSamplingData::GetB(size_t index) 163 { 163 { 164 if (index < fB->size()) << 164 if (index < b->size()) 165 return (*fB)[index]; << 165 return (*b)[index]; 166 else 166 else 167 return 0; 167 return 0; 168 } 168 } 169 169 170 //....oooOO0OOooo........oooOO0OOooo........oo 170 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.. 171 G4double G4PenelopeSamplingData::SampleValue(G 171 G4double G4PenelopeSamplingData::SampleValue(G4double maxRand) 172 { 172 { 173 //One passes here a random number in (0,1). 173 //One passes here a random number in (0,1). 174 //Notice: it possible that is between (0,b) 174 //Notice: it possible that is between (0,b) with b<1 175 size_t points = GetNumberOfStoredPoints(); 175 size_t points = GetNumberOfStoredPoints(); 176 176 177 size_t itn = (size_t) (maxRand*(points-1)); 177 size_t itn = (size_t) (maxRand*(points-1)); 178 size_t i = (*fITTL)[itn]; << 178 size_t i = (*ITTL)[itn]; 179 size_t j = (*fITTU)[itn]; << 179 size_t j = (*ITTU)[itn]; 180 180 181 while ((j-i) > 1) 181 while ((j-i) > 1) 182 { 182 { 183 size_t k = (i+j)/2; 183 size_t k = (i+j)/2; 184 if (maxRand > (*fPAC)[k]) << 184 if (maxRand > (*pac)[k]) 185 i = k; 185 i = k; 186 else 186 else 187 j = k; 187 j = k; 188 } 188 } 189 189 190 //Sampling from the rational inverse cumulat 190 //Sampling from the rational inverse cumulative distribution 191 G4double result = 0; 191 G4double result = 0; 192 192 193 G4double rr = maxRand - (*fPAC)[i]; << 193 G4double rr = maxRand - (*pac)[i]; 194 if (rr > 1e-16) 194 if (rr > 1e-16) 195 { 195 { 196 G4double d = (*fPAC)[i+1]-(*fPAC)[i]; << 196 G4double d = (*pac)[i+1]-(*pac)[i]; 197 result = (*fX)[i]+ << 197 result = (*x)[i]+ 198 ((1.0+(*fA)[i]+(*fB)[i])*d*rr/ << 198 ((1.0+(*a)[i]+(*b)[i])*d*rr/ 199 (d*d+((*fA)[i]*d+(*fB)[i]*rr)*rr))*((*fX)[i << 199 (d*d+((*a)[i]*d+(*b)[i]*rr)*rr))*((*x)[i+1]-(*x)[i]); 200 } 200 } 201 else 201 else 202 result = (*fX)[i]; << 202 result = (*x)[i]; 203 203 204 return result; 204 return result; 205 } 205 } 206 206