| 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 hadronic/Hadr01/src/Histo.cc <<
27 /// \brief Implementation of the Histo class <<
28 // <<
29 // <<
30 //-------------------------------------------- 26 //---------------------------------------------------------------------------
31 // 27 //
32 // ClassName: Histo - Generic histogram/ntup 28 // ClassName: Histo - Generic histogram/ntuple manager class
33 // 29 //
34 // 30 //
35 // Author: V.Ivanchenko 30.10.03 31 // Author: V.Ivanchenko 30.10.03
36 // 32 //
37 //-------------------------------------------- 33 //----------------------------------------------------------------------------
38 // 34 //
39 35
40 //....oooOO0OOooo........oooOO0OOooo........oo << 36 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
41 //....oooOO0OOooo........oooOO0OOooo........oo << 37 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
42 38
43 #include "Histo.hh" 39 #include "Histo.hh"
44 40
>> 41 #ifdef G4ANALYSIS_USE
>> 42 #include <memory> // for the auto_ptr(T>
>> 43 #include <AIDA/AIDA.h>
45 #include "HistoMessenger.hh" 44 #include "HistoMessenger.hh"
>> 45 #endif
46 46
47 #include "G4RootAnalysisManager.hh" << 47 //#include <iomanip>
48 48
49 //....oooOO0OOooo........oooOO0OOooo........oo << 49 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
50 //....oooOO0OOooo........oooOO0OOooo........oo << 50 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
51 51
52 Histo::Histo() << 52 Histo::Histo(G4int ver)
53 : fManager(nullptr), <<
54 fHistName("test"), <<
55 fTupleName("tuple"), <<
56 fTupleTitle("test"), <<
57 fNHisto(0), <<
58 fVerbose(0), <<
59 fDefaultAct(true), <<
60 fHistoActive(false), <<
61 fNtupleActive(false) <<
62 { 53 {
63 fMessenger = new HistoMessenger(this); << 54 verbose = ver;
>> 55 histName = "histo";
>> 56 histType = "hbook";
>> 57 option = "--noErrors uncompress";
>> 58 nHisto = 0;
>> 59 defaultAct = 1;
>> 60 tupleName = "tuple.hbook";
>> 61 tupleId = "100";
>> 62 tupleList = "";
>> 63 ntup = 0;
>> 64 messenger = 0;
>> 65
>> 66 #ifdef G4ANALYSIS_USE
>> 67 messenger = new HistoMessenger(this);
>> 68 #endif
64 } 69 }
65 70
66 //....oooOO0OOooo........oooOO0OOooo........oo << 71 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
67 72
68 Histo::~Histo() 73 Histo::~Histo()
69 { 74 {
70 delete fMessenger; << 75 #ifdef G4ANALYSIS_USE
>> 76 for(G4int i=0; i<nHisto; i++) {
>> 77 if(histo[i]) delete histo[i];
>> 78 }
>> 79 delete messenger;
>> 80 #endif
71 } 81 }
72 82
73 //....oooOO0OOooo........oooOO0OOooo........oo << 83 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
74 84
75 void Histo::Book() << 85 void Histo::book()
76 { 86 {
77 if (!(fHistoActive || fNtupleActive)) { << 87 #ifdef G4ANALYSIS_USE
78 return; << 88 G4cout << "### Histo books " << nHisto << " histograms " << G4endl;
79 } << 89 // Creating the analysis factory
>> 90 std::auto_ptr< AIDA::IAnalysisFactory > af( AIDA_createAnalysisFactory() );
>> 91 if(verbose>0)
>> 92 G4cout<<"HIsto books analysis factory ......... "<<G4endl;
>> 93
>> 94 // Creating the tree factory
>> 95 std::auto_ptr< AIDA::ITreeFactory > tf( af->createTreeFactory() );
>> 96 if(verbose>0)
>> 97 G4cout<<"Histo books tree factory ......... "<<G4endl;
80 98
81 // Always creating analysis manager << 99 G4String histExt = "";
82 fManager = G4RootAnalysisManager::Instance() << 100 char* path = getenv("PHYSLIST");
>> 101 if (path) histExt = "_" + G4String(path);
>> 102
>> 103 G4String histDir = "";
>> 104 path = getenv("HISTODIR");
>> 105 if (path) histDir = G4String(path) + "/";
83 106
84 // Creating a tree mapped to a new hbook fil 107 // Creating a tree mapped to a new hbook file.
85 G4String nam = fHistName + ".root"; << 108 G4String name = histDir + histName + histExt + "." + histType;
>> 109 tree = tf->create(name, histType, false, true, option);
>> 110 G4cout << "Histo: tree store : " << tree->storeName() << G4endl;
86 111
87 // Open file histogram file << 112 // Creating a histogram factory, whose histograms will be handled by the tree
88 if (!fManager->OpenFile(nam)) { << 113 std::auto_ptr< AIDA::IHistogramFactory > hf(af->createHistogramFactory( *tree ));
89 G4cout << "Histo::Book: ERROR open file <" <<
90 fHistoActive = false; <<
91 fNtupleActive = false; <<
92 return; <<
93 } <<
94 G4cout << "### Histo::Save: Opended file <" <<
95 << G4endl; <<
96 114
97 // Creating an 1-dimensional histograms in t 115 // Creating an 1-dimensional histograms in the root directory of the tree
98 for (G4int i = 0; i < fNHisto; ++i) { << 116 for(G4int i=0; i<nHisto; i++) {
99 if (fActive[i]) { << 117 if(active[i]) {
100 G4String ss = "h" + fIds[i]; << 118 if(verbose>0)
101 fHisto[i] = fManager->CreateH1(ss, fTitl << 119 G4cout<<" I am in book: histogram "<< i << " id= " << ids[i] <<G4endl;
102 if (fVerbose > 0) { << 120
103 G4cout << "Created histogram #" << i < << 121 histo[i] = hf->createHistogram1D(ids[i], tittles[i], bins[i], xmin[i], xmax[i]);
104 << fTitles[i] << G4endl; <<
105 } <<
106 } 122 }
107 } 123 }
108 // Creating a tuple factory, whose tuples wi 124 // Creating a tuple factory, whose tuples will be handled by the tree
109 if (fNtupleActive) { << 125 if(tupleList != "") {
110 fManager->CreateNtuple(fTupleName, fTupleT << 126 if(verbose>0)
111 G4int i; << 127 G4cout<<"Histo books tuple factory for "<<tupleName <<G4endl;
112 G4int n = fNtupleI.size(); << 128 std::auto_ptr< AIDA::ITupleFactory > tpf( af->createTupleFactory( *tree ) );
113 for (i = 0; i < n; ++i) { << 129 ntup = tpf->create(tupleId, tupleName, tupleList);
114 if (fTupleI[i] == -1) { <<
115 fTupleI[i] = fManager->CreateNtupleICo <<
116 } <<
117 } <<
118 n = fNtupleF.size(); <<
119 for (i = 0; i < n; ++i) { <<
120 if (fTupleF[i] == -1) { <<
121 fTupleF[i] = fManager->CreateNtupleFCo <<
122 } <<
123 } <<
124 n = fNtupleD.size(); <<
125 for (i = 0; i < n; ++i) { <<
126 if (fTupleD[i] == -1) { <<
127 fTupleD[i] = fManager->CreateNtupleDCo <<
128 } <<
129 } <<
130 } 130 }
131 } << 131 #endif
>> 132 }
132 133
133 //....oooOO0OOooo........oooOO0OOooo........oo << 134 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
134 135
135 void Histo::Save() << 136 void Histo::save()
136 { 137 {
137 if (!(fHistoActive || fNtupleActive)) { << 138 #ifdef G4ANALYSIS_USE
138 return; <<
139 } <<
140 <<
141 // Write histogram file 139 // Write histogram file
142 if (!fManager->Write()) { << 140 tree->commit();
143 G4cout << "Histo::Save: FATAL ERROR writin << 141 G4cout << "Closing the tree..." << G4endl;
144 exit(1); << 142 tree->close();
145 } << 143 G4cout << "Histograms and Ntuples are saved" << G4endl;
146 if (fVerbose > 0) { << 144 for(G4int i=0; i<nHisto; i++) {
147 G4cout << "### Histo::Save: Histograms and << 145 if(histo[i]) histo[i]->reset();
148 } 146 }
149 if (fManager->CloseFile() && fVerbose > 0) { << 147 #endif
150 G4cout << "Histo::Save: File is closed" << <<
151 } <<
152 fManager->Clear(); <<
153 } 148 }
154 149
155 //....oooOO0OOooo........oooOO0OOooo........oo << 150 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
156 151
157 void Histo::Add1D(const G4String& id, const G4 << 152 void Histo::reset()
158 G4double u) <<
159 { 153 {
160 if (fVerbose > 0) { << 154 #ifdef G4ANALYSIS_USE
161 G4cout << "Histo::Add1D: New histogram wil << 155 for(G4int i=0; i<nHisto; i++) {
162 << " " << x1 << " " << x2 << " " << 156 if(histo[i]) histo[i]->reset();
163 } <<
164 ++fNHisto; <<
165 x1 /= u; <<
166 x2 /= u; <<
167 fActive.push_back(fDefaultAct); <<
168 fBins.push_back(nb); <<
169 fXmin.push_back(x1); <<
170 fXmax.push_back(x2); <<
171 fUnit.push_back(u); <<
172 fIds.push_back(id); <<
173 fTitles.push_back(name); <<
174 fHisto.push_back(-1); <<
175 } <<
176 <<
177 //....oooOO0OOooo........oooOO0OOooo........oo <<
178 <<
179 void Histo::SetHisto1D(G4int i, G4int nb, G4do <<
180 { <<
181 if (i >= 0 && i < fNHisto) { <<
182 if (fVerbose > 0) { <<
183 G4cout << "Histo::SetHisto1D: #" << i << <<
184 << G4endl; <<
185 } <<
186 fBins[i] = nb; <<
187 fXmin[i] = x1; <<
188 fXmax[i] = x2; <<
189 fUnit[i] = u; <<
190 fActive[i] = true; <<
191 fHistoActive = true; <<
192 } <<
193 else { <<
194 G4cout << "Histo::SetHisto1D: WARNING! wro <<
195 } 157 }
>> 158 #endif
196 } 159 }
197 160
198 //....oooOO0OOooo........oooOO0OOooo........oo << 161 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
199 162
200 void Histo::Activate(G4int i, G4bool val) << 163 void Histo::setFileType(const G4String& nam)
201 { 164 {
202 if (fVerbose > 1) { << 165 if(nam == "hbook" || nam == "root" || nam == "aida") histType = nam;
203 G4cout << "Histo::Activate: Histogram: #" << 166 else if(nam == "XML" || nam == "xml") histType = "aida";
204 } <<
205 if (i >= 0 && i < fNHisto) { <<
206 fActive[i] = val; <<
207 if (val) { <<
208 fHistoActive = true; <<
209 } <<
210 } <<
211 } 167 }
212 168
213 //....oooOO0OOooo........oooOO0OOooo........oo << 169 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
214 170
215 void Histo::Fill(G4int i, G4double x, G4double << 171 void Histo::add1D(const G4String& id, const G4String& name, G4int nb,
>> 172 G4double x1, G4double x2, G4double u)
216 { 173 {
217 if (!fHistoActive) { << 174 if(nHisto > 0) {
218 return; << 175 for(G4int i=0; i<nHisto; i++) {
219 } << 176 if(ids[i] == id) return;
220 if (fVerbose > 1) { <<
221 G4cout << "Histo::Fill: Histogram: #" << i <<
222 } <<
223 if (i >= 0 && i < fNHisto) { <<
224 if (fActive[i]) { <<
225 fManager->FillH1(fHisto[i], x / fUnit[i] <<
226 } 177 }
227 } 178 }
228 else { <<
229 G4cout << "Histo::Fill: WARNING! wrong his <<
230 } <<
231 } <<
232 179
233 //....oooOO0OOooo........oooOO0OOooo........oo << 180 if(verbose > 0)
>> 181 G4cout << "New histogram will be booked: #" << id << " <" << name
>> 182 << " " << nb << " " << x1 << " " << x2 << " " << u
>> 183 << G4endl;
234 184
235 void Histo::ScaleH1(G4int i, G4double x) << 185 nHisto++;
236 { << 186 x1 /= u;
237 if (!fHistoActive) { << 187 x2 /= u;
238 return; << 188 active.push_back(defaultAct);
239 } << 189 bins.push_back(nb);
240 if (fVerbose > 0) { << 190 xmin.push_back(x1);
241 G4cout << "Histo::Scale: Histogram: #" << << 191 xmax.push_back(x2);
242 } << 192 unit.push_back(u);
243 if (i >= 0 && i < fNHisto) { << 193 ids.push_back(id);
244 if (fActive[i]) { << 194 tittles.push_back(name);
245 fManager->GetH1(fHisto[i])->scale(x); << 195 histo.push_back(0);
>> 196 }
>> 197
>> 198 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
>> 199
>> 200 void Histo::setHisto1D(G4int i, G4int nb, G4double x1, G4double x2, G4double u)
>> 201 {
>> 202 if(i>=0 && i<nHisto) {
>> 203 if(verbose > 0)
>> 204 {
>> 205 G4cout << "Update histogram: #" << i
>> 206 << " " << nb << " " << x1 << " " << x2 << " " << u
>> 207 << G4endl;
246 } 208 }
247 } << 209 bins[i] = nb;
248 else { << 210 xmin[i] = x1;
249 G4cout << "Histo::Scale: WARNING! wrong hi << 211 xmax[i] = x2;
>> 212 unit[i] = u;
>> 213 } else {
>> 214 G4cout << "Histo::setHisto1D: WARNING! wrong histogram index " << i << G4endl;
>> 215 G4cout << "nHisto = " << nHisto << G4endl;
250 } 216 }
251 } 217 }
252 218
253 //....oooOO0OOooo........oooOO0OOooo........oo << 219 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
254 220
255 void Histo::AddTuple(const G4String& w1) << 221 void Histo::fill(G4int i, G4double x, G4double w)
256 { 222 {
257 fTupleTitle = w1; << 223 if(verbose > 1)
>> 224 G4cout << "fill histogram: #" << i << " at x= " << x
>> 225 << " weight= " << w << " unit= " << unit[i]
>> 226 << G4endl;
>> 227 #ifdef G4ANALYSIS_USE
>> 228 if(i>=0 && i<nHisto) {
>> 229 histo[i]->fill((float)(x/unit[i]), (float)w);
>> 230 } else {
>> 231 G4cout << "Histo::fill: WARNING! wrong histogram index " << i << G4endl;
>> 232 }
>> 233 #endif
258 } 234 }
259 235
260 //....oooOO0OOooo........oooOO0OOooo........oo << 236 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
261 237
262 void Histo::AddTupleI(const G4String& w1) << 238 void Histo::scale(G4int i, G4double x)
263 { 239 {
264 fNtupleActive = true; << 240 if(verbose > 0)
265 fNtupleI.push_back(w1); << 241 G4cout << "Scale histogram: #" << i << " by factor " << x << G4endl;
266 fTupleI.push_back(-1); <<
267 } <<
268 <<
269 //....oooOO0OOooo........oooOO0OOooo........oo <<
270 242
271 void Histo::AddTupleF(const G4String& w1) << 243 #ifdef G4ANALYSIS_USE
272 { << 244 if(i>=0 && i<nHisto) {
273 fNtupleActive = true; << 245 histo[i]->scale(x);
274 fNtupleF.push_back(w1); << 246 } else {
275 fTupleF.push_back(-1); << 247 G4cout << "Histo::scale: WARNING! wrong histogram index " << i << G4endl;
>> 248 }
>> 249 #endif
276 } 250 }
277 251
278 //....oooOO0OOooo........oooOO0OOooo........oo << 252 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
279 253
280 void Histo::AddTupleD(const G4String& w1) << 254 void Histo::addTuple(const G4String& w1, const G4String& w2, const G4String& w3)
281 { 255 {
282 fNtupleActive = true; << 256 tupleId = w1;
283 fNtupleD.push_back(w1); << 257 tupleName = w2;
284 fTupleD.push_back(-1); << 258 tupleList = w3;
285 } <<
286 <<
287 //....oooOO0OOooo........oooOO0OOooo........oo <<
288 259
289 void Histo::FillTupleI(G4int i, G4int x) << 260 G4cout<<" addTuple: Id "<<w1<<" Name "<<w2<<" List "<<w3<<G4endl;
290 { <<
291 if (!fNtupleActive) { <<
292 return; <<
293 } <<
294 G4int n = fNtupleI.size(); <<
295 if (i >= 0 && i < n) { <<
296 if (fVerbose > 1) { <<
297 G4cout << "Histo::FillTupleI: i= " << i <<
298 << "> = " << x << G4endl; <<
299 } <<
300 fManager->FillNtupleIColumn(fTupleI[i], x) <<
301 } <<
302 else { <<
303 G4cout << "Histo::FillTupleI: WARNING! wro <<
304 } <<
305 } 261 }
306 262
307 //....oooOO0OOooo........oooOO0OOooo........oo << 263 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
308 264
309 void Histo::FillTupleF(G4int i, G4float x) << 265 void Histo::fillTuple(const G4String& parname, G4double x)
310 { 266 {
311 if (!fNtupleActive) { << 267 if(verbose > 1)
312 return; << 268 G4cout << "fill tuple by parameter <" << parname << "> = " << x << G4endl;
313 } << 269
314 G4int n = fNtupleF.size(); << 270 #ifdef G4ANALYSIS_USE
315 if (i >= 0 && i < n) { << 271 if(ntup) ntup->fill(ntup->findColumn(parname), (float)x);
316 if (fVerbose > 1) { << 272 #endif
317 G4cout << "Histo::FillTupleF: i= " << i <<
318 << "> = " << x << G4endl; <<
319 } <<
320 fManager->FillNtupleFColumn(fTupleF[i], x) <<
321 } <<
322 else { <<
323 G4cout << "Histo::FillTupleF: WARNING! wro <<
324 } <<
325 } 273 }
326 274
327 //....oooOO0OOooo........oooOO0OOooo........oo << 275 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
328 276
329 void Histo::FillTupleD(G4int i, G4double x) << 277 void Histo::addRow()
330 { 278 {
331 if (!fNtupleActive) { << 279 #ifdef G4ANALYSIS_USE
332 return; << 280 if(ntup) ntup->addRow();
333 } << 281 #endif
334 G4int n = fNtupleD.size(); << 282 }
335 if (i >= 0 && i < n) { <<
336 if (fVerbose > 1) { <<
337 G4cout << "Histo::FillTupleD: i= " << i <<
338 << "> = " << x << G4endl; <<
339 } <<
340 fManager->FillNtupleDColumn(fTupleD[i], x) <<
341 } <<
342 else { <<
343 G4cout << "Histo::FillTupleD: WARNING! wro <<
344 } <<
345 } <<
346 283
347 //....oooOO0OOooo........oooOO0OOooo........oo << 284 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
348 << 285 void Histo::print(G4int i)
349 void Histo::AddRow() <<
350 { 286 {
351 if (!fNtupleActive) { << 287 G4cout<<"### Histogram "<<i<<" ###"<<G4endl;
352 return; << 288 #ifdef G4ANALYSIS_USE
353 } << 289 if(i>=0 && i<nHisto) {
354 fManager->AddNtupleRow(); << 290 G4double step = (xmax[i] - xmin[i])/G4double( bins[i]);
355 } << 291 G4double x = xmin[i] - step*0.5;
>> 292 G4double y, maxX=0, maxY=0;
>> 293 G4int maxJ=0;
356 294
357 //....oooOO0OOooo........oooOO0OOooo........oo << 295 for(G4int j=0; j<bins[i]; j++) {
>> 296 x += step;
>> 297 y = histo[i]->binHeight(j);
>> 298 if(maxY < y) {maxY = y; maxX = x; maxJ = j;}
358 299
359 void Histo::SetFileName(const G4String& nam) << 300 G4cout<<x<<" "<<y<<G4endl;
360 { << 301 }
361 fHistName = nam; << 302 G4cout<<" maxJ "<<maxJ<<" maxX "<<maxX<<" maxY "<<maxY<<G4endl;
362 fHistoActive = true; << 303 }
>> 304 #endif
363 } 305 }
>> 306 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
364 307
365 //....oooOO0OOooo........oooOO0OOooo........oo <<
366 308