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