| 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