| 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 /// \file electromagnetic/TestEm17/src/HistoMa << 26 // $Id: HistoManager.cc,v 1.4 2008/09/26 20:15:04 maire Exp $
27 /// \brief Implementation of the HistoManager << 27 // GEANT4 tag $Name: geant4-09-02-patch-04 $
28 // <<
29 // <<
30 // 28 //
31 //....oooOO0OOooo........oooOO0OOooo........oo 29 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
32 //....oooOO0OOooo........oooOO0OOooo........oo 30 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
33 31
34 #include "HistoManager.hh" 32 #include "HistoManager.hh"
35 <<
36 #include "HistoMessenger.hh" 33 #include "HistoMessenger.hh"
37 <<
38 #include "G4PhysicalConstants.hh" <<
39 #include "G4SystemOfUnits.hh" <<
40 #include "G4UnitsTable.hh" 34 #include "G4UnitsTable.hh"
41 35
>> 36 #ifdef G4ANALYSIS_USE
>> 37 #include "AIDA/AIDA.h"
>> 38 #endif
>> 39
42 //....oooOO0OOooo........oooOO0OOooo........oo 40 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
43 41
44 HistoManager::HistoManager() : fHistoMessenger << 42 HistoManager::HistoManager()
>> 43 :af(0),tree(0),hf(0),factoryOn(false)
45 { 44 {
>> 45 #ifdef G4ANALYSIS_USE
>> 46 // Creating the analysis factory
>> 47 af = AIDA_createAnalysisFactory();
>> 48 if(!af) {
>> 49 G4cout << " HistoManager::HistoManager() :"
>> 50 << " problem creating the AIDA analysis factory."
>> 51 << G4endl;
>> 52 }
>> 53 #endif
>> 54
46 fileName[0] = "testem17"; 55 fileName[0] = "testem17";
47 factoryOn = false; << 56 fileType = "root";
48 fNbHist = 0; << 57 fileOption = "--noErrors export=root uncompress";
49 <<
50 // histograms 58 // histograms
51 for (G4int k = 0; k < kMaxHisto; k++) { << 59 for (G4int k=0; k<MaxHisto; k++) {
52 fHistId[k] = 0; << 60 histo[k] = 0;
53 fHistPt[k] = 0; << 61 exist[k] = false;
54 fExist[k] = false; << 62 Unit[k] = 1.0;
55 fUnit[k] = 1.0; << 63 Width[k] = 1.0;
56 fWidth[k] = 1.0; <<
57 fAscii[k] = false; <<
58 } 64 }
59 65
60 fHistoMessenger = new HistoMessenger(this); << 66 histoMessenger = new HistoMessenger(this);
61 } 67 }
62 68
63 //....oooOO0OOooo........oooOO0OOooo........oo 69 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
64 70
65 HistoManager::~HistoManager() 71 HistoManager::~HistoManager()
66 { 72 {
67 delete fHistoMessenger; << 73 delete histoMessenger;
>> 74
>> 75 #ifdef G4ANALYSIS_USE
>> 76 delete af;
>> 77 #endif
68 } 78 }
69 79
70 //....oooOO0OOooo........oooOO0OOooo........oo 80 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
71 81
72 void HistoManager::Book() << 82 void HistoManager::book()
73 { << 83 {
74 // if no histos, do nothing << 84 #ifdef G4ANALYSIS_USE
75 if (fNbHist == 0) return; << 85 if(!af) return;
76 << 86
77 // Create or get analysis manager << 87 // Creating a tree mapped to an hbook file.
78 // The choice of analysis technology is done << 88 fileName[1] = fileName[0] + "." + fileType;
79 // in HistoManager.hh << 89 G4bool readOnly = false;
80 G4AnalysisManager* analysisManager = G4Analy << 90 G4bool createNew = true;
81 analysisManager->SetDefaultFileType("root"); << 91 AIDA::ITreeFactory* tf = af->createTreeFactory();
82 analysisManager->SetVerboseLevel(0); << 92 tree = tf->create(fileName[1], fileType, readOnly, createNew, fileOption);
83 G4String extension = analysisManager->GetFil << 93 delete tf;
84 fileName[1] = fileName[0] + "." + extension; << 94 if(!tree) {
85 << 95 G4cout << "HistoManager::book() :"
86 // Open an output file << 96 << " problem creating the AIDA tree with "
87 // << 97 << " storeName = " << fileName[1]
88 G4bool fileOpen = analysisManager->OpenFile( << 98 << " storeType = " << fileType
89 if (!fileOpen) { << 99 << " readOnly = " << readOnly
90 G4cout << "\n---> HistoManager::book(): ca << 100 << " createNew = " << createNew
>> 101 << " options = " << fileOption
>> 102 << G4endl;
91 return; 103 return;
92 } 104 }
93 105
94 // create selected histograms << 106 // Creating a histogram factory, whose histograms will be handled by the tree
95 // << 107 hf = af->createHistogramFactory(*tree);
96 analysisManager->SetFirstHistoId(1); <<
97 108
98 for (G4int k = 0; k < kMaxHisto; k++) { << 109 // create selected histograms
99 if (fExist[k]) { << 110 for (G4int k=0; k<MaxHisto; k++) {
100 fHistId[k] = analysisManager->CreateH1(f << 111 if (exist[k]) {
101 fHistPt[k] = analysisManager->GetH1(fHis << 112 histo[k] = hf->createHistogram1D( Label[k], Title[k],
>> 113 Nbins[k], Vmin[k], Vmax[k]);
102 factoryOn = true; 114 factoryOn = true;
103 } 115 }
104 } 116 }
>> 117 if(factoryOn)
>> 118 G4cout << "\n----> Histogram Tree is opened in " << fileName[1] << G4endl;
105 119
106 if (factoryOn) G4cout << "\n----> Histogram << 120 #endif
107 } 121 }
108 122
109 //....oooOO0OOooo........oooOO0OOooo........oo 123 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
110 124
111 void HistoManager::Save() << 125 void HistoManager::save()
112 { 126 {
>> 127 #ifdef G4ANALYSIS_USE
113 if (factoryOn) { 128 if (factoryOn) {
114 G4cout << "\n----> HistoManager::save " << << 129 tree->commit(); // Writing the histograms to the file
115 G4AnalysisManager* analysisManager = G4Ana << 130 tree->close(); // and closing the tree (and the file)
116 analysisManager->Write(); << 131 G4cout << "\n----> Histogram Tree is saved in " << fileName[1] << G4endl;
117 analysisManager->CloseFile(); << 132
118 SaveAscii(); // Write fAscii file, if any << 133 delete hf;
119 G4cout << "\n----> Histograms are saved in << 134 delete tree;
120 << 135 tree = 0;
121 analysisManager->Clear(); <<
122 factoryOn = false; 136 factoryOn = false;
123 } 137 }
>> 138 #endif
124 } 139 }
125 140
126 //....oooOO0OOooo........oooOO0OOooo........oo 141 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
127 142
128 void HistoManager::FillHisto(G4int ih, G4doubl 143 void HistoManager::FillHisto(G4int ih, G4double e, G4double weight)
129 { 144 {
130 if (ih > kMaxHisto) { << 145 if (ih > MaxHisto) {
131 G4cout << "---> warning from HistoManager: 146 G4cout << "---> warning from HistoManager::FillHisto() : histo " << ih
132 << "does not fExist; e= " << e << " << 147 << "does not exist; e= " << e << " w= " << weight << G4endl;
133 return; 148 return;
134 } 149 }
135 << 150 #ifdef G4ANALYSIS_USE
136 if (fHistPt[ih]) fHistPt[ih]->fill(e / fUnit << 151 if(exist[ih]) histo[ih]->fill(e/Unit[ih], weight);
>> 152 #endif
137 } 153 }
138 154
139 //....oooOO0OOooo........oooOO0OOooo........oo 155 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
140 156
141 void HistoManager::SetHisto(G4int ih, G4int nb << 157 void HistoManager::SetHisto(G4int ih,
142 const G4String& un << 158 G4int nbins, G4double valmin, G4double valmax, const G4String& unit)
143 { 159 {
144 if (ih > kMaxHisto) { << 160 if (ih > MaxHisto) {
145 G4cout << "---> warning from HistoManager: << 161 G4cout << "---> warning from HistoManager::SetHisto() : histo " << ih
146 << G4endl; << 162 << "does not exist" << G4endl;
147 return; 163 return;
148 } 164 }
149 << 165
150 const G4String id[] = {"0", "1", "2", "3" << 166 const G4String id[] = { "0", "1", "2", "3", "4", "5", "6", "7"};
151 "10", "11", "12", "13 << 167 const G4String title[] =
152 "20", "21", "22", "23 << 168 { "dummy", //0
153 const G4String title[] = { << 169 "log10(Etransfert/Emu) Ionization", //1
154 "dummy", // 0 << 170 "log10(Etransfert/Emu) Pair", //2
155 "log10(Eloss/Emu) muIonization", // 1 << 171 "log10(Etransfert/Emu) Brems", //3
156 "log10(Eloss/Emu) muPair", // 2 << 172 "log10(Etransfert/Emu) Nuclear", //4
157 "log10(Eloss/Emu) muBrems", // 3 << 173 "log10(Etransfert/Emu) Ionization", //5
158 "log10(Eloss/Emu) muNuclear", // 4 << 174 "log10(Etransfert/Emu) Pair", //6
159 "log10(Eloss/Emu) hIonization", // 5 << 175 "log10(Etransfert/Emu) Brems" //7
160 "log10(Eloss/Emu) hPair", // 6 << 176 };
161 "log10(Eloss/Emu) hBrems", // 7 <<
162 "log10(Eloss/Emu) muToMuonPair", // 8 <<
163 "dummy", // 9 <<
164 "dummy", // 10 <<
165 "log10(Eloss/Emu) muIonization", // 11 <<
166 "log10(Eloss/Emu) muPair", // 12 <<
167 "log10(Eloss/Emu) muBrems", // 13 <<
168 "log10(Eloss/Emu) muNuclear", // 14 <<
169 "dummy", // 15 <<
170 "dummy", // 16 <<
171 "dummy", // 17 <<
172 "log10(Eloss/Emu) muToMuonPair", // 18 <<
173 "dummy", // 19 <<
174 "dummy", // 20 <<
175 "CS(1/mm) vs ekin muIonisation", // 21 <<
176 "CS(1/mm) vs ekin muPair", // 22 <<
177 "CS(1/mm) vs ekin muBrems", // 23 <<
178 "CS(1/mm) vs ekin muNuclear", // 24 <<
179 "dummy", // 25 <<
180 "dummy", // 26 <<
181 "dummy", // 27 <<
182 "CS(1/mm) vs ekin muToMuonPair", // 28 <<
183 "dummy" // 29 <<
184 }; <<
185 177
186 G4String titl = title[ih]; 178 G4String titl = title[ih];
187 fUnit[ih] = 1.; << 179 G4double vmin = valmin, vmax = valmax;
>> 180 Unit[ih] = 1.;
188 181
189 if (unit != "none") { 182 if (unit != "none") {
190 titl = title[ih] + " (" + unit + ")"; 183 titl = title[ih] + " (" + unit + ")";
191 fUnit[ih] = G4UnitDefinition::GetValueOf(u << 184 Unit[ih] = G4UnitDefinition::GetValueOf(unit);
>> 185 vmin = valmin/Unit[ih]; vmax = valmax/Unit[ih];
192 } 186 }
193 187
194 fExist[ih] = true; << 188 exist[ih] = true;
195 fLabel[ih] = "h" + id[ih]; << 189 Label[ih] = id[ih];
196 fTitle[ih] = titl; << 190 Title[ih] = titl;
197 fNbins[ih] = nbins; << 191 Nbins[ih] = nbins;
198 fVmin[ih] = vmin; << 192 Vmin[ih] = vmin;
199 fVmax[ih] = vmax; << 193 Vmax[ih] = vmax;
200 fWidth[ih] = fUnit[ih] * (vmax - vmin) / nbi << 194 Width[ih] = (valmax-valmin)/nbins;
201 195
202 fNbHist++; << 196 G4cout << "----> SetHisto " << ih << ": " << titl << "; "
>> 197 << nbins << " bins from "
>> 198 << vmin << " " << unit << " to " << vmax << " " << unit << G4endl;
203 199
204 G4cout << "----> SetHisto " << ih << ": " << <<
205 << " " << unit << " to " << vmax << " <<
206 } 200 }
207 201
208 //....oooOO0OOooo........oooOO0OOooo........oo 202 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
209 203
210 void HistoManager::Normalize(G4int ih, G4doubl << 204 void HistoManager::RemoveHisto(G4int ih)
211 { 205 {
212 if (ih >= kMaxHisto) { << 206 if (ih > MaxHisto) {
213 G4cout << "---> warning from HistoManager: << 207 G4cout << "---> warning from HistoManager::RemoveHisto() : histo " << ih
214 << G4endl; << 208 << "does not exist" << G4endl;
215 return; 209 return;
216 } 210 }
217 211
218 if (fHistPt[ih]) fHistPt[ih]->scale(fac); << 212 histo[ih] = 0; exist[ih] = false;
219 } 213 }
220 214
221 //....oooOO0OOooo........oooOO0OOooo........oo 215 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
222 216
223 void HistoManager::PrintHisto(G4int ih) <<
224 { <<
225 if (ih < kMaxHisto) { <<
226 fAscii[ih] = true; <<
227 fAscii[0] = true; <<
228 } <<
229 else <<
230 G4cout << "---> warning from HistoManager: <<
231 << G4endl; <<
232 } <<
233 <<
234 //....oooOO0OOooo........oooOO0OOooo........oo <<
235 <<
236 #include <fstream> <<
237 <<
238 void HistoManager::SaveAscii() <<
239 { <<
240 if (!fAscii[0]) return; <<
241 <<
242 G4String name = fileName[0] + ".ascii"; <<
243 std::ofstream File(name, std::ios::out); <<
244 if (!File) { <<
245 G4cout << "\n---> HistoManager::saveAscii( <<
246 return; <<
247 } <<
248 <<
249 File.setf(std::ios::scientific, std::ios::fl <<
250 <<
251 // write selected histograms <<
252 for (G4int ih = 0; ih < kMaxHisto; ih++) { <<
253 if (fHistPt[ih] && fAscii[ih]) { <<
254 File << "\n 1D histogram " << ih << ": <<
255 << G4endl; <<
256 <<
257 for (G4int iBin = 0; iBin < fNbins[ih]; <<
258 File << " " << iBin << "\t" << fHistP <<
259 << fHistPt[ih]->bin_height(iBin) <<
260 } <<
261 } <<
262 } <<
263 } <<
264 <<
265 //....oooOO0OOooo........oooOO0OOooo........oo <<
266 217