Geant4 Cross Reference |
1 // 1 2 // ******************************************* 3 // * License and Disclaimer 4 // * 5 // * The Geant4 software is copyright of th 6 // * the Geant4 Collaboration. It is provided 7 // * conditions of the Geant4 Software License 8 // * LICENSE and available at http://cern.ch/ 9 // * include a list of copyright holders. 10 // * 11 // * Neither the authors of this software syst 12 // * institutes,nor the agencies providing fin 13 // * work make any representation or warran 14 // * regarding this software system or assum 15 // * use. Please see the license in the file 16 // * for the full disclaimer and the limitatio 17 // * 18 // * This code implementation is the result 19 // * technical work of the GEANT4 collaboratio 20 // * By using, copying, modifying or distri 21 // * any work based on the software) you ag 22 // * use in resulting scientific publicati 23 // * acceptance of all terms of the Geant4 Sof 24 // ******************************************* 25 // 26 // 27 // Author: Luciano Pandola 28 // 29 // History: 30 // -------- 31 // 09 Dec 2009 L Pandola First implementa 32 // 33 #include "G4PenelopeSamplingData.hh" 34 35 //....oooOO0OOooo........oooOO0OOooo........oo 36 G4PenelopeSamplingData::G4PenelopeSamplingData 37 fNP(nPoints) 38 { 39 //create vectors 40 fX = new G4DataVector(); 41 fPAC = new G4DataVector(); 42 fA = new G4DataVector(); 43 fB = new G4DataVector(); 44 fITTL = new std::vector<size_t>; 45 fITTU = new std::vector<size_t>; 46 } 47 48 //....oooOO0OOooo........oooOO0OOooo........oo 49 G4PenelopeSamplingData::~G4PenelopeSamplingDat 50 { 51 if (fX) delete fX; 52 if (fPAC) delete fPAC; 53 if (fA) delete fA; 54 if (fB) delete fB; 55 if (fITTL) delete fITTL; 56 if (fITTU) delete fITTU; 57 } 58 59 //....oooOO0OOooo........oooOO0OOooo........oo 60 61 size_t G4PenelopeSamplingData::GetNumberOfStor 62 { 63 size_t points = fX->size(); 64 65 //check everything is all right 66 if (fPAC->size() != points || fA->size() != 67 fB->size() != points || fITTL->size() != 68 fITTU->size() != points) 69 { 70 G4ExceptionDescription ed; 71 ed << "Data vectors look to have differe 72 G4Exception("G4PenelopeSamplingData::Get 73 FatalException,ed); 74 } 75 return points; 76 } 77 78 //....oooOO0OOooo........oooOO0OOooo........oo 79 void G4PenelopeSamplingData::Clear() 80 { 81 if (fX) delete fX; 82 if (fPAC) delete fPAC; 83 if (fA) delete fA; 84 if (fB) delete fB; 85 if (fITTL) delete fITTL; 86 if (fITTU) delete fITTU; 87 //create vectors 88 fX = new G4DataVector(); 89 fPAC = new G4DataVector(); 90 fA = new G4DataVector(); 91 fB = new G4DataVector(); 92 fITTL = new std::vector<size_t>; 93 fITTU = new std::vector<size_t>; 94 } 95 96 //....oooOO0OOooo........oooOO0OOooo........oo 97 void G4PenelopeSamplingData::AddPoint(G4double 98 size_t ITTL0,size_t ITTU0) 99 { 100 fX->push_back(x0); 101 fPAC->push_back(pac0); 102 fA->push_back(a0); 103 fB->push_back(b0); 104 fITTL->push_back(ITTL0); 105 fITTU->push_back(ITTU0); 106 107 //check how many points we do have now 108 size_t nOfPoints = GetNumberOfStoredPoints() 109 110 if (nOfPoints > ((size_t)fNP)) 111 { 112 G4cout << "G4PenelopeSamplingData::AddPo 113 G4cout << "WARNING: Up to now there are 114 G4cout << "while the anticipated (declar 115 } 116 return; 117 } 118 119 //....oooOO0OOooo........oooOO0OOooo........oo 120 void G4PenelopeSamplingData::DumpTable() 121 { 122 123 G4cout << "********************************* 124 G4cout << GetNumberOfStoredPoints() << " poi 125 G4cout << "********************************* 126 for (size_t i=0;i<GetNumberOfStoredPoints(); 127 { 128 G4cout << i << " " << (*fX)[i] << " " << 129 (*fB)[i] << " " << (*fITTL)[i] << " " << (*f 130 } 131 G4cout << "********************************* 132 } 133 134 //....oooOO0OOooo........oooOO0OOooo........oo 135 G4double G4PenelopeSamplingData::GetX(size_t i 136 { 137 if (index < fX->size()) 138 return (*fX)[index]; 139 else 140 return 0; 141 } 142 143 //....oooOO0OOooo........oooOO0OOooo........oo 144 G4double G4PenelopeSamplingData::GetPAC(size_t 145 { 146 if (index < fPAC->size()) 147 return (*fPAC)[index]; 148 else 149 return 0; 150 } 151 152 //....oooOO0OOooo........oooOO0OOooo........oo 153 G4double G4PenelopeSamplingData::GetA(size_t i 154 { 155 if (index < fA->size()) 156 return (*fA)[index]; 157 else 158 return 0; 159 } 160 161 //....oooOO0OOooo........oooOO0OOooo........oo 162 G4double G4PenelopeSamplingData::GetB(size_t i 163 { 164 if (index < fB->size()) 165 return (*fB)[index]; 166 else 167 return 0; 168 } 169 170 //....oooOO0OOooo........oooOO0OOooo........oo 171 G4double G4PenelopeSamplingData::SampleValue(G 172 { 173 //One passes here a random number in (0,1). 174 //Notice: it possible that is between (0,b) 175 size_t points = GetNumberOfStoredPoints(); 176 177 size_t itn = (size_t) (maxRand*(points-1)); 178 size_t i = (*fITTL)[itn]; 179 size_t j = (*fITTU)[itn]; 180 181 while ((j-i) > 1) 182 { 183 size_t k = (i+j)/2; 184 if (maxRand > (*fPAC)[k]) 185 i = k; 186 else 187 j = k; 188 } 189 190 //Sampling from the rational inverse cumulat 191 G4double result = 0; 192 193 G4double rr = maxRand - (*fPAC)[i]; 194 if (rr > 1e-16) 195 { 196 G4double d = (*fPAC)[i+1]-(*fPAC)[i]; 197 result = (*fX)[i]+ 198 ((1.0+(*fA)[i]+(*fB)[i])*d*rr/ 199 (d*d+((*fA)[i]*d+(*fB)[i]*rr)*rr))*((*fX)[i 200 } 201 else 202 result = (*fX)[i]; 203 204 return result; 205 } 206