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