source: trunk/source/processes/hadronic/models/high_energy/src/G4HEPlot.cc@ 1350

Last change on this file since 1350 was 1340, checked in by garnier, 15 years ago

update ti head
File size: 11.1 KB
Line 
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//
27// $Id: G4HEPlot.cc,v 1.11 2006/06/29 20:30:36 gunter Exp $
28// GEANT4 tag $Name: geant4-09-03-ref-09 $
29//
30//
31
32#include "globals.hh"
33#include "G4ios.hh"
34
35//
36// G4 Gheisha friend class G4GHEPlot
37// last modified: H. Fesefeldt 02-July--1998
38
39#include "G4HEPlot.hh"
40
41
42void
43G4HEPlot::Init( G4int nbin, G4double xstart, G4double xbin )
44 {
45 Xbin = xbin;
46 Nbin = nbin;
47 Xstart = xstart;
48 Entries = 0;
49 EntriesOverflow = 0;
50 EntriesUnderflow = 0;
51 Weight = 0.;
52 WeightOverflow = 0.;
53 WeightUnderflow = 0.;
54 Xvalue = new G4double[Nbin];
55 Yvalue = new G4double[Nbin];
56 for(G4int i=0; i<Nbin; i++)
57 {
58 Xvalue[i] = Xstart + i*Xbin;
59 Yvalue[i] = 0.;
60 }
61 return;
62 }
63
64void
65G4HEPlot::Add( G4double s1, G4double s2,
66 const G4HEPlot & p1, const G4HEPlot & p2 )
67 {
68 if(p1.Nbin != p2.Nbin)
69 {
70 G4cout << "G4HEPlot::Add: Plots must have same number of bins !" << G4endl;
71 return;
72 }
73 if(Nbin != p1.Nbin)
74 {
75 G4cout << "G4HEPlot::Add: Plot must be initialized before using it !" << G4endl;
76 return;
77 }
78
79 for(G4int i=0; i<Nbin; i++)
80 {
81 Yvalue[i] = s1*p1.Yvalue[i] + s2*p2.Yvalue[i];
82 }
83 Entries = p1.Entries + p2.Entries;
84 EntriesOverflow = p1.EntriesOverflow + p2.EntriesOverflow;
85 EntriesUnderflow = p1.EntriesUnderflow + p2.EntriesUnderflow;
86 Weight = s1*p1.Weight + s2*p2.Weight;
87 WeightUnderflow = s1*p1.WeightUnderflow + s2*p2.WeightUnderflow;
88 WeightOverflow = s1*p1.WeightOverflow + s2*p2.WeightOverflow;
89 return;
90 }
91
92void
93G4HEPlot::Multiply( G4double s1, G4double s2,
94 const G4HEPlot & p1, const G4HEPlot & p2 )
95 {
96 if(p1.Nbin != p2.Nbin)
97 {
98 G4cout << "G4HEPlot::Multiply: Plots must have same number of bins !" << G4endl;
99 return;
100 }
101 if(Nbin != p1.Nbin)
102 {
103 G4cout << "G4HEPlot::Multiply: Plot must be initialized before using it !" << G4endl;
104 return;
105 }
106
107 for(G4int i=0; i<Nbin; i++)
108 {
109 Yvalue[i] = s1*p1.Yvalue[i]* s2*p2.Yvalue[i];
110 }
111 Entries = p1.Entries * p2.Entries;
112 EntriesOverflow = p1.EntriesOverflow * p2.EntriesOverflow;
113 EntriesUnderflow = p1.EntriesUnderflow * p2.EntriesUnderflow;
114 Weight = s1*p1.Weight * s2*p2.Weight;
115 WeightUnderflow = s1*p1.WeightUnderflow * s2*p2.WeightUnderflow;
116 WeightOverflow = s1*p1.WeightOverflow * s2*p2.WeightOverflow;
117 return;
118 }
119
120void
121G4HEPlot::Divide( G4double s1, G4double s2,
122 const G4HEPlot & p1, const G4HEPlot & p2 )
123 {
124 if(p1.Nbin != p2.Nbin)
125 {
126 G4cout << "G4HEPlot::Divide: Plots must have same number of bins !" << G4endl;
127 return;
128 }
129 if(Nbin != p1.Nbin)
130 {
131 G4cout << "G4HEPlot::Divide: Plot must be defined before using it !" << G4endl;
132 return;
133 }
134
135 for(G4int i=0; i<Nbin; i++)
136 {
137 if(p2.Yvalue[i] == 0.)
138 {
139 Yvalue[i] = 0.;
140 }
141 else
142 {
143 Yvalue[i] = (s1*p1.Yvalue[i]) / (s2*p2.Yvalue[i]);
144 }
145 }
146 if(p2.Entries > 0)
147 Entries = p1.Entries / p2.Entries;
148 if(p2.EntriesOverflow > 0)
149 EntriesOverflow = p1.EntriesOverflow / p2.EntriesOverflow;
150 if(p2.EntriesUnderflow > 0)
151 EntriesUnderflow = p1.EntriesUnderflow / p2.EntriesUnderflow;
152 if(p2.Weight != 0.)
153 Weight = s1*p1.Weight / s2*p2.Weight;
154 if(p2.WeightUnderflow != 0.)
155 WeightUnderflow = s1*p1.WeightUnderflow / s2*p2.WeightUnderflow;
156 if(p2.WeightOverflow != 0.)
157 WeightOverflow = s1*p1.WeightOverflow / s2*p2.WeightOverflow;
158 return;
159 }
160
161void
162G4HEPlot::Scale( G4double s, const G4HEPlot & p)
163 {
164 if(Nbin != p.Nbin)
165 {
166 G4cout << "G4HEPlot::Add: Plot must be defined before using it !" << G4endl;
167 return;
168 }
169
170 for(G4int i=0; i<Nbin; i++)
171 {
172 Yvalue[i] = s*p.Yvalue[i];
173 }
174 Entries = p.Entries ;
175 EntriesOverflow = p.EntriesOverflow;
176 EntriesUnderflow = p.EntriesUnderflow;
177 Weight = s*p.Weight;
178 WeightUnderflow = s*p.WeightUnderflow;
179 WeightOverflow = s*p.WeightOverflow;
180 return;
181 }
182
183void
184G4HEPlot::XScale(G4double a, G4double b)
185 {
186 Xstart = b + Xstart;
187 Xbin = a*Xbin;
188 for(G4int i=0; i<Nbin; i++)
189 {
190 Xvalue[i] = Xstart + i*Xbin;
191 }
192 }
193
194void
195G4HEPlot::Shift(G4int nshift)
196 {
197 G4int i;
198 if(nshift < 0)
199 {
200 for(i=0;i<Nbin-nshift;i++)
201 {
202 Yvalue[i] = Yvalue[i-nshift];
203 }
204 for(i=Nbin-nshift;i<Nbin;i++)
205 {
206 Yvalue[i] = 0.;
207 }
208 }
209 if(nshift > 0)
210 {
211 for(i=Nbin-1;i>nshift-1;i--)
212 {
213 Yvalue[i] = Yvalue[i-nshift];
214 }
215 for(i=0;i<nshift;i++)
216 {
217 Yvalue[i] = 0.;
218 }
219 return;
220 }
221 }
222void
223G4HEPlot::Log( G4double s, const G4HEPlot & p)
224 {
225 if(Nbin != p.Nbin)
226 {
227 G4cout << "G4HEPlot::Log: Plot must be defined before using it !" << G4endl;
228 return;
229 }
230
231 for(G4int i=0; i<Nbin; i++)
232 {
233 if(s*p.Yvalue[i] <= 0.) Yvalue[i] = 0.;
234 else Yvalue[i] = std::log10(s*p.Yvalue[i]);
235 }
236
237 Entries = p.Entries ;
238 EntriesOverflow = p.EntriesOverflow;
239 EntriesUnderflow = p.EntriesUnderflow;
240 if(p.Weight > 0)
241 Weight = std::log10(s*p.Weight);
242 if(p.WeightUnderflow > 0)
243 WeightUnderflow = std::log10(s*p.WeightUnderflow);
244 if(p.WeightOverflow > 0)
245 WeightOverflow = std::log10(s*p.WeightOverflow);
246 return;
247 }
248
249void
250G4HEPlot::Sqrt(G4double s, const G4HEPlot & p)
251 {
252 if(Nbin != p.Nbin)
253 {
254 G4cout << " G4HEPlot::Sqrt: Plot must be defined before using it !" << G4endl;
255 return;
256 }
257 for (G4int i=0; i<Nbin; i++)
258 {
259 if(s*p.Yvalue[i] <= 0.) Yvalue[i] = 0.;
260 else Yvalue[i] = std::sqrt(s*p.Yvalue[i]);
261 }
262 Entries = p.Entries;
263 EntriesOverflow = p.EntriesOverflow;
264 EntriesUnderflow = p.EntriesUnderflow;
265 if(s*p.Weight > 0.) Weight = std::sqrt(s*p.Weight);
266 if(s*p.WeightUnderflow > 0.) WeightUnderflow = std::sqrt(s*p.WeightUnderflow);
267 if(s*p.WeightOverflow > 0.) WeightOverflow = std::sqrt(s*p.WeightOverflow);
268 return;
269 }
270
271void
272G4HEPlot::Reset()
273 {
274 for(G4int i=0; i<Nbin; i++)
275 {
276 Yvalue[i] = 0.;
277 }
278 Entries = 0;
279 EntriesOverflow = 0;
280 EntriesUnderflow = 0;
281 Weight = 0.;
282 WeightOverflow = 0.;
283 WeightUnderflow = 0.;
284 return;
285 }
286
287void
288G4HEPlot::LinearFit(G4double& a, G4double& b)
289 {
290 G4double a11 = 0.;
291 G4double a12 = 0.;
292 G4double a22 = 0.;
293 G4double b1 = 0.;
294 G4double b2 = 0.;
295 for(G4int i = 0; i<Nbin; i++)
296 {
297 if(Yvalue[i] != 0.)
298 {
299 a11 += (Xvalue[i]+Xbin/2.)*(Xvalue[i]+Xbin/2.);
300 a12 += (Xvalue[i]+Xbin/2.);
301 a22 += 1.;
302 b1 += Yvalue[i]*(Xvalue[i]+Xbin/2.);
303 b2 += Yvalue[i];
304 }
305 }
306 b = (b1*a12-b2*a11)/(a12*a12-a11*a22);
307 a = (b1 - b*a12)/a11;
308 return;
309 }
310
311void
312G4HEPlot::Fill(G4double x, G4double weight)
313 {
314 G4int i = (int)((x - Xstart)/Xbin);
315 if(i < 0)
316 {
317 EntriesUnderflow++;
318 WeightUnderflow += weight;
319 }
320 else if(i >= Nbin)
321 {
322 EntriesOverflow++;
323 WeightOverflow += weight;
324 }
325 else
326 {
327 Entries++;
328 Weight += weight;
329 Yvalue[i] += weight;
330 }
331 return;
332 }
333
334void
335G4HEPlot::Print( G4String name, G4int iplot)
336 {
337 G4cout << name << iplot << " = new TH1F(" << '"' << iplot
338 << '"' << "," << '"' << '"' << "," << Nbin << ","
339 << Xstart << "," << Xstart+Nbin*Xbin << ");" << G4endl;
340 G4cout << "for(I=1; I<=" << Nbin <<"; I++) {" << G4endl;
341 G4cout << "Y[I] = 0.;" << G4endl;
342 G4cout << "}" << G4endl;
343 for (G4int i=0; i<Nbin; i++)
344 {
345 G4cout << "Y[" << i+1 << "] = " << Yvalue[i] << ";" << G4endl;
346 }
347 G4cout << "XA = " << Xstart << ";" << G4endl;
348 G4cout << "STEP2 = " << Xbin/2. << ";" << G4endl;
349 G4cout << "for(I=1; I<=" << Nbin << "; I++) {" << G4endl;
350 G4cout << "XA = XA + STEP2;" << G4endl;
351 G4cout << name << iplot << "->Fill(XA,Y[I]);" << G4endl;
352 G4cout << "XA = XA + STEP2;" << G4endl;
353 G4cout << "}" << G4endl;
354 return;
355 }
356
357void
358G4HEPlot::DumpToFile(G4int aPlot, G4String aName)
359 {
360 FILE* fz;
361 fz = fopen(aName, "a");
362 fprintf(fz,"%3d %3d %12.5e %12.5e \n",aPlot, Nbin, Xstart, Xbin);
363 G4double xval = Xstart - Xbin/2.;
364 for(G4int i=0; i<Nbin; i++)
365 {
366 xval += Xbin;
367 fprintf(fz,"%12.5e %12.5e \n", xval, Yvalue[i]);
368 }
369 fclose(fz);
370 }
371
372void
373G4HEPlot::GetFromFile(G4int aPlot, G4String aName)
374 {
375 FILE* fz;
376 long int ip, nb;
377 G4double xs,xb,x,y;
378 if((fz = fopen(aName, "r")) != NULL)
379 {
380 while(fscanf(fz,"%ld %ld %le %le", &ip, &nb, &xs, &xb)==4)
381 {
382 if(ip == aPlot)
383 {
384 if(Nbin == 0) Init(nb, xs, xb);
385 }
386 for (G4int i=0; i<nb; i++)
387 {
388 fscanf(fz,"%le %le", &x, &y);
389 if(ip == aPlot) Fill(x,y);
390 }
391 if(ip == aPlot) break;
392 }
393 fclose(fz);
394 return;
395 }
396 else
397 {
398 G4cout << " File " << aName << " not found " << G4endl;
399 }
400 return;
401 }
402
403
404
Note: See TracBrowser for help on using the repository browser.