ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4LEPTSDistribution.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4LEPTSDistribution.cc
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 }