Geant4  10.03.p01
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
G4LEPTSDistribution.cc
Go to the documentation of this file.
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 #include "G4LEPTSDistribution.hh"
27 
28 
30 {
31 }
32 
33 
35 
36  G4int eB, out, out2;
37  float float_data1,float_data2;
38  G4double sum, esum;
39  FILE * fp;
40 
41  for (eB=0;eB<10000;eB++){
42  E[eB]=0.0;
43  f[eB]=0.0;
44  F[eB]=0.0;
45  eF[eB]=0.0;
46  }
47 
48  if ((fp=fopen(fileName.c_str(), "r"))==NULL){
49  //G4cout << "Error reading " << fileName << G4endl;
50  NoBins = 0;
51  bFileFound = false;
52  return;
53  }
54  else{
55  bFileFound = true;
56  // G4cout << "Read Distro (" << fileName << ") " << G4endl;
57  out=1;
58  eB=1;
59  while (out==1){
60  out = fscanf(fp,"%f \n",&float_data1);
61  out2 = fscanf(fp,"%f \n",&float_data2);
62  if (out==1 && out2==1){
63  E[eB]=(G4double)float_data1;
64  f[eB]=(G4double)float_data2;
65  eB++;
66  }
67  }
68 
69  fclose(fp);
70  }
71 
72  NoBins=eB-1; //=1272+1 or 9607+1;
73 
74  if( NoBins >= NMAX )
75  printf("ERROR !!!! Eloss NoBins= %d \n", NoBins);
76 
77  sum=0.0;
78  esum=0.0;
79  for (eB=0;eB<=NoBins;eB++) {
80  if( f[eB] > 0) {
81  sum+=f[eB];
82  esum+=E[eB]*f[eB];
83  }
84  F[eB]=sum;
85  eF[eB]=esum;
86  }
87 
88  // if( verboseLevel >= 1 ) G4cout << "Norm: " << F[NoBins] << " NoBins: "<< NoBins << G4endl;
89 
90  for (eB=0;eB<=NoBins;eB++) {
91  eF[eB] = eF[eB]/F[eB];
92  F[eB] = F[eB]/F[NoBins];
93  }
94  //for (eB=0;eB<=NoBins;eB++)
95  //G4cout << "eff " << E[eB] << " " << f[eB] << " " << F[eB] << "\n";
96 }
97 
99 {
100 
101  G4int eB, out, out2;
102  float float_data1,float_data2;
103  G4double sum, esum;
104 
105  for (eB=0;eB<10000;eB++){
106  E[eB]=0.0;
107  f[eB]=0.0;
108  F[eB]=0.0;
109  eF[eB]=0.0;
110  }
111 
112  bFileFound = true;
113  out=1;
114  eB=1;
115 
116  for( G4int id = 0; id < nData; id++ ){
117  out = fscanf(fp,"%f \n",&float_data1);
118  out2 = fscanf(fp,"%f \n",&float_data2);
119  if (out==1 && out2==1){
120  E[eB]=(G4double)float_data1;
121  f[eB]=(G4double)float_data2;
122  eB++;
123  }else{
124  return 1;
125  }
126  }
127 
128  NoBins=eB-1; //=1272+1 or 9607+1;
129 
130  if( NoBins >= NMAX )
131  printf("ERROR !!!! Eloss NoBins= %d \n", NoBins);
132 
133  sum=0.0;
134  esum=0.0;
135  for (eB=0;eB<=NoBins;eB++) {
136  if( f[eB] > 0) {
137  sum+=f[eB];
138  esum+=E[eB]*f[eB];
139  }
140  F[eB]=sum;
141  eF[eB]=esum;
142  }
143 
144  //if( verboseLevel >= 1 ) G4cout << "Norm: " << F[NoBins] << " NoBins: "<< NoBins << G4endl;
145 
146  for (eB=0;eB<=NoBins;eB++) {
147  eF[eB] = eF[eB]/F[eB];
148  F[eB] = F[eB]/F[NoBins];
149  }
150  //for (eB=0;eB<=NoBins;eB++)
151  //G4cout << "eff " << E[eB] << " " << f[eB] << " " << F[eB] << "\n";
152 
153  return 0;
154 }
155 
156 
157 
159 // Sample Energy from Cumulative distr. G4interval [eMin, eMax]
160 
161  if( eMin > eMax) return 0.0;
162 
163  G4int i,j,k=0, iMin, iMax;
164 
165  i=0; j=NoBins;
166  while ((j-i)>1) {
167  k=(i+j)/2;
168  if( E[k] < eMax ) i=k;
169  else j=k;
170  }
171  iMax = i;
172 
173  i=0; j=NoBins;
174  while ((j-i)>1) {
175  k=(i+j)/2;
176  if( E[k] < eMin ) i=k;
177  else j=k;
178  }
179  iMin = i;
180 
181  G4double rnd = F[iMin] + (F[iMax] - F[iMin]) * G4UniformRand();
182 
183  i=0; j=NoBins;
184  while ((j-i)>1) {
185  k=(i+j)/2;
186  if( F[k]<rnd) i=k;
187  else j=k;
188  }
189 
190  G4double Sampled = E[k];
191 
192  if( Sampled < eMin) Sampled = eMin;
193  else if( Sampled > eMax) Sampled = eMax;
194 
195  return Sampled;
196 }
int G4int
Definition: G4Types.hh:78
#define G4UniformRand()
Definition: Randomize.hh:97
void ReadFile(G4String fileName)
G4double Sample(G4double, G4double)
bool G4bool
Definition: G4Types.hh:79
#define NMAX
Definition: adler32.cc:15
double G4double
Definition: G4Types.hh:76