Geant4_10
G4H1ToolsManager.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 // $Id: G4H1ToolsManager.cc 70604 2013-06-03 11:27:06Z ihrivnac $
27 
28 // Author: Ivana Hrivnacova, 18/06/2013 (ivana@ipno.in2p3.fr)
29 
30 #include "G4H1ToolsManager.hh"
31 #include "G4HnManager.hh"
33 #include "G4AnalysisUtilities.hh"
34 
35 #include "tools/histo/h1d"
36 
37 #include <fstream>
38 
39 using namespace G4Analysis;
40 
41 //
42 // Constructors, destructor
43 //
44 
45 //_____________________________________________________________________________
47  : G4VH1Manager(state),
48  fH1Vector(),
49  fH1NameIdMap()
50 {
51 }
52 
53 //_____________________________________________________________________________
55 {
56  std::vector<tools::histo::h1d*>::iterator it;
57  for (it = fH1Vector.begin(); it != fH1Vector.end(); it++ ) {
58  delete (*it);
59  }
60 }
61 
62 //
63 // Utility functions
64 //
65 
66 namespace {
67 
68 //_____________________________________________________________________________
69 void UpdateH1Information(G4HnInformation* information,
70  const G4String& unitName,
71  const G4String& fcnName,
72  G4BinScheme binScheme)
73 {
74  G4double unit = GetUnitValue(unitName);
75  G4Fcn fcn = GetFunction(fcnName);
76  information->fXUnitName = unitName;
77  information->fYUnitName = unitName;
78  information->fXFcnName = fcnName;
79  information->fYFcnName = fcnName;
80  information->fXUnit = unit;
81  information->fYUnit = unit;
82  information->fXFcn = fcn;
83  information->fYFcn = fcn;
84  information->fXBinScheme = binScheme;
85  information->fYBinScheme = binScheme;
86 }
87 
88 //_____________________________________________________________________________
89 void AddH1Annotation(tools::histo::h1d* h1d,
90  const G4String& unitName,
91  const G4String& fcnName)
92 {
93  G4String axisTitle;
94  UpdateTitle(axisTitle, unitName, fcnName);
95  h1d->add_annotation(tools::histo::key_axis_x_title(), axisTitle);
96 }
97 
98 //_____________________________________________________________________________
99 tools::histo::h1d* CreateToolsH1(const G4String& title,
100  G4int nbins, G4double xmin, G4double xmax,
101  const G4String& fcnName,
102  const G4String& binSchemeName)
103 {
104  G4Fcn fcn = GetFunction(fcnName);
105  G4BinScheme binScheme = GetBinScheme(binSchemeName);
106 
107  if ( binScheme != kLogBinScheme ) {
108  if ( binScheme == kUserBinScheme ) {
109  // This should never happen, but let's make sure about it
110  // by issuing a warning
111  G4ExceptionDescription description;
112  description
113  << " User binning scheme setting was ignored." << G4endl
114  << " Linear binning will be applied with given (nbins, xmin, xmax) values";
115  G4Exception("G4H1ToolsManager::CreateH1",
116  "Analysis_W013", JustWarning, description);
117  }
118  return new tools::histo::h1d(title, nbins, fcn(xmin), fcn(xmax));
119  }
120  else {
121  // Compute edges
122  std::vector<G4double> edges;
123  ComputeEdges(nbins, xmin, xmax, fcn, binScheme, edges);
124  return new tools::histo::h1d(title, edges);
125  }
126 }
127 
128 //_____________________________________________________________________________
129 tools::histo::h1d* CreateToolsH1(const G4String& title,
130  const std::vector<G4double>& edges,
131  const G4String& fcnName)
132 {
133  G4Fcn fcn = GetFunction(fcnName);
134 
135  // Apply function
136  std::vector<G4double> newEdges;
137  ComputeEdges(edges, fcn, newEdges);
138 
139  return new tools::histo::h1d(title, newEdges);
140 }
141 
142 //_____________________________________________________________________________
143 void ConfigureToolsH1(tools::histo::h1d* h1d,
144  G4int nbins, G4double xmin, G4double xmax,
145  const G4String& fcnName,
146  const G4String& binSchemeName)
147 {
148  G4Fcn fcn = GetFunction(fcnName);
149  G4BinScheme binScheme = GetBinScheme(binSchemeName);
150 
151  if ( binScheme != kLogBinScheme ) {
152  if ( binScheme == kUserBinScheme ) {
153  // This should never happen, but let's make sure about it
154  // by issuing a warning
155  G4ExceptionDescription description;
156  description
157  << " User binning scheme setting was ignored." << G4endl
158  << " Linear binning will be applied with given (nbins, xmin, xmax) values";
159  G4Exception("G4H1ToolsManager::SetH1",
160  "Analysis_W013", JustWarning, description);
161  }
162  h1d->configure(nbins, fcn(xmin), fcn(xmax));
163  }
164  else {
165  // Compute bins
166  std::vector<G4double> edges;
167  ComputeEdges(nbins, xmin, xmax, fcn, binScheme, edges);
168  h1d->configure(edges);
169  }
170 }
171 
172 //_____________________________________________________________________________
173 void ConfigureToolsH1(tools::histo::h1d* h1d,
174  const std::vector<G4double>& edges,
175  const G4String& fcnName)
176 {
177  // Apply function to edges
178  G4Fcn fcn = GetFunction(fcnName);
179  std::vector<G4double> newEdges;
180  ComputeEdges(edges, fcn, newEdges);
181 
182  h1d->configure(newEdges);
183 }
184 
185 }
186 
187 //
188 // private methods
189 //
190 
191 //_____________________________________________________________________________
192 tools::histo::h1d* G4H1ToolsManager::GetH1InFunction(G4int id,
193  G4String functionName, G4bool warn,
194  G4bool onlyIfActive) const
195 {
196  G4int index = id - fFirstId;
197  if ( index < 0 || index >= G4int(fH1Vector.size()) ) {
198  if ( warn) {
199  G4String inFunction = "G4H1ToolsManager::";
200  inFunction += functionName;
201  G4ExceptionDescription description;
202  description << " " << "histogram " << id << " does not exist.";
203  G4Exception(inFunction, "Analysis_W007", JustWarning, description);
204  }
205  return 0;
206  }
207 
208  // Do not return histogram if inactive
209  if ( fState.GetIsActivation() && onlyIfActive && ( ! fHnManager->GetActivation(id) ) ) {
210  return 0;
211  }
212 
213  return fH1Vector[index];
214 }
215 
216 //_____________________________________________________________________________
217 void G4H1ToolsManager::AddH1Information(const G4String& name,
218  const G4String& unitName,
219  const G4String& fcnName,
220  G4BinScheme binScheme) const
221 {
222  G4double unit = GetUnitValue(unitName);
223  G4Fcn fcn = GetFunction(fcnName);
224  fHnManager->AddH1Information(name, unitName, fcnName, unit, fcn, binScheme);
225 }
226 
227 //_____________________________________________________________________________
228 G4int G4H1ToolsManager::RegisterToolsH1(tools::histo::h1d* h1d,
229  const G4String& name)
230 {
231  G4int index = fH1Vector.size();
232  fH1Vector.push_back(h1d);
233 
234  fLockFirstId = true;
235  fH1NameIdMap[name] = index + fFirstId;
236  return index + fFirstId;
237 }
238 
239 //
240 // protected methods
241 //
242 
243 //_____________________________________________________________________________
245  G4int nbins, G4double xmin, G4double xmax,
246  const G4String& unitName, const G4String& fcnName,
247  const G4String& binSchemeName)
248 {
249 #ifdef G4VERBOSE
250  if ( fState.GetVerboseL4() )
251  fState.GetVerboseL4()->Message("create", "H1", name);
252 #endif
253  tools::histo::h1d* h1d
254  = CreateToolsH1(title, nbins, xmin, xmax, fcnName, binSchemeName);
255 
256  // Add annotation
257  AddH1Annotation(h1d, unitName, fcnName);
258 
259  // Save H1 information
260  G4BinScheme binScheme = GetBinScheme(binSchemeName);
261  AddH1Information(name, unitName, fcnName, binScheme);
262 
263  // Register histogram
264  G4int id = RegisterToolsH1(h1d, name);
265 
266 #ifdef G4VERBOSE
267  if ( fState.GetVerboseL2() )
268  fState.GetVerboseL2()->Message("create", "H1", name);
269 #endif
270  return id;
271 }
272 
273 //_____________________________________________________________________________
275  const std::vector<G4double>& edges,
276  const G4String& unitName, const G4String& fcnName)
277 {
278 #ifdef G4VERBOSE
279  if ( fState.GetVerboseL4() )
280  fState.GetVerboseL4()->Message("create", "H1", name);
281 #endif
282  tools::histo::h1d* h1d
283  = CreateToolsH1(title, edges, fcnName);
284 
285  // Add annotation
286  AddH1Annotation(h1d, unitName, fcnName);
287 
288  // Save H1 information
289  AddH1Information(name, unitName, fcnName, kUserBinScheme);
290 
291  // Register histogram
292  G4int id = RegisterToolsH1(h1d, name);
293 
294 #ifdef G4VERBOSE
295  if ( fState.GetVerboseL2() )
296  fState.GetVerboseL2()->Message("create", "H1", name);
297 #endif
298  return id;
299 }
300 
301 //_____________________________________________________________________________
303  G4int nbins, G4double xmin, G4double xmax,
304  const G4String& unitName, const G4String& fcnName,
305  const G4String& binSchemeName)
306 {
307  tools::histo::h1d* h1d = GetH1InFunction(id, "SetH1", false, false);
308  if ( ! h1d ) return false;
309 
311 #ifdef G4VERBOSE
312  if ( fState.GetVerboseL4() )
313  fState.GetVerboseL4()->Message("configure", "H1", info->fName);
314 #endif
315 
316  // Configure tools h1
317  ConfigureToolsH1(h1d, nbins, xmin, xmax, fcnName, binSchemeName);
318 
319  // Add annotation
320  AddH1Annotation(h1d, unitName, fcnName);
321 
322  // Update information
323  G4BinScheme binScheme = GetBinScheme(binSchemeName);
324  UpdateH1Information(info, unitName, fcnName, binScheme);
325 
326  // Set activation
327  fHnManager->SetActivation(id, true);
328 
329  return true;
330 }
331 
332 //_____________________________________________________________________________
334  const std::vector<G4double>& edges,
335  const G4String& unitName,
336  const G4String& fcnName)
337 {
338  tools::histo::h1d* h1d = GetH1InFunction(id, "SetH1", false, false);
339  if ( ! h1d ) return false;
340 
342 #ifdef G4VERBOSE
343  if ( fState.GetVerboseL4() )
344  fState.GetVerboseL4()->Message("configure", "H1", info->fName);
345 #endif
346 
347  // Configure tools h1
348  ConfigureToolsH1(h1d, edges, fcnName);
349 
350  // Add annotation
351  AddH1Annotation(h1d, unitName, fcnName);
352 
353  // Update information
354  UpdateH1Information(info, unitName, fcnName, kUserBinScheme);
355 
356  // Set activation
357  fHnManager->SetActivation(id, true);
358 
359  return true;
360 }
361 
362 
363 //_____________________________________________________________________________
365 {
366  tools::histo::h1d* h1d = GetH1InFunction(id, "ScaleH1", false, false);
367  if ( ! h1d ) return false;
368 
369  return h1d->scale(factor);
370 }
371 
372 //_____________________________________________________________________________
374 {
375  tools::histo::h1d* h1d = GetH1InFunction(id, "FillH1", true, false);
376  if ( ! h1d ) return false;
377 
378  if ( fState.GetIsActivation() && ( ! fHnManager->GetActivation(id) ) ) {
379  //G4cout << "Skipping FillH1 for " << id << G4endl;
380  return false;
381  }
382 
384  h1d->fill(info->fXFcn(value/info->fXUnit), weight);
385 #ifdef G4VERBOSE
386  if ( fState.GetVerboseL4() ) {
387  G4ExceptionDescription description;
388  description << " id " << id << " value " << value;
389  fState.GetVerboseL4()->Message("fill", "H1", description);
390  }
391 #endif
392  return true;
393 }
394 
395 //_____________________________________________________________________________
397 {
398  std::map<G4String, G4int>::const_iterator it = fH1NameIdMap.find(name);
399  if ( it == fH1NameIdMap.end() ) {
400  if ( warn) {
401  G4String inFunction = "G4H1ToolsManager::GetH1Id";
402  G4ExceptionDescription description;
403  description << " " << "histogram " << name << " does not exist.";
404  G4Exception(inFunction, "Analysis_W007", JustWarning, description);
405  }
406  return -1;
407  }
408  return it->second;
409 }
410 
411 //_____________________________________________________________________________
413 {
414  tools::histo::h1d* h1d = GetH1InFunction(id, "GetH1Nbins");
415  if ( ! h1d ) return 0;
416 
417  return h1d->axis().bins();
418 }
419 
420 //_____________________________________________________________________________
422 {
423 // Returns xmin value with applied unit and histogram function
424 
425  tools::histo::h1d* h1d = GetH1InFunction(id, "GetH1Xmin");
426  if ( ! h1d ) return 0;
427 
428  G4HnInformation* info = fHnManager->GetHnInformation(id, "GetH1Xmin");
429  return info->fXFcn(h1d->axis().lower_edge()*info->fXUnit);
430 }
431 
432 //_____________________________________________________________________________
434 {
435  tools::histo::h1d* h1d = GetH1InFunction(id, "GetH1Xmax");
436  if ( ! h1d ) return 0;
437 
438  G4HnInformation* info = fHnManager->GetHnInformation(id, "GetH1Xmax");
439  return info->fXFcn(h1d->axis().upper_edge()*info->fXUnit);
440 }
441 
442 //_____________________________________________________________________________
444 {
445  tools::histo::h1d* h1d = GetH1InFunction(id, "GetH1XWidth", true, false);
446  if ( ! h1d ) return 0;
447 
448  G4int nbins = h1d->axis().bins();
449  if ( ! nbins ) {
450  G4ExceptionDescription description;
451  description << " nbins = 0 (for h1 id = " << id << ").";
452  G4Exception("G4H1ToolsManager::GetH1Width",
453  "Analysis_W014", JustWarning, description);
454  return 0;
455  }
456 
457  G4HnInformation* info = fHnManager->GetHnInformation(id, "GetH1XWidth");
458  return ( info->fXFcn(h1d->axis().upper_edge())
459  - info->fXFcn(h1d->axis().lower_edge()))*info->fXUnit/nbins;
460 }
461 
462 //_____________________________________________________________________________
464 {
465  tools::histo::h1d* h1d = GetH1InFunction(id, "SetH1Title");
466  if ( ! h1d ) return false;
467 
468  return h1d->set_title(title);
469 }
470 
471 //_____________________________________________________________________________
473 {
474  tools::histo::h1d* h1d = GetH1InFunction(id, "SetH1XAxisTitle");
475  if ( ! h1d ) return false;
476 
477  h1d->add_annotation(tools::histo::key_axis_x_title(), title);
478  return true;
479 }
480 
481 //_____________________________________________________________________________
483 {
484  tools::histo::h1d* h1d = GetH1InFunction(id, "SetH1YAxisTitle");
485  if ( ! h1d ) return false;
486 
487  h1d->add_annotation(tools::histo::key_axis_y_title(), title);
488  return true;
489 }
490 
491 //_____________________________________________________________________________
493 {
494  tools::histo::h1d* h1d = GetH1InFunction(id, "GetH1Title");
495  if ( ! h1d ) return "";
496 
497  return h1d->title();
498 }
499 
500 
501 //_____________________________________________________________________________
503 {
504  tools::histo::h1d* h1d = GetH1InFunction(id, "GetH1XAxisTitle");
505  if ( ! h1d ) return "";
506 
507  G4String title;
508  G4bool result = h1d->annotation(tools::histo::key_axis_x_title(), title);
509  if ( ! result ) {
510  G4ExceptionDescription description;
511  description << " Failed to get x_axis title for h1 id = " << id << ").";
512  G4Exception("G4H1ToolsManager::GetH1XAxisTitle",
513  "Analysis_W014", JustWarning, description);
514  return "";
515  }
516 
517  return title;
518 }
519 
520 //_____________________________________________________________________________
522 {
523  tools::histo::h1d* h1d = GetH1InFunction(id, "GetH1YAxisTitle");
524  if ( ! h1d ) return "";
525 
526  G4String title;
527  G4bool result = h1d->annotation(tools::histo::key_axis_y_title(), title);
528  if ( ! result ) {
529  G4ExceptionDescription description;
530  description << " Failed to get y_axis title for h1 id = " << id << ").";
531  G4Exception("G4H1ToolsManager::GetH1YAxisTitle",
532  "Analysis_W014", JustWarning, description);
533  return "";
534  }
535 
536  return title;
537 }
538 
539 //_____________________________________________________________________________
541 {
542 // Write selected objects on ASCII file
543 // According to the implementation by Michel Maire, originally in
544 // extended examples.
545 
546  // h1 histograms
547  for ( G4int i=0; i<G4int(fH1Vector.size()); ++i ) {
548  G4int id = i + fFirstId;
549  G4HnInformation* info = fHnManager->GetHnInformation(id,"WriteOnAscii");
550  // skip writing if activation is enabled and H1 is inactivated
551  if ( ! info->fAscii ) continue;
552  tools::histo::h1d* h1 = fH1Vector[i];
553 
554 #ifdef G4VERBOSE
555  if ( fState.GetVerboseL3() )
556  fState.GetVerboseL3()->Message("write on ascii", "h1d", info->fName);
557 #endif
558 
559  output << "\n 1D histogram " << id << ": " << h1->title()
560  << "\n \n \t X \t\t Y" << G4endl;
561 
562  for (G4int j=0; j< G4int(h1->axis().bins()); ++j) {
563  output << " " << j << "\t"
564  << h1->axis().bin_center(j) << "\t"
565  << h1->bin_height(j) << G4endl;
566  }
567  }
568 
569  return true;
570 }
571 
572 //
573 // public methods
574 //
575 
576 //_____________________________________________________________________________
578  const std::vector<tools::histo::h1d*>& h1Vector)
579 {
580 #ifdef G4VERBOSE
581  if ( fState.GetVerboseL4() )
582  fState.GetVerboseL4()->Message("merge", "all h1", "");
583 #endif
584  std::vector<tools::histo::h1d*>::const_iterator itw = h1Vector.begin();
585  std::vector<tools::histo::h1d*>::iterator it;
586  for (it = fH1Vector.begin(); it != fH1Vector.end(); it++ ) {
587  (*it)->add(*(*itw++));
588  }
589 #ifdef G4VERBOSE
590  if ( fState.GetVerboseL1() )
591  fState.GetVerboseL1()->Message("merge", "all h1", "");
592 #endif
593 }
594 
595 //_____________________________________________________________________________
597 {
598 // Reset histograms and ntuple
599 
600  G4bool finalResult = true;
601 
602  std::vector<tools::histo::h1d*>::iterator it;
603  for (it = fH1Vector.begin(); it != fH1Vector.end(); it++ ) {
604  G4bool result = (*it)->reset();
605  if ( ! result ) finalResult = false;
606  }
607 
608  return finalResult;
609 }
610 
611 //_____________________________________________________________________________
613 {
614  return ! fH1Vector.size();
615 }
616 
617 //_____________________________________________________________________________
619  G4bool onlyIfActive) const
620 {
621  return GetH1InFunction(id, "GetH1", warn, onlyIfActive);
622 }
623 
void Message(const G4String &action, const G4String &object, const G4String &objectName, G4bool success=true) const
virtual G4bool SetH1YAxisTitle(G4int id, const G4String &title)
virtual G4bool SetH1XAxisTitle(G4int id, const G4String &title)
virtual ~G4H1ToolsManager()
std::ostringstream G4ExceptionDescription
Definition: globals.hh:76
virtual G4bool SetH1(G4int id, G4int nbins, G4double xmin, G4double xmax, const G4String &unitName="none", const G4String &fcnName="none", const G4String &binSchemeName="linear")
Int_t index
Definition: macro.C:9
TH1F * h1
Definition: plot.C:43
virtual G4double GetH1Xmin(G4int id) const
void AddH1Information(const G4String &name, const G4String &unitName, const G4String &fcnName, G4double unit, G4Fcn fx, G4BinScheme binScheme)
Definition: G4HnManager.cc:57
tools::histo::h1d * GetH1(G4int id, G4bool warn=true, G4bool onlyIfActive=true) const
G4double(* G4Fcn)(G4double)
Definition: G4Fcn.hh:36
G4double G4NeutronHPJENDLHEData::G4double result
virtual G4bool SetH1Title(G4int id, const G4String &title)
const XML_Char * name
Definition: expat.h:151
double weight
Definition: plottest35.C:25
int G4int
Definition: G4Types.hh:78
G4HnInformation * GetHnInformation(G4int id, G4String functionName="", G4bool warn=true) const
Definition: G4HnManager.cc:86
G4bool GetActivation(G4int id) const
Definition: G4HnManager.cc:209
const G4AnalysisVerbose * GetVerboseL2() const
void UpdateTitle(G4String &title, const G4String &unitName, const G4String &fcnName)
virtual G4String GetH1XAxisTitle(G4int id) const
const G4AnalysisVerbose * GetVerboseL3() const
virtual G4int GetH1Nbins(G4int id) const
const G4AnalysisVerbose * GetVerboseL4() const
bool G4bool
Definition: G4Types.hh:79
virtual G4int GetH1Id(const G4String &name, G4bool warn=true) const
virtual G4bool WriteOnAscii(std::ofstream &output)
G4double GetUnitValue(const G4String &unit)
virtual G4int CreateH1(const G4String &name, const G4String &title, G4int nbins, G4double xmin, G4double xmax, const G4String &unitName="none", const G4String &fcnName="none", const G4String &binScheme="linear")
G4BinScheme fXBinScheme
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4HnManager * fHnManager
virtual G4double GetH1Width(G4int id) const
G4BinScheme fYBinScheme
G4Fcn GetFunction(const G4String &fcnName)
Definition: G4Fcn.cc:36
subroutine title(NA, NB, NCA, NCB)
Definition: dpm25nuc7.f:1744
virtual G4String GetH1YAxisTitle(G4int id) const
G4H1ToolsManager(const G4AnalysisManagerState &state)
G4BinScheme GetBinScheme(const G4String &binSchemeName)
Definition: G4BinScheme.cc:36
TH1D * h1d
Definition: berger.C:21
const XML_Char XML_Encoding * info
Definition: expat.h:530
const XML_Char int const XML_Char * value
Definition: expat.h:331
#define G4endl
Definition: G4ios.hh:61
void ComputeEdges(G4int nbins, G4double xmin, G4double xmax, G4Fcn fcn, G4BinScheme, std::vector< G4double > &edges)
Definition: G4BinScheme.cc:56
G4BinScheme
Definition: G4BinScheme.hh:40
void AddH1Vector(const std::vector< tools::histo::h1d * > &h1Vector)
virtual G4double GetH1Xmax(G4int id) const
double G4double
Definition: G4Types.hh:76
virtual G4bool ScaleH1(G4int id, G4double factor)
virtual G4String GetH1Title(G4int id) const
virtual G4bool FillH1(G4int id, G4double value, G4double weight=1.0)
const G4AnalysisVerbose * GetVerboseL1() const
void SetActivation(G4bool activation)
Definition: G4HnManager.cc:140
G4bool IsEmpty() const
const G4AnalysisManagerState & fState