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 // G4LEPTSDistribution << 27 // << 28 // Author: Pedro Arce (CIEMAT), 2014 << 29 // ------------------------------------------- << 30 << 31 #include "G4LEPTSDistribution.hh" 26 #include "G4LEPTSDistribution.hh" 32 27 33 #include <stdio.h> << 34 #include <iostream> << 35 28 36 void G4LEPTSDistribution::ReadFile(const G4Str << 29 G4LEPTSDistribution::G4LEPTSDistribution() 37 { 30 { >> 31 } >> 32 >> 33 >> 34 void G4LEPTSDistribution::ReadFile(G4String fileName) { >> 35 38 G4int eB, out, out2; 36 G4int eB, out, out2; 39 G4float float_data1,float_data2; << 37 float float_data1,float_data2; 40 G4double sum, esum; 38 G4double sum, esum; 41 FILE * fp; 39 FILE * fp; 42 40 43 for (eB=0; eB<10000; ++eB) << 41 for (eB=0;eB<10000;eB++){ 44 { << 45 E[eB]=0.0; 42 E[eB]=0.0; 46 f[eB]=0.0; 43 f[eB]=0.0; 47 F[eB]=0.0; 44 F[eB]=0.0; 48 eF[eB]=0.0; 45 eF[eB]=0.0; 49 } 46 } 50 47 51 if ((fp=fopen(fileName.c_str(), "r"))==nullp << 48 if ((fp=fopen(fileName.c_str(), "r"))==NULL){ 52 { << 49 //G4cout << "Error reading " << fileName << G4endl; 53 NoBins = 0; 50 NoBins = 0; 54 bFileFound = false; 51 bFileFound = false; 55 return; 52 return; 56 } 53 } 57 << 54 else{ 58 bFileFound = true; << 55 bFileFound = true; 59 << 56 // G4cout << "Read Distro (" << fileName << ") " << G4endl; 60 out=1; << 57 out=1; 61 eB=1; << 58 eB=1; 62 while (out==1) << 59 while (out==1){ 63 { << 60 out = fscanf(fp,"%f \n",&float_data1); 64 out = fscanf(fp,"%f \n",&float_data1); << 61 out2 = fscanf(fp,"%f \n",&float_data2); 65 out2 = fscanf(fp,"%f \n",&float_data2); << 62 if (out==1 && out2==1){ 66 if (out==1 && out2==1) << 63 E[eB]=(G4double)float_data1; 67 { << 64 f[eB]=(G4double)float_data2; 68 E[eB]=(G4double)float_data1; << 65 eB++; 69 f[eB]=(G4double)float_data2; << 66 } 70 ++eB; << 71 } 67 } 72 } << 73 68 74 fclose(fp); << 69 fclose(fp); >> 70 } 75 71 76 NoBins=eB-1; //=1272+1 or 9607+1; 72 NoBins=eB-1; //=1272+1 or 9607+1; 77 73 78 if( NoBins >= NMAX ) 74 if( NoBins >= NMAX ) 79 { << 75 printf("ERROR !!!! Eloss NoBins= %d \n", NoBins); 80 std::ostringstream message; << 81 message << "ERROR !!!! Eloss NoBins = " < << 82 G4Exception("G4LEPTSDistribution::ReadFile << 83 FatalException, message); << 84 } << 85 76 86 sum=0.0; 77 sum=0.0; 87 esum=0.0; 78 esum=0.0; 88 for (eB=0; eB<=NoBins; ++eB) << 79 for (eB=0;eB<=NoBins;eB++) { 89 { << 80 if( f[eB] > 0) { 90 if( f[eB] > 0) << 91 { << 92 sum+=f[eB]; 81 sum+=f[eB]; 93 esum+=E[eB]*f[eB]; 82 esum+=E[eB]*f[eB]; 94 } 83 } 95 F[eB]=sum; 84 F[eB]=sum; 96 eF[eB]=esum; 85 eF[eB]=esum; 97 } 86 } 98 87 99 for (eB=0; eB<=NoBins; ++eB) << 88 // if( verboseLevel >= 1 ) G4cout << "Norm: " << F[NoBins] << " NoBins: "<< NoBins << G4endl; 100 { << 89 >> 90 for (eB=0;eB<=NoBins;eB++) { 101 eF[eB] = eF[eB]/F[eB]; 91 eF[eB] = eF[eB]/F[eB]; 102 F[eB] = F[eB]/F[NoBins]; 92 F[eB] = F[eB]/F[NoBins]; 103 } 93 } >> 94 //for (eB=0;eB<=NoBins;eB++) >> 95 //G4cout << "eff " << E[eB] << " " << f[eB] << " " << F[eB] << "\n"; 104 } 96 } 105 97 106 G4bool G4LEPTSDistribution::ReadFile( FILE* fp 98 G4bool G4LEPTSDistribution::ReadFile( FILE* fp, G4int nData ) 107 { 99 { 108 100 109 G4int eB, out, out2; 101 G4int eB, out, out2; 110 G4float float_data1,float_data2; << 102 float float_data1,float_data2; 111 G4double sum, esum; 103 G4double sum, esum; 112 104 113 for (eB=0; eB<10000; ++eB) << 105 for (eB=0;eB<10000;eB++){ 114 { << 115 E[eB]=0.0; 106 E[eB]=0.0; 116 f[eB]=0.0; 107 f[eB]=0.0; 117 F[eB]=0.0; 108 F[eB]=0.0; 118 eF[eB]=0.0; 109 eF[eB]=0.0; 119 } 110 } 120 111 121 bFileFound = true; 112 bFileFound = true; 122 out=1; 113 out=1; 123 eB=1; 114 eB=1; 124 115 125 for( G4int id = 0; id < nData; ++id ) << 116 for( G4int id = 0; id < nData; id++ ){ 126 { << 127 out = fscanf(fp,"%f \n",&float_data1); 117 out = fscanf(fp,"%f \n",&float_data1); 128 out2 = fscanf(fp,"%f \n",&float_data2); 118 out2 = fscanf(fp,"%f \n",&float_data2); 129 if (out==1 && out2==1){ 119 if (out==1 && out2==1){ 130 E[eB]=(G4double)float_data1; 120 E[eB]=(G4double)float_data1; 131 f[eB]=(G4double)float_data2; 121 f[eB]=(G4double)float_data2; 132 ++eB; << 122 eB++; 133 } << 123 }else{ 134 else << 124 return 1; 135 { << 136 return true; << 137 } 125 } 138 } 126 } 139 127 140 NoBins=eB-1; //=1272+1 or 9607+1; 128 NoBins=eB-1; //=1272+1 or 9607+1; 141 129 142 if( NoBins >= NMAX ) 130 if( NoBins >= NMAX ) 143 { << 131 printf("ERROR !!!! Eloss NoBins= %d \n", NoBins); 144 std::ostringstream message; << 145 message << "ERROR !!!! Eloss NoBins = " < << 146 G4Exception("G4LEPTSDistribution::ReadFile << 147 FatalException, message); << 148 } << 149 132 150 sum=0.0; 133 sum=0.0; 151 esum=0.0; 134 esum=0.0; 152 for (eB=0; eB<=NoBins; ++eB) << 135 for (eB=0;eB<=NoBins;eB++) { 153 { << 136 if( f[eB] > 0) { 154 if( f[eB] > 0) << 155 { << 156 sum+=f[eB]; 137 sum+=f[eB]; 157 esum+=E[eB]*f[eB]; 138 esum+=E[eB]*f[eB]; 158 } 139 } 159 F[eB]=sum; 140 F[eB]=sum; 160 eF[eB]=esum; 141 eF[eB]=esum; 161 } 142 } >> 143 >> 144 //if( verboseLevel >= 1 ) G4cout << "Norm: " << F[NoBins] << " NoBins: "<< NoBins << G4endl; 162 145 163 for (eB=0; eB<=NoBins; ++eB) << 146 for (eB=0;eB<=NoBins;eB++) { 164 { << 165 eF[eB] = eF[eB]/F[eB]; 147 eF[eB] = eF[eB]/F[eB]; 166 F[eB] = F[eB]/F[NoBins]; 148 F[eB] = F[eB]/F[NoBins]; 167 } 149 } >> 150 //for (eB=0;eB<=NoBins;eB++) >> 151 //G4cout << "eff " << E[eB] << " " << f[eB] << " " << F[eB] << "\n"; 168 152 169 return false; << 153 return 0; 170 } 154 } 171 155 172 G4double G4LEPTSDistribution::Sample( G4double << 156 173 { << 157 174 // Sample Energy from Cumulative distr. G4in << 158 G4double G4LEPTSDistribution::Sample( G4double eMin, G4double eMax) { >> 159 // Sample Energy from Cumulative distr. G4interval [eMin, eMax] 175 160 176 if( eMin > eMax) return 0.0; 161 if( eMin > eMax) return 0.0; 177 162 178 G4int i,j,k=0, iMin, iMax; 163 G4int i,j,k=0, iMin, iMax; 179 164 180 i=0; j=NoBins; 165 i=0; j=NoBins; 181 while ((j-i)>1) << 166 while ((j-i)>1) { 182 { << 183 k=(i+j)/2; 167 k=(i+j)/2; 184 if( E[k] < eMax ) i=k; 168 if( E[k] < eMax ) i=k; 185 else j=k; 169 else j=k; 186 } 170 } 187 iMax = i; 171 iMax = i; 188 172 189 i=0; j=NoBins; 173 i=0; j=NoBins; 190 while ((j-i)>1) << 174 while ((j-i)>1) { 191 { << 192 k=(i+j)/2; 175 k=(i+j)/2; 193 if( E[k] < eMin ) i=k; 176 if( E[k] < eMin ) i=k; 194 else j=k; 177 else j=k; 195 } 178 } 196 iMin = i; 179 iMin = i; 197 180 198 G4double rnd = F[iMin] + (F[iMax] - F[iMin]) 181 G4double rnd = F[iMin] + (F[iMax] - F[iMin]) * G4UniformRand(); 199 182 200 i=0; j=NoBins; 183 i=0; j=NoBins; 201 while ((j-i)>1) << 184 while ((j-i)>1) { 202 { << 203 k=(i+j)/2; 185 k=(i+j)/2; 204 if( F[k]<rnd) i=k; 186 if( F[k]<rnd) i=k; 205 else j=k; 187 else j=k; 206 } 188 } 207 189 208 G4double Sampled = E[k]; 190 G4double Sampled = E[k]; 209 191 210 if( Sampled < eMin) Sampled = eMin; 192 if( Sampled < eMin) Sampled = eMin; 211 else if( Sampled > eMax) Sampled = eMax; 193 else if( Sampled > eMax) Sampled = eMax; 212 194 213 return Sampled; 195 return Sampled; 214 } 196 } 215 197